PID Climate (#885)

* PID Climate

* Add sensor for debugging PID output value

* Add dump_config, use percent

* Add more observable values

* Update

* Set target temperature

* Add autotuner

* Add algorithm explanation

* Add autotuner action, update controller

* Add simulator

* Format

* Change defaults

* Updates

esphome/components/pid/climate.py

@@ -0,0 +1,79 @@
+import esphome.codegen as cg
+import esphome.config_validation as cv
+from esphome import automation
+from esphome.components import climate, sensor, output
+from esphome.const import CONF_ID, CONF_SENSOR
+
+pid_ns = cg.esphome_ns.namespace('pid')
+PIDClimate = pid_ns.class_('PIDClimate', climate.Climate, cg.Component)
+PIDAutotuneAction = pid_ns.class_('PIDAutotuneAction', automation.Action)
+
+CONF_DEFAULT_TARGET_TEMPERATURE = 'default_target_temperature'
+
+CONF_KP = 'kp'
+CONF_KI = 'ki'
+CONF_KD = 'kd'
+CONF_CONTROL_PARAMETERS = 'control_parameters'
+CONF_COOL_OUTPUT = 'cool_output'
+CONF_HEAT_OUTPUT = 'heat_output'
+CONF_NOISEBAND = 'noiseband'
+CONF_POSITIVE_OUTPUT = 'positive_output'
+CONF_NEGATIVE_OUTPUT = 'negative_output'
+CONF_MIN_INTEGRAL = 'min_integral'
+CONF_MAX_INTEGRAL = 'max_integral'
+
+CONFIG_SCHEMA = cv.All(climate.CLIMATE_SCHEMA.extend({
+    cv.GenerateID(): cv.declare_id(PIDClimate),
+    cv.Required(CONF_SENSOR): cv.use_id(sensor.Sensor),
+    cv.Required(CONF_DEFAULT_TARGET_TEMPERATURE): cv.temperature,
+    cv.Optional(CONF_COOL_OUTPUT): cv.use_id(output.FloatOutput),
+    cv.Optional(CONF_HEAT_OUTPUT): cv.use_id(output.FloatOutput),
+    cv.Required(CONF_CONTROL_PARAMETERS): cv.Schema({
+        cv.Required(CONF_KP): cv.float_,
+        cv.Optional(CONF_KI, default=0.0): cv.float_,
+        cv.Optional(CONF_KD, default=0.0): cv.float_,
+        cv.Optional(CONF_MIN_INTEGRAL, default=-1): cv.float_,
+        cv.Optional(CONF_MAX_INTEGRAL, default=1): cv.float_,
+    }),
+}), cv.has_at_least_one_key(CONF_COOL_OUTPUT, CONF_HEAT_OUTPUT))
+
+
+def to_code(config):
+    var = cg.new_Pvariable(config[CONF_ID])
+    yield cg.register_component(var, config)
+    yield climate.register_climate(var, config)
+
+    sens = yield cg.get_variable(config[CONF_SENSOR])
+    cg.add(var.set_sensor(sens))
+
+    if CONF_COOL_OUTPUT in config:
+        out = yield cg.get_variable(config[CONF_COOL_OUTPUT])
+        cg.add(var.set_cool_output(out))
+    if CONF_HEAT_OUTPUT in config:
+        out = yield cg.get_variable(config[CONF_HEAT_OUTPUT])
+        cg.add(var.set_heat_output(out))
+    params = config[CONF_CONTROL_PARAMETERS]
+    cg.add(var.set_kp(params[CONF_KP]))
+    cg.add(var.set_ki(params[CONF_KI]))
+    cg.add(var.set_kd(params[CONF_KD]))
+    if CONF_MIN_INTEGRAL in params:
+        cg.add(var.set_min_integral(params[CONF_MIN_INTEGRAL]))
+    if CONF_MAX_INTEGRAL in params:
+        cg.add(var.set_max_integral(params[CONF_MAX_INTEGRAL]))
+
+    cg.add(var.set_default_target_temperature(config[CONF_DEFAULT_TARGET_TEMPERATURE]))
+
+
+@automation.register_action('climate.pid.autotune', PIDAutotuneAction, automation.maybe_simple_id({
+    cv.Required(CONF_ID): cv.use_id(PIDClimate),
+    cv.Optional(CONF_NOISEBAND, default=0.25): cv.float_,
+    cv.Optional(CONF_POSITIVE_OUTPUT, default=1.0): cv.possibly_negative_percentage,
+    cv.Optional(CONF_NEGATIVE_OUTPUT, default=-1.0): cv.possibly_negative_percentage,
+}))
+def esp8266_set_frequency_to_code(config, action_id, template_arg, args):
+    paren = yield cg.get_variable(config[CONF_ID])
+    var = cg.new_Pvariable(action_id, template_arg, paren)
+    cg.add(var.set_noiseband(config[CONF_NOISEBAND]))
+    cg.add(var.set_positive_output(config[CONF_POSITIVE_OUTPUT]))
+    cg.add(var.set_negative_output(config[CONF_NEGATIVE_OUTPUT]))
+    yield var

esphome/components/pid/pid_autotuner.cpp

@@ -0,0 +1,358 @@
+#include "pid_autotuner.h"
+#include "esphome/core/log.h"
+
+namespace esphome {
+namespace pid {
+
+static const char *TAG = "pid.autotune";
+
+/*
+ * # PID Autotuner
+ *
+ * Autotuning of PID parameters is a very interesting topic. There has been
+ * a lot of research over the years to create algorithms that can efficiently determine
+ * suitable starting PID parameters.
+ *
+ * The most basic approach is the Ziegler-Nichols method, which can determine good PID parameters
+ * in a manual process:
+ *  - Set ki, kd to zero.
+ *  - Increase kp until the output oscillates *around* the setpoint. This value kp is called the
+ *    "ultimate gain" K_u.
+ *  - Additionally, record the period of the observed oscillation as P_u (also called T_u).
+ *  - suitable PID parameters are then: kp=0.6*K_u, ki=1.2*K_u/P_u, kd=0.075*K_u*P_u (additional variants of
+ *    these "magic" factors exist as well [2]).
+ *
+ * Now we'd like to automate that process to get K_u and P_u without the user. So we'd like to somehow
+ * make the observed variable oscillate. One observation is that in many applications of PID controllers
+ * the observed variable has some amount of "delay" to the output value (think heating an object, it will
+ * take a few seconds before the sensor can sense the change of temperature) [3].
+ *
+ * It turns out one way to induce such an oscillation is by using a really dumb heating controller:
+ * When the observed value is below the setpoint, heat at 100%. If it's below, cool at 100% (or disable heating).
+ * We call this the "RelayFunction" - the class is responsible for making the observed value oscillate around the
+ * setpoint. We actually use a hysteresis filter (like the bang bang controller) to make the process immune to
+ * noise in the input data, but the math is the same [1].
+ *
+ * Next, now that we have induced an oscillation, we want to measure the frequency (or period) of oscillation.
+ * This is what "OscillationFrequencyDetector" is for: it records zerocrossing events (when the observed value
+ * crosses the setpoint). From that data, we can determine the average oscillating period. This is the P_u of the
+ * ZN-method.
+ *
+ * Finally, we need to determine K_u, the ultimate gain. It turns out we can calculate this based on the amplitude of
+ * oscillation ("induced amplitude `a`) as described in [1]:
+ *   K_u = (4d) / (πa)
+ * where d is the magnitude of the relay function (in range -d to +d).
+ * To measure `a`, we look at the current phase the relay function is in - if it's in the "heating" phase, then we
+ * expect the lowest temperature (=highest error) to be found in the phase because the peak will always happen slightly
+ * after the relay function has changed state (assuming a delay-dominated process).
+ *
+ * Finally, we use some heuristics to determine if the data we've received so far is good:
+ *  - First, of course we must have enough data to calculate the values.
+ *  - The ZC events need to happen at a relatively periodic rate. If the heating/cooling speeds are very different,
+ *    I've observed the ZN parameters are not very useful.
+ *  - The induced amplitude should not deviate too much. If the amplitudes deviate too much this means there has
+ *    been some outside influence (or noise) on the system, and the measured amplitude values are not reliable.
+ *
+ * There are many ways this method can be improved, but on my simulation data the current method already produces very
+ * good results. Some ideas for future improvements:
+ *  - Relay Function improvements:
+ *    - Integrator, Preload, Saturation Relay ([1])
+ *  - Use phase of measured signal relative to relay function.
+ *  - Apply PID parameters from ZN, but continuously tweak them in a second step.
+ *
+ * [1]: https://warwick.ac.uk/fac/cross_fac/iatl/reinvention/archive/volume5issue2/hornsey/
+ * [2]: http://www.mstarlabs.com/control/znrule.html
+ * [3]: https://www.academia.edu/38620114/SEBORG_3rd_Edition_Process_Dynamics_and_Control
+ */
+
+PIDAutotuner::PIDAutotuneResult PIDAutotuner::update(float setpoint, float process_variable) {
+  PIDAutotuner::PIDAutotuneResult res;
+  if (this->state_ == AUTOTUNE_SUCCEEDED) {
+    res.result_params = this->get_ziegler_nichols_pid_();
+    return res;
+  }
+
+  if (!isnan(this->setpoint_) && this->setpoint_ != setpoint) {
+    ESP_LOGW(TAG, "Setpoint changed during autotune! The result will not be accurate!");
+  }
+  this->setpoint_ = setpoint;
+
+  float error = setpoint - process_variable;
+  const uint32_t now = millis();
+
+  float output = this->relay_function_.update(error);
+  this->frequency_detector_.update(now, error);
+  this->amplitude_detector_.update(error, this->relay_function_.state);
+  res.output = output;
+
+  if (!this->frequency_detector_.has_enough_data() || !this->amplitude_detector_.has_enough_data()) {
+    // not enough data for calculation yet
+    ESP_LOGV(TAG, "  Not enough data yet for aututuner");
+    return res;
+  }
+
+  bool zc_symmetrical = this->frequency_detector_.is_increase_decrease_symmetrical();
+  bool amplitude_convergent = this->frequency_detector_.is_increase_decrease_symmetrical();
+  if (!zc_symmetrical || !amplitude_convergent) {
+    // The frequency/amplitude is not fully accurate yet, try to wait
+    // until the fault clears, or terminate after a while anyway
+    if (zc_symmetrical) {
+      ESP_LOGVV(TAG, "  ZC is not symmetrical");
+    }
+    if (amplitude_convergent) {
+      ESP_LOGVV(TAG, "  Amplitude is not convergent");
+    }
+    uint32_t phase = this->relay_function_.phase_count;
+    ESP_LOGVV(TAG, "  Phase %u, enough=%u", phase, enough_data_phase_);
+
+    if (this->enough_data_phase_ == 0) {
+      this->enough_data_phase_ = phase;
+    } else if (phase - this->enough_data_phase_ <= 6) {
+      // keep trying for at least 6 more phases
+      return res;
+    } else {
+      // proceed to calculating PID parameters
+      // warning will be shown in "Checks" section
+    }
+  }
+
+  ESP_LOGI(TAG, "PID Autotune finished!");
+
+  float osc_ampl = this->amplitude_detector_.get_mean_oscillation_amplitude();
+  float d = (this->relay_function_.output_positive - this->relay_function_.output_negative) / 2.0f;
+  ESP_LOGVV(TAG, "  Relay magnitude: %f", d);
+  this->ku_ = 4.0f * d / float(M_PI * osc_ampl);
+  this->pu_ = this->frequency_detector_.get_mean_oscillation_period();
+
+  this->state_ = AUTOTUNE_SUCCEEDED;
+  res.result_params = this->get_ziegler_nichols_pid_();
+  this->dump_config();
+
+  return res;
+}
+void PIDAutotuner::dump_config() {
+  ESP_LOGI(TAG, "PID Autotune:");
+  if (this->state_ == AUTOTUNE_SUCCEEDED) {
+    ESP_LOGI(TAG, "  State: Succeeded!");
+    bool has_issue = false;
+    if (!this->amplitude_detector_.is_amplitude_convergent()) {
+      ESP_LOGW(TAG, "  Could not reliable determine oscillation amplitude, PID parameters may be inaccurate!");
+      ESP_LOGW(TAG, "    Please make sure you eliminate all outside influences on the measured temperature.");
+      has_issue = true;
+    }
+    if (!this->frequency_detector_.is_increase_decrease_symmetrical()) {
+      ESP_LOGW(TAG, "  Oscillation Frequency is not symmetrical. PID parameters may be inaccurate!");
+      ESP_LOGW(
+          TAG,
+          "    This is usually because the heat and cool processes do not change the temperature at the same rate.");
+      ESP_LOGW(TAG,
+               "    Please try reducing the positive_output value (or increase negative_output in case of a cooler)");
+      has_issue = true;
+    }
+    if (!has_issue) {
+      ESP_LOGI(TAG, "  All checks passed!");
+    }
+
+    auto fac = get_ziegler_nichols_pid_();
+    ESP_LOGI(TAG, "  Calculated PID parameters (\"Ziegler-Nichols PID\" rule):");
+    ESP_LOGI(TAG, " ");
+    ESP_LOGI(TAG, "  control_parameters:");
+    ESP_LOGI(TAG, "    kp: %.5f", fac.kp);
+    ESP_LOGI(TAG, "    ki: %.5f", fac.ki);
+    ESP_LOGI(TAG, "    kd: %.5f", fac.kd);
+    ESP_LOGI(TAG, " ");
+    ESP_LOGI(TAG, "  Please copy these values into your YAML configuration! They will reset on the next reboot.");
+
+    ESP_LOGV(TAG, "  Oscillation Period: %f", this->frequency_detector_.get_mean_oscillation_period());
+    ESP_LOGV(TAG, "  Oscillation Amplitude: %f", this->amplitude_detector_.get_mean_oscillation_amplitude());
+    ESP_LOGV(TAG, "  Ku: %f, Pu: %f", this->ku_, this->pu_);
+
+    ESP_LOGD(TAG, "  Alternative Rules:");
+    // http://www.mstarlabs.com/control/znrule.html
+    print_rule_("Ziegler-Nichols PI", 0.45f, 0.54f, 0.0f);
+    print_rule_("Pessen Integral PID", 0.7f, 1.75f, 0.105f);
+    print_rule_("Some Overshoot PID", 0.333f, 0.667f, 0.111f);
+    print_rule_("No Overshoot PID", 0.2f, 0.4f, 0.0625f);
+  }
+
+  if (this->state_ == AUTOTUNE_RUNNING) {
+    ESP_LOGI(TAG, "  Autotune is still running!");
+    ESP_LOGD(TAG, "  Status: Trying to reach %.2f °C", setpoint_ - relay_function_.current_target_error());
+    ESP_LOGD(TAG, "  Stats so far:");
+    ESP_LOGD(TAG, "    Phases: %u", relay_function_.phase_count);
+    ESP_LOGD(TAG, "    Detected %u zero-crossings", frequency_detector_.zerocrossing_intervals.size());  // NOLINT
+    ESP_LOGD(TAG, "    Current Phase Min: %.2f, Max: %.2f", amplitude_detector_.phase_min,
+             amplitude_detector_.phase_max);
+  }
+}
+PIDAutotuner::PIDResult PIDAutotuner::calculate_pid_(float kp_factor, float ki_factor, float kd_factor) {
+  float kp = kp_factor * ku_;
+  float ki = ki_factor * ku_ / pu_;
+  float kd = kd_factor * ku_ * pu_;
+  return {
+      .kp = kp,
+      .ki = ki,
+      .kd = kd,
+  };
+}
+void PIDAutotuner::print_rule_(const char *name, float kp_factor, float ki_factor, float kd_factor) {
+  auto fac = calculate_pid_(kp_factor, ki_factor, kd_factor);
+  ESP_LOGD(TAG, "    Rule '%s':", name);
+  ESP_LOGD(TAG, "      kp: %.5f, ki: %.5f, kd: %.5f", fac.kp, fac.ki, fac.kd);
+}
+
+// ================== RelayFunction ==================
+float PIDAutotuner::RelayFunction::update(float error) {
+  if (this->state == RELAY_FUNCTION_INIT) {
+    bool pos = error > this->noiseband;
+    state = pos ? RELAY_FUNCTION_POSITIVE : RELAY_FUNCTION_NEGATIVE;
+  }
+  bool change = false;
+  if (this->state == RELAY_FUNCTION_POSITIVE && error < -this->noiseband) {
+    // Positive hysteresis reached, change direction
+    this->state = RELAY_FUNCTION_NEGATIVE;
+    change = true;
+  } else if (this->state == RELAY_FUNCTION_NEGATIVE && error > this->noiseband) {
+    // Negative hysteresis reached, change direction
+    this->state = RELAY_FUNCTION_POSITIVE;
+    change = true;
+  }
+
+  float output = state == RELAY_FUNCTION_POSITIVE ? output_positive : output_negative;
+  if (change) {
+    this->phase_count++;
+    ESP_LOGV(TAG, "Autotune: Turning output to %.1f%%", output * 100);
+  }
+
+  return output;
+}
+
+// ================== OscillationFrequencyDetector ==================
+void PIDAutotuner::OscillationFrequencyDetector::update(uint32_t now, float error) {
+  if (this->state == FREQUENCY_DETECTOR_INIT) {
+    bool pos = error > this->noiseband;
+    state = pos ? FREQUENCY_DETECTOR_POSITIVE : FREQUENCY_DETECTOR_NEGATIVE;
+  }
+
+  bool had_crossing = false;
+  if (this->state == FREQUENCY_DETECTOR_POSITIVE && error < -this->noiseband) {
+    this->state = FREQUENCY_DETECTOR_NEGATIVE;
+    had_crossing = true;
+  } else if (this->state == FREQUENCY_DETECTOR_NEGATIVE && error > this->noiseband) {
+    this->state = FREQUENCY_DETECTOR_POSITIVE;
+    had_crossing = true;
+  }
+
+  if (had_crossing) {
+    // Had crossing above hysteresis threshold, record
+    ESP_LOGV(TAG, "Autotune: Detected Zero-Cross at %u", now);
+    if (this->last_zerocross != 0) {
+      uint32_t dt = now - this->last_zerocross;
+      ESP_LOGV(TAG, "  dt: %u", dt);
+      this->zerocrossing_intervals.push_back(dt);
+    }
+    this->last_zerocross = now;
+  }
+}
+bool PIDAutotuner::OscillationFrequencyDetector::has_enough_data() const {
+  // Do we have enough data in this detector to generate PID values?
+  return this->zerocrossing_intervals.size() >= 2;
+}
+float PIDAutotuner::OscillationFrequencyDetector::get_mean_oscillation_period() const {
+  // Get the mean oscillation period in seconds
+  // Only call if has_enough_data() has returned true.
+  float sum = 0.0f;
+  for (uint32_t v : this->zerocrossing_intervals)
+    sum += v;
+  // zerocrossings are each half-period, multiply by 2
+  float mean_value = sum / this->zerocrossing_intervals.size();
+  // divide by 1000 to get seconds, multiply by two because zc happens two times per period
+  float mean_period = mean_value / 1000 * 2;
+  return mean_period;
+}
+bool PIDAutotuner::OscillationFrequencyDetector::is_increase_decrease_symmetrical() const {
+  // Check if increase/decrease of process value was symmetrical
+  // If the process value increases much faster than it decreases, the generated PID values will
+  // not be very good and the function output values need to be adjusted
+  // Happens for example with a well-insulated heating element.
+  // We calculate this based on the zerocrossing interval.
+  if (zerocrossing_intervals.empty())
+    return false;
+  uint32_t max_interval = zerocrossing_intervals[0];
+  uint32_t min_interval = zerocrossing_intervals[0];
+  for (uint32_t interval : zerocrossing_intervals) {
+    max_interval = std::max(max_interval, interval);
+    min_interval = std::min(min_interval, interval);
+  }
+  float ratio = min_interval / float(max_interval);
+  return ratio >= 0.66;
+}
+
+// ================== OscillationAmplitudeDetector ==================
+void PIDAutotuner::OscillationAmplitudeDetector::update(float error,
+                                                        PIDAutotuner::RelayFunction::RelayFunctionState relay_state) {
+  if (relay_state != last_relay_state) {
+    if (last_relay_state == RelayFunction::RELAY_FUNCTION_POSITIVE) {
+      // Transitioned from positive error to negative error.
+      // The positive error peak must have been in previous segment (180° shifted)
+      // record phase_max
+      this->phase_maxs.push_back(phase_max);
+      ESP_LOGV(TAG, "Autotune: Phase Max: %f", phase_max);
+    } else if (last_relay_state == RelayFunction::RELAY_FUNCTION_NEGATIVE) {
+      // Transitioned from negative error to positive error.
+      // The negative error peak must have been in previous segment (180° shifted)
+      // record phase_min
+      this->phase_mins.push_back(phase_min);
+      ESP_LOGV(TAG, "Autotune: Phase Min: %f", phase_min);
+    }
+    // reset phase values for next phase
+    this->phase_min = error;
+    this->phase_max = error;
+  }
+  this->last_relay_state = relay_state;
+
+  this->phase_min = std::min(this->phase_min, error);
+  this->phase_max = std::max(this->phase_max, error);
+
+  // Check arrays sizes, we keep at most 7 items (6 datapoints is enough, and data at beginning might not
+  // have been stabilized)
+  if (this->phase_maxs.size() > 7)
+    this->phase_maxs.erase(this->phase_maxs.begin());
+  if (this->phase_mins.size() > 7)
+    this->phase_mins.erase(this->phase_mins.begin());
+}
+bool PIDAutotuner::OscillationAmplitudeDetector::has_enough_data() const {
+  // Return if we have enough data to generate PID parameters
+  // The first phase is not very useful if the setpoint is not set to the starting process value
+  // So discard first phase. Otherwise we need at least two phases.
+  return std::min(phase_mins.size(), phase_maxs.size()) >= 3;
+}
+float PIDAutotuner::OscillationAmplitudeDetector::get_mean_oscillation_amplitude() const {
+  float total_amplitudes = 0;
+  size_t total_amplitudes_n = 0;
+  for (int i = 1; i < std::min(phase_mins.size(), phase_maxs.size()) - 1; i++) {
+    total_amplitudes += std::abs(phase_maxs[i] - phase_mins[i + 1]);
+    total_amplitudes_n++;
+  }
+  float mean_amplitude = total_amplitudes / total_amplitudes_n;
+  // Amplitude is measured from center, divide by 2
+  return mean_amplitude / 2.0f;
+}
+bool PIDAutotuner::OscillationAmplitudeDetector::is_amplitude_convergent() const {
+  // Check if oscillation amplitude is convergent
+  // We implement this by checking global extrema against average amplitude
+  if (this->phase_mins.empty() || this->phase_maxs.empty())
+    return false;
+
+  float global_max = phase_maxs[0], global_min = phase_mins[0];
+  for (auto v : this->phase_mins)
+    global_min = std::min(global_min, v);
+  for (auto v : this->phase_maxs)
+    global_max = std::min(global_max, v);
+  float global_amplitude = (global_max - global_min) / 2.0f;
+  float mean_amplitude = this->get_mean_oscillation_amplitude();
+  return (mean_amplitude - global_amplitude) / (global_amplitude) < 0.05f;
+}
+
+}  // namespace pid
+}  // namespace esphome

esphome/components/pid/pid_autotuner.h

@@ -0,0 +1,110 @@
+#pragma once
+
+#include "esphome/core/component.h"
+#include "esphome/core/optional.h"
+#include "pid_controller.h"
+#include "pid_simulator.h"
+
+namespace esphome {
+namespace pid {
+
+class PIDAutotuner {
+ public:
+  struct PIDResult {
+    float kp;
+    float ki;
+    float kd;
+  };
+  struct PIDAutotuneResult {
+    float output;
+    optional<PIDResult> result_params;
+  };
+
+  void config(float output_min, float output_max) {
+    relay_function_.output_negative = std::max(relay_function_.output_negative, output_min);
+    relay_function_.output_positive = std::min(relay_function_.output_positive, output_max);
+  }
+  PIDAutotuneResult update(float setpoint, float process_variable);
+  bool is_finished() const { return state_ != AUTOTUNE_RUNNING; }
+
+  void dump_config();
+
+  void set_noiseband(float noiseband) {
+    relay_function_.noiseband = noiseband;
+    // ZC detector uses 1/4 the noiseband of relay function (noise suppression)
+    frequency_detector_.noiseband = noiseband / 4;
+  }
+  void set_output_positive(float output_positive) { relay_function_.output_positive = output_positive; }
+  void set_output_negative(float output_negative) { relay_function_.output_negative = output_negative; }
+
+ protected:
+  struct RelayFunction {
+    float update(float error);
+
+    float current_target_error() const {
+      if (state == RELAY_FUNCTION_INIT)
+        return 0;
+      if (state == RELAY_FUNCTION_POSITIVE)
+        return -noiseband;
+      return noiseband;
+    }
+
+    enum RelayFunctionState {
+      RELAY_FUNCTION_INIT,
+      RELAY_FUNCTION_POSITIVE,
+      RELAY_FUNCTION_NEGATIVE,
+    } state = RELAY_FUNCTION_INIT;
+    float noiseband = 0.5;
+    float output_positive = 1;
+    float output_negative = -1;
+    uint32_t phase_count = 0;
+  } relay_function_;
+  struct OscillationFrequencyDetector {
+    void update(uint32_t now, float error);
+
+    bool has_enough_data() const;
+
+    float get_mean_oscillation_period() const;
+
+    bool is_increase_decrease_symmetrical() const;
+
+    enum FrequencyDetectorState {
+      FREQUENCY_DETECTOR_INIT,
+      FREQUENCY_DETECTOR_POSITIVE,
+      FREQUENCY_DETECTOR_NEGATIVE,
+    } state;
+    float noiseband = 0.05;
+    uint32_t last_zerocross{0};
+    std::vector<uint32_t> zerocrossing_intervals;
+  } frequency_detector_;
+  struct OscillationAmplitudeDetector {
+    void update(float error, RelayFunction::RelayFunctionState relay_state);
+
+    bool has_enough_data() const;
+
+    float get_mean_oscillation_amplitude() const;
+
+    bool is_amplitude_convergent() const;
+
+    float phase_min = NAN;
+    float phase_max = NAN;
+    std::vector<float> phase_mins;
+    std::vector<float> phase_maxs;
+    RelayFunction::RelayFunctionState last_relay_state = RelayFunction::RELAY_FUNCTION_INIT;
+  } amplitude_detector_;
+  PIDResult calculate_pid_(float kp_factor, float ki_factor, float kd_factor);
+  void print_rule_(const char *name, float kp_factor, float ki_factor, float kd_factor);
+  PIDResult get_ziegler_nichols_pid_() { return calculate_pid_(0.6f, 1.2f, 0.075f); }
+
+  uint32_t enough_data_phase_ = 0;
+  float setpoint_ = NAN;
+  enum State {
+    AUTOTUNE_RUNNING,
+    AUTOTUNE_SUCCEEDED,
+  } state_ = AUTOTUNE_RUNNING;
+  float ku_;
+  float pu_;
+};
+
+}  // namespace pid
+}  // namespace esphome

esphome/components/pid/pid_climate.cpp

@@ -0,0 +1,152 @@
+#include "pid_climate.h"
+#include "esphome/core/log.h"
+
+namespace esphome {
+namespace pid {
+
+static const char *TAG = "pid.climate";
+
+void PIDClimate::setup() {
+  this->sensor_->add_on_state_callback([this](float state) {
+    // only publish if state/current temperature has changed in two digits of precision
+    this->do_publish_ = roundf(state * 100) != roundf(this->current_temperature * 100);
+    this->current_temperature = state;
+    this->update_pid_();
+  });
+  this->current_temperature = this->sensor_->state;
+  // restore set points
+  auto restore = this->restore_state_();
+  if (restore.has_value()) {
+    restore->to_call(this).perform();
+  } else {
+    // restore from defaults, change_away handles those for us
+    this->mode = climate::CLIMATE_MODE_AUTO;
+    this->target_temperature = this->default_target_temperature_;
+  }
+}
+void PIDClimate::control(const climate::ClimateCall &call) {
+  if (call.get_mode().has_value())
+    this->mode = *call.get_mode();
+  if (call.get_target_temperature().has_value())
+    this->target_temperature = *call.get_target_temperature();
+
+  // If switching to non-auto mode, set output immediately
+  if (this->mode != climate::CLIMATE_MODE_AUTO)
+    this->handle_non_auto_mode_();
+
+  this->publish_state();
+}
+climate::ClimateTraits PIDClimate::traits() {
+  auto traits = climate::ClimateTraits();
+  traits.set_supports_current_temperature(true);
+  traits.set_supports_auto_mode(true);
+  traits.set_supports_two_point_target_temperature(false);
+  traits.set_supports_cool_mode(this->supports_cool_());
+  traits.set_supports_heat_mode(this->supports_heat_());
+  traits.set_supports_action(true);
+  return traits;
+}
+void PIDClimate::dump_config() {
+  LOG_CLIMATE("", "PID Climate", this);
+  ESP_LOGCONFIG(TAG, "  Control Parameters:");
+  ESP_LOGCONFIG(TAG, "    kp: %.5f, ki: %.5f, kd: %.5f", controller_.kp, controller_.ki, controller_.kd);
+
+  if (this->autotuner_ != nullptr) {
+    this->autotuner_->dump_config();
+  }
+}
+void PIDClimate::write_output_(float value) {
+  this->output_value_ = value;
+
+  // first ensure outputs are off (both outputs not active at the same time)
+  if (this->supports_cool_() && value >= 0)
+    this->cool_output_->set_level(0.0f);
+  if (this->supports_heat_() && value <= 0)
+    this->heat_output_->set_level(0.0f);
+
+  // value < 0 means cool, > 0 means heat
+  if (this->supports_cool_() && value < 0)
+    this->cool_output_->set_level(std::min(1.0f, -value));
+  if (this->supports_heat_() && value > 0)
+    this->heat_output_->set_level(std::min(1.0f, value));
+
+  // Update action variable for user feedback what's happening
+  climate::ClimateAction new_action;
+  if (this->supports_cool_() && value < 0)
+    new_action = climate::CLIMATE_ACTION_COOLING;
+  else if (this->supports_heat_() && value > 0)
+    new_action = climate::CLIMATE_ACTION_HEATING;
+  else if (this->mode == climate::CLIMATE_MODE_OFF)
+    new_action = climate::CLIMATE_ACTION_OFF;
+  else
+    new_action = climate::CLIMATE_ACTION_IDLE;
+
+  if (new_action != this->action) {
+    this->action = new_action;
+    this->do_publish_ = true;
+  }
+  this->pid_computed_callback_.call();
+}
+void PIDClimate::handle_non_auto_mode_() {
+  // in non-auto mode, switch directly to appropriate action
+  //  - HEAT mode / COOL mode -> Output at ±100%
+  //  - OFF mode -> Output at 0%
+  if (this->mode == climate::CLIMATE_MODE_HEAT) {
+    this->write_output_(1.0);
+  } else if (this->mode == climate::CLIMATE_MODE_COOL) {
+    this->write_output_(-1.0);
+  } else if (this->mode == climate::CLIMATE_MODE_OFF) {
+    this->write_output_(0.0);
+  } else {
+    assert(false);
+  }
+}
+void PIDClimate::update_pid_() {
+  float value;
+  if (isnan(this->current_temperature) || isnan(this->target_temperature)) {
+    // if any control parameters are nan, turn off all outputs
+    value = 0.0;
+  } else {
+    // Update PID controller irrespective of current mode, to not mess up D/I terms
+    // In non-auto mode, we just discard the output value
+    value = this->controller_.update(this->target_temperature, this->current_temperature);
+
+    // Check autotuner
+    if (this->autotuner_ != nullptr && !this->autotuner_->is_finished()) {
+      auto res = this->autotuner_->update(this->target_temperature, this->current_temperature);
+      if (res.result_params.has_value()) {
+        this->controller_.kp = res.result_params->kp;
+        this->controller_.ki = res.result_params->ki;
+        this->controller_.kd = res.result_params->kd;
+        // keep autotuner instance so that subsequent dump_configs will print the long result message.
+      } else {
+        value = res.output;
+        if (mode != climate::CLIMATE_MODE_AUTO) {
+          ESP_LOGW(TAG, "For PID autotuner you need to set AUTO (also called heat/cool) mode!");
+        }
+      }
+    }
+  }
+
+  if (this->mode != climate::CLIMATE_MODE_AUTO) {
+    this->handle_non_auto_mode_();
+  } else {
+    this->write_output_(value);
+  }
+
+  if (this->do_publish_)
+    this->publish_state();
+}
+void PIDClimate::start_autotune(std::unique_ptr<PIDAutotuner> &&autotune) {
+  this->autotuner_ = std::move(autotune);
+  float min_value = this->supports_cool_() ? -1.0f : 0.0f;
+  float max_value = this->supports_heat_() ? 1.0f : 0.0f;
+  this->autotuner_->config(min_value, max_value);
+  this->set_interval("autotune-progress", 10000, [this]() {
+    if (this->autotuner_ != nullptr && !this->autotuner_->is_finished())
+      this->autotuner_->dump_config();
+  });
+}
+
+}  // namespace pid
+}  // namespace esphome

esphome/components/pid/pid_climate.h

@@ -0,0 +1,94 @@
+#pragma once
+
+#include "esphome/core/component.h"
+#include "esphome/core/helpers.h"
+#include "esphome/core/automation.h"
+#include "esphome/components/climate/climate.h"
+#include "esphome/components/sensor/sensor.h"
+#include "esphome/components/output/float_output.h"
+#include "pid_controller.h"
+#include "pid_autotuner.h"
+
+namespace esphome {
+namespace pid {
+
+class PIDClimate : public climate::Climate, public Component {
+ public:
+  PIDClimate() = default;
+  void setup() override;
+  void dump_config() override;
+
+  void set_sensor(sensor::Sensor *sensor) { sensor_ = sensor; }
+  void set_cool_output(output::FloatOutput *cool_output) { cool_output_ = cool_output; }
+  void set_heat_output(output::FloatOutput *heat_output) { heat_output_ = heat_output; }
+  void set_kp(float kp) { controller_.kp = kp; }
+  void set_ki(float ki) { controller_.ki = ki; }
+  void set_kd(float kd) { controller_.kd = kd; }
+  void set_min_integral(float min_integral) { controller_.min_integral = min_integral; }
+  void set_max_integral(float max_integral) { controller_.max_integral = max_integral; }
+
+  float get_output_value() const { return output_value_; }
+  float get_error_value() const { return controller_.error; }
+  float get_proportional_term() const { return controller_.proportional_term; }
+  float get_integral_term() const { return controller_.integral_term; }
+  float get_derivative_term() const { return controller_.derivative_term; }
+  void add_on_pid_computed_callback(std::function<void()> &&callback) {
+    pid_computed_callback_.add(std::move(callback));
+  }
+  void set_default_target_temperature(float default_target_temperature) {
+    default_target_temperature_ = default_target_temperature;
+  }
+  void start_autotune(std::unique_ptr<PIDAutotuner> &&autotune);
+
+ protected:
+  /// Override control to change settings of the climate device.
+  void control(const climate::ClimateCall &call) override;
+  /// Return the traits of this controller.
+  climate::ClimateTraits traits() override;
+
+  void update_pid_();
+
+  bool supports_cool_() const { return this->cool_output_ != nullptr; }
+  bool supports_heat_() const { return this->heat_output_ != nullptr; }
+
+  void write_output_(float value);
+  void handle_non_auto_mode_();
+
+  /// The sensor used for getting the current temperature
+  sensor::Sensor *sensor_;
+  output::FloatOutput *cool_output_ = nullptr;
+  output::FloatOutput *heat_output_ = nullptr;
+  PIDController controller_;
+  /// Output value as reported by the PID controller, for PIDClimateSensor
+  float output_value_;
+  CallbackManager<void()> pid_computed_callback_;
+  float default_target_temperature_;
+  std::unique_ptr<PIDAutotuner> autotuner_;
+  bool do_publish_ = false;
+};
+
+template<typename... Ts> class PIDAutotuneAction : public Action<Ts...> {
+ public:
+  PIDAutotuneAction(PIDClimate *parent) : parent_(parent) {}
+
+  void play(Ts... x) {
+    auto tuner = make_unique<PIDAutotuner>();
+    tuner->set_noiseband(this->noiseband_);
+    tuner->set_output_negative(this->negative_output_);
+    tuner->set_output_positive(this->positive_output_);
+    this->parent_->start_autotune(std::move(tuner));
+  }
+
+  void set_noiseband(float noiseband) { noiseband_ = noiseband; }
+  void set_positive_output(float positive_output) { positive_output_ = positive_output; }
+  void set_negative_output(float negative_output) { negative_output_ = negative_output; }
+
+ protected:
+  float noiseband_;
+  float positive_output_;
+  float negative_output_;
+  PIDClimate *parent_;
+};
+
+}  // namespace pid
+}  // namespace esphome

esphome/components/pid/pid_controller.h

@@ -0,0 +1,79 @@
+#pragma once
+
+#include "esphome/core/esphal.h"
+
+namespace esphome {
+namespace pid {
+
+struct PIDController {
+  float update(float setpoint, float process_value) {
+    // e(t) ... error at timestamp t
+    // r(t) ... setpoint
+    // y(t) ... process value (sensor reading)
+    // u(t) ... output value
+
+    float dt = calculate_relative_time_();
+
+    // e(t) := r(t) - y(t)
+    error = setpoint - process_value;
+
+    // p(t) := K_p * e(t)
+    proportional_term = kp * error;
+
+    // i(t) := K_i * \int_{0}^{t} e(t) dt
+    accumulated_integral_ += error * dt * ki;
+    // constrain accumulated integral value
+    if (!isnan(min_integral) && accumulated_integral_ < min_integral)
+      accumulated_integral_ = min_integral;
+    if (!isnan(max_integral) && accumulated_integral_ > max_integral)
+      accumulated_integral_ = max_integral;
+    integral_term = accumulated_integral_;
+
+    // d(t) := K_d * de(t)/dt
+    float derivative = 0.0f;
+    if (dt != 0.0f)
+      derivative = (error - previous_error_) / dt;
+    previous_error_ = error;
+    derivative_term = kd * derivative;
+
+    // u(t) := p(t) + i(t) + d(t)
+    return proportional_term + integral_term + derivative_term;
+  }
+
+  /// Proportional gain K_p.
+  float kp = 0;
+  /// Integral gain K_i.
+  float ki = 0;
+  /// Differential gain K_d.
+  float kd = 0;
+
+  float min_integral = NAN;
+  float max_integral = NAN;
+
+  // Store computed values in struct so that values can be monitored through sensors
+  float error;
+  float proportional_term;
+  float integral_term;
+  float derivative_term;
+
+ protected:
+  float calculate_relative_time_() {
+    uint32_t now = millis();
+    uint32_t dt = now - this->last_time_;
+    if (last_time_ == 0) {
+      last_time_ = now;
+      return 0.0f;
+    }
+    last_time_ = now;
+    return dt / 1000.0f;
+  }
+
+  /// Error from previous update used for derivative term
+  float previous_error_ = 0;
+  /// Accumulated integral value
+  float accumulated_integral_ = 0;
+  uint32_t last_time_ = 0;
+};
+
+}  // namespace pid
+}  // namespace esphome

esphome/components/pid/pid_simulator.h

@@ -0,0 +1,75 @@
+#pragma once
+
+#include "esphome/core/component.h"
+#include "esphome/core/helpers.h"
+#include "esphome/components/sensor/sensor.h"
+#include "esphome/components/output/float_output.h"
+
+namespace esphome {
+namespace pid {
+
+class PIDSimulator : public PollingComponent, public output::FloatOutput {
+ public:
+  PIDSimulator() : PollingComponent(1000) {}
+
+  float surface = 1;                     /// surface area in m²
+  float mass = 3;                        /// mass of simulated object in kg
+  float temperature = 21;                /// current temperature of object in °C
+  float efficiency = 0.98;               /// heating efficiency, 1 is 100% efficient
+  float thermal_conductivity = 15;       /// thermal conductivity of surface are in W/(m*K), here: steel
+  float specific_heat_capacity = 4.182;  /// specific heat capacity of mass in kJ/(kg*K), here: water
+  float heat_power = 500;                /// Heating power in W
+  float ambient_temperature = 20;        /// Ambient temperature in °C
+  float update_interval = 1;             /// The simulated updated interval in seconds
+  std::vector<float> delayed_temps;      /// storage of past temperatures for delaying temperature reading
+  size_t delay_cycles = 15;              /// how many update cycles to delay the output
+  float output_value = 0.0;              /// Current output value of heating element
+  sensor::Sensor *sensor = new sensor::Sensor();
+
+  float delta_t(float power) {
+    // P = Q / t
+    // Q = c * m * 𝚫t
+    // 𝚫t = (P*t) / (c*m)
+    float c = this->specific_heat_capacity;
+    float t = this->update_interval;
+    float p = power / 1000;  //  in kW
+    float m = this->mass;
+    return (p * t) / (c * m);
+  }
+
+  float update_temp() {
+    float value = clamp(output_value, 0.0f, 1.0f);
+
+    // Heat
+    float power = value * heat_power * efficiency;
+    temperature += this->delta_t(power);
+
+    // Cool
+    // Q = k_w * A * (T_mass - T_ambient)
+    // P = Q / t
+    float dt = temperature - ambient_temperature;
+    float cool_power = (thermal_conductivity * surface * dt) / update_interval;
+    temperature -= this->delta_t(cool_power);
+
+    // Delay temperature readings
+    delayed_temps.push_back(temperature);
+    if (delayed_temps.size() > delay_cycles)
+      delayed_temps.erase(delayed_temps.begin());
+    float prev_temp = this->delayed_temps[0];
+    float alpha = 0.1f;
+    float ret = (1 - alpha) * prev_temp + alpha * prev_temp;
+    return ret;
+  }
+
+  void setup() override { sensor->publish_state(this->temperature); }
+  void update() override {
+    float new_temp = this->update_temp();
+    sensor->publish_state(new_temp);
+  }
+
+ protected:
+  void write_state(float state) override { this->output_value = state; }
+};
+
+}  // namespace pid
+}  // namespace esphome

esphome/components/pid/sensor/__init__.py

@@ -0,0 +1,36 @@
+import esphome.codegen as cg
+import esphome.config_validation as cv
+from esphome.components import sensor
+from esphome.const import CONF_ID, UNIT_PERCENT, ICON_GAUGE, CONF_TYPE
+from ..climate import pid_ns, PIDClimate
+
+PIDClimateSensor = pid_ns.class_('PIDClimateSensor', sensor.Sensor, cg.Component)
+PIDClimateSensorType = pid_ns.enum('PIDClimateSensorType')
+
+PID_CLIMATE_SENSOR_TYPES = {
+    'RESULT': PIDClimateSensorType.PID_SENSOR_TYPE_RESULT,
+    'ERROR': PIDClimateSensorType.PID_SENSOR_TYPE_ERROR,
+    'PROPORTIONAL': PIDClimateSensorType.PID_SENSOR_TYPE_PROPORTIONAL,
+    'INTEGRAL': PIDClimateSensorType.PID_SENSOR_TYPE_INTEGRAL,
+    'DERIVATIVE': PIDClimateSensorType.PID_SENSOR_TYPE_DERIVATIVE,
+    'HEAT': PIDClimateSensorType.PID_SENSOR_TYPE_HEAT,
+    'COOL': PIDClimateSensorType.PID_SENSOR_TYPE_COOL,
+}
+
+CONF_CLIMATE_ID = 'climate_id'
+CONFIG_SCHEMA = sensor.sensor_schema(UNIT_PERCENT, ICON_GAUGE, 1).extend({
+    cv.GenerateID(): cv.declare_id(PIDClimateSensor),
+    cv.GenerateID(CONF_CLIMATE_ID): cv.use_id(PIDClimate),
+
+    cv.Required(CONF_TYPE): cv.enum(PID_CLIMATE_SENSOR_TYPES, upper=True),
+}).extend(cv.COMPONENT_SCHEMA)
+
+
+def to_code(config):
+    parent = yield cg.get_variable(config[CONF_CLIMATE_ID])
+    var = cg.new_Pvariable(config[CONF_ID])
+    yield sensor.register_sensor(var, config)
+    yield cg.register_component(var, config)
+
+    cg.add(var.set_parent(parent))
+    cg.add(var.set_type(config[CONF_TYPE]))

esphome/components/pid/sensor/pid_climate_sensor.cpp

@@ -0,0 +1,47 @@
+#include "pid_climate_sensor.h"
+#include "esphome/core/log.h"
+#include "esphome/core/helpers.h"
+
+namespace esphome {
+namespace pid {
+
+static const char *TAG = "pid.sensor";
+
+void PIDClimateSensor::setup() {
+  this->parent_->add_on_pid_computed_callback([this]() { this->update_from_parent_(); });
+  this->update_from_parent_();
+}
+void PIDClimateSensor::update_from_parent_() {
+  float value;
+  switch (this->type_) {
+    case PID_SENSOR_TYPE_RESULT:
+      value = this->parent_->get_output_value();
+      break;
+    case PID_SENSOR_TYPE_ERROR:
+      value = this->parent_->get_error_value();
+      break;
+    case PID_SENSOR_TYPE_PROPORTIONAL:
+      value = this->parent_->get_proportional_term();
+      break;
+    case PID_SENSOR_TYPE_INTEGRAL:
+      value = this->parent_->get_integral_term();
+      break;
+    case PID_SENSOR_TYPE_DERIVATIVE:
+      value = this->parent_->get_derivative_term();
+      break;
+    case PID_SENSOR_TYPE_HEAT:
+      value = clamp(this->parent_->get_output_value(), 0.0f, 1.0f);
+      break;
+    case PID_SENSOR_TYPE_COOL:
+      value = clamp(-this->parent_->get_output_value(), 0.0f, 1.0f);
+      break;
+    default:
+      value = NAN;
+      break;
+  }
+  this->publish_state(value * 100.0f);
+}
+void PIDClimateSensor::dump_config() { LOG_SENSOR("", "PID Climate Sensor", this); }
+
+}  // namespace pid
+}  // namespace esphome

esphome/components/pid/sensor/pid_climate_sensor.h

@@ -0,0 +1,34 @@
+#pragma once
+
+#include "esphome/core/component.h"
+#include "esphome/components/pid/pid_climate.h"
+
+namespace esphome {
+namespace pid {
+
+enum PIDClimateSensorType {
+  PID_SENSOR_TYPE_RESULT,
+  PID_SENSOR_TYPE_ERROR,
+  PID_SENSOR_TYPE_PROPORTIONAL,
+  PID_SENSOR_TYPE_INTEGRAL,
+  PID_SENSOR_TYPE_DERIVATIVE,
+  PID_SENSOR_TYPE_HEAT,
+  PID_SENSOR_TYPE_COOL,
+};
+
+class PIDClimateSensor : public sensor::Sensor, public Component {
+ public:
+  void setup() override;
+  void set_parent(PIDClimate *parent) { parent_ = parent; }
+  void set_type(PIDClimateSensorType type) { type_ = type; }
+
+  void dump_config() override;
+
+ protected:
+  void update_from_parent_();
+  PIDClimate *parent_;
+  PIDClimateSensorType type_;
+};
+
+}  // namespace pid
+}  // namespace esphome