ESPHome  2022.8.0
pid_autotuner.cpp
Go to the documentation of this file.
1 #include "pid_autotuner.h"
2 #include "esphome/core/log.h"
3 
4 namespace esphome {
5 namespace pid {
6 
7 static const char *const TAG = "pid.autotune";
8 
9 /*
10  * # PID Autotuner
11  *
12  * Autotuning of PID parameters is a very interesting topic. There has been
13  * a lot of research over the years to create algorithms that can efficiently determine
14  * suitable starting PID parameters.
15  *
16  * The most basic approach is the Ziegler-Nichols method, which can determine good PID parameters
17  * in a manual process:
18  * - Set ki, kd to zero.
19  * - Increase kp until the output oscillates *around* the setpoint. This value kp is called the
20  * "ultimate gain" K_u.
21  * - Additionally, record the period of the observed oscillation as P_u (also called T_u).
22  * - 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
23  * these "magic" factors exist as well [2]).
24  *
25  * Now we'd like to automate that process to get K_u and P_u without the user. So we'd like to somehow
26  * make the observed variable oscillate. One observation is that in many applications of PID controllers
27  * the observed variable has some amount of "delay" to the output value (think heating an object, it will
28  * take a few seconds before the sensor can sense the change of temperature) [3].
29  *
30  * It turns out one way to induce such an oscillation is by using a really dumb heating controller:
31  * When the observed value is below the setpoint, heat at 100%. If it's below, cool at 100% (or disable heating).
32  * We call this the "RelayFunction" - the class is responsible for making the observed value oscillate around the
33  * setpoint. We actually use a hysteresis filter (like the bang bang controller) to make the process immune to
34  * noise in the input data, but the math is the same [1].
35  *
36  * Next, now that we have induced an oscillation, we want to measure the frequency (or period) of oscillation.
37  * This is what "OscillationFrequencyDetector" is for: it records zerocrossing events (when the observed value
38  * crosses the setpoint). From that data, we can determine the average oscillating period. This is the P_u of the
39  * ZN-method.
40  *
41  * Finally, we need to determine K_u, the ultimate gain. It turns out we can calculate this based on the amplitude of
42  * oscillation ("induced amplitude `a`) as described in [1]:
43  * K_u = (4d) / (πa)
44  * where d is the magnitude of the relay function (in range -d to +d).
45  * To measure `a`, we look at the current phase the relay function is in - if it's in the "heating" phase, then we
46  * expect the lowest temperature (=highest error) to be found in the phase because the peak will always happen slightly
47  * after the relay function has changed state (assuming a delay-dominated process).
48  *
49  * Finally, we use some heuristics to determine if the data we've received so far is good:
50  * - First, of course we must have enough data to calculate the values.
51  * - The ZC events need to happen at a relatively periodic rate. If the heating/cooling speeds are very different,
52  * I've observed the ZN parameters are not very useful.
53  * - The induced amplitude should not deviate too much. If the amplitudes deviate too much this means there has
54  * been some outside influence (or noise) on the system, and the measured amplitude values are not reliable.
55  *
56  * There are many ways this method can be improved, but on my simulation data the current method already produces very
57  * good results. Some ideas for future improvements:
58  * - Relay Function improvements:
59  * - Integrator, Preload, Saturation Relay ([1])
60  * - Use phase of measured signal relative to relay function.
61  * - Apply PID parameters from ZN, but continuously tweak them in a second step.
62  *
63  * [1]: https://warwick.ac.uk/fac/cross_fac/iatl/reinvention/archive/volume5issue2/hornsey/
64  * [2]: http://www.mstarlabs.com/control/znrule.html
65  * [3]: https://www.academia.edu/38620114/SEBORG_3rd_Edition_Process_Dynamics_and_Control
66  */
67 
68 PIDAutotuner::PIDAutotuneResult PIDAutotuner::update(float setpoint, float process_variable) {
70  if (this->state_ == AUTOTUNE_SUCCEEDED) {
72  return res;
73  }
74 
75  if (!std::isnan(this->setpoint_) && this->setpoint_ != setpoint) {
76  ESP_LOGW(TAG, "Setpoint changed during autotune! The result will not be accurate!");
77  }
78  this->setpoint_ = setpoint;
79 
80  float error = setpoint - process_variable;
81  const uint32_t now = millis();
82 
83  float output = this->relay_function_.update(error);
84  this->frequency_detector_.update(now, error);
85  this->amplitude_detector_.update(error, this->relay_function_.state);
86  res.output = output;
87 
89  // not enough data for calculation yet
90  ESP_LOGV(TAG, " Not enough data yet for aututuner");
91  return res;
92  }
93 
94  bool zc_symmetrical = this->frequency_detector_.is_increase_decrease_symmetrical();
95  bool amplitude_convergent = this->frequency_detector_.is_increase_decrease_symmetrical();
96  if (!zc_symmetrical || !amplitude_convergent) {
97  // The frequency/amplitude is not fully accurate yet, try to wait
98  // until the fault clears, or terminate after a while anyway
99  if (zc_symmetrical) {
100  ESP_LOGVV(TAG, " ZC is not symmetrical");
101  }
102  if (amplitude_convergent) {
103  ESP_LOGVV(TAG, " Amplitude is not convergent");
104  }
105  uint32_t phase = this->relay_function_.phase_count;
106  ESP_LOGVV(TAG, " Phase %u, enough=%u", phase, enough_data_phase_);
107 
108  if (this->enough_data_phase_ == 0) {
109  this->enough_data_phase_ = phase;
110  } else if (phase - this->enough_data_phase_ <= 6) {
111  // keep trying for at least 6 more phases
112  return res;
113  } else {
114  // proceed to calculating PID parameters
115  // warning will be shown in "Checks" section
116  }
117  }
118 
119  ESP_LOGI(TAG, "PID Autotune finished!");
120 
121  float osc_ampl = this->amplitude_detector_.get_mean_oscillation_amplitude();
122  float d = (this->relay_function_.output_positive - this->relay_function_.output_negative) / 2.0f;
123  ESP_LOGVV(TAG, " Relay magnitude: %f", d);
124  this->ku_ = 4.0f * d / float(M_PI * osc_ampl);
126 
127  this->state_ = AUTOTUNE_SUCCEEDED;
129  this->dump_config();
130 
131  return res;
132 }
134  ESP_LOGI(TAG, "PID Autotune:");
135  if (this->state_ == AUTOTUNE_SUCCEEDED) {
136  ESP_LOGI(TAG, " State: Succeeded!");
137  bool has_issue = false;
139  ESP_LOGW(TAG, " Could not reliable determine oscillation amplitude, PID parameters may be inaccurate!");
140  ESP_LOGW(TAG, " Please make sure you eliminate all outside influences on the measured temperature.");
141  has_issue = true;
142  }
144  ESP_LOGW(TAG, " Oscillation Frequency is not symmetrical. PID parameters may be inaccurate!");
145  ESP_LOGW(
146  TAG,
147  " This is usually because the heat and cool processes do not change the temperature at the same rate.");
148  ESP_LOGW(TAG,
149  " Please try reducing the positive_output value (or increase negative_output in case of a cooler)");
150  has_issue = true;
151  }
152  if (!has_issue) {
153  ESP_LOGI(TAG, " All checks passed!");
154  }
155 
156  auto fac = get_ziegler_nichols_pid_();
157  ESP_LOGI(TAG, " Calculated PID parameters (\"Ziegler-Nichols PID\" rule):");
158  ESP_LOGI(TAG, " ");
159  ESP_LOGI(TAG, " control_parameters:");
160  ESP_LOGI(TAG, " kp: %.5f", fac.kp);
161  ESP_LOGI(TAG, " ki: %.5f", fac.ki);
162  ESP_LOGI(TAG, " kd: %.5f", fac.kd);
163  ESP_LOGI(TAG, " ");
164  ESP_LOGI(TAG, " Please copy these values into your YAML configuration! They will reset on the next reboot.");
165 
166  ESP_LOGV(TAG, " Oscillation Period: %f", this->frequency_detector_.get_mean_oscillation_period());
167  ESP_LOGV(TAG, " Oscillation Amplitude: %f", this->amplitude_detector_.get_mean_oscillation_amplitude());
168  ESP_LOGV(TAG, " Ku: %f, Pu: %f", this->ku_, this->pu_);
169 
170  ESP_LOGD(TAG, " Alternative Rules:");
171  // http://www.mstarlabs.com/control/znrule.html
172  print_rule_("Ziegler-Nichols PI", 0.45f, 0.54f, 0.0f);
173  print_rule_("Pessen Integral PID", 0.7f, 1.75f, 0.105f);
174  print_rule_("Some Overshoot PID", 0.333f, 0.667f, 0.111f);
175  print_rule_("No Overshoot PID", 0.2f, 0.4f, 0.0625f);
176  }
177 
178  if (this->state_ == AUTOTUNE_RUNNING) {
179  ESP_LOGI(TAG, " Autotune is still running!");
180  ESP_LOGD(TAG, " Status: Trying to reach %.2f °C", setpoint_ - relay_function_.current_target_error());
181  ESP_LOGD(TAG, " Stats so far:");
182  ESP_LOGD(TAG, " Phases: %u", relay_function_.phase_count);
183  ESP_LOGD(TAG, " Detected %u zero-crossings", frequency_detector_.zerocrossing_intervals.size()); // NOLINT
184  ESP_LOGD(TAG, " Current Phase Min: %.2f, Max: %.2f", amplitude_detector_.phase_min,
186  }
187 }
188 PIDAutotuner::PIDResult PIDAutotuner::calculate_pid_(float kp_factor, float ki_factor, float kd_factor) {
189  float kp = kp_factor * ku_;
190  float ki = ki_factor * ku_ / pu_;
191  float kd = kd_factor * ku_ * pu_;
192  return {
193  .kp = kp,
194  .ki = ki,
195  .kd = kd,
196  };
197 }
198 void PIDAutotuner::print_rule_(const char *name, float kp_factor, float ki_factor, float kd_factor) {
199  auto fac = calculate_pid_(kp_factor, ki_factor, kd_factor);
200  ESP_LOGD(TAG, " Rule '%s':", name);
201  ESP_LOGD(TAG, " kp: %.5f, ki: %.5f, kd: %.5f", fac.kp, fac.ki, fac.kd);
202 }
203 
204 // ================== RelayFunction ==================
206  if (this->state == RELAY_FUNCTION_INIT) {
207  bool pos = error > this->noiseband;
208  state = pos ? RELAY_FUNCTION_POSITIVE : RELAY_FUNCTION_NEGATIVE;
209  }
210  bool change = false;
211  if (this->state == RELAY_FUNCTION_POSITIVE && error < -this->noiseband) {
212  // Positive hysteresis reached, change direction
213  this->state = RELAY_FUNCTION_NEGATIVE;
214  change = true;
215  } else if (this->state == RELAY_FUNCTION_NEGATIVE && error > this->noiseband) {
216  // Negative hysteresis reached, change direction
217  this->state = RELAY_FUNCTION_POSITIVE;
218  change = true;
219  }
220 
221  float output = state == RELAY_FUNCTION_POSITIVE ? output_positive : output_negative;
222  if (change) {
223  this->phase_count++;
224  ESP_LOGV(TAG, "Autotune: Turning output to %.1f%%", output * 100);
225  }
226 
227  return output;
228 }
229 
230 // ================== OscillationFrequencyDetector ==================
231 void PIDAutotuner::OscillationFrequencyDetector::update(uint32_t now, float error) {
232  if (this->state == FREQUENCY_DETECTOR_INIT) {
233  bool pos = error > this->noiseband;
234  state = pos ? FREQUENCY_DETECTOR_POSITIVE : FREQUENCY_DETECTOR_NEGATIVE;
235  }
236 
237  bool had_crossing = false;
238  if (this->state == FREQUENCY_DETECTOR_POSITIVE && error < -this->noiseband) {
239  this->state = FREQUENCY_DETECTOR_NEGATIVE;
240  had_crossing = true;
241  } else if (this->state == FREQUENCY_DETECTOR_NEGATIVE && error > this->noiseband) {
242  this->state = FREQUENCY_DETECTOR_POSITIVE;
243  had_crossing = true;
244  }
245 
246  if (had_crossing) {
247  // Had crossing above hysteresis threshold, record
248  ESP_LOGV(TAG, "Autotune: Detected Zero-Cross at %u", now);
249  if (this->last_zerocross != 0) {
250  uint32_t dt = now - this->last_zerocross;
251  ESP_LOGV(TAG, " dt: %u", dt);
252  this->zerocrossing_intervals.push_back(dt);
253  }
254  this->last_zerocross = now;
255  }
256 }
258  // Do we have enough data in this detector to generate PID values?
259  return this->zerocrossing_intervals.size() >= 2;
260 }
262  // Get the mean oscillation period in seconds
263  // Only call if has_enough_data() has returned true.
264  float sum = 0.0f;
265  for (uint32_t v : this->zerocrossing_intervals)
266  sum += v;
267  // zerocrossings are each half-period, multiply by 2
268  float mean_value = sum / this->zerocrossing_intervals.size();
269  // divide by 1000 to get seconds, multiply by two because zc happens two times per period
270  float mean_period = mean_value / 1000 * 2;
271  return mean_period;
272 }
274  // Check if increase/decrease of process value was symmetrical
275  // If the process value increases much faster than it decreases, the generated PID values will
276  // not be very good and the function output values need to be adjusted
277  // Happens for example with a well-insulated heating element.
278  // We calculate this based on the zerocrossing interval.
279  if (zerocrossing_intervals.empty())
280  return false;
281  uint32_t max_interval = zerocrossing_intervals[0];
282  uint32_t min_interval = zerocrossing_intervals[0];
283  for (uint32_t interval : zerocrossing_intervals) {
284  max_interval = std::max(max_interval, interval);
285  min_interval = std::min(min_interval, interval);
286  }
287  float ratio = min_interval / float(max_interval);
288  return ratio >= 0.66;
289 }
290 
291 // ================== OscillationAmplitudeDetector ==================
294  if (relay_state != last_relay_state) {
295  if (last_relay_state == RelayFunction::RELAY_FUNCTION_POSITIVE) {
296  // Transitioned from positive error to negative error.
297  // The positive error peak must have been in previous segment (180° shifted)
298  // record phase_max
299  this->phase_maxs.push_back(phase_max);
300  ESP_LOGV(TAG, "Autotune: Phase Max: %f", phase_max);
301  } else if (last_relay_state == RelayFunction::RELAY_FUNCTION_NEGATIVE) {
302  // Transitioned from negative error to positive error.
303  // The negative error peak must have been in previous segment (180° shifted)
304  // record phase_min
305  this->phase_mins.push_back(phase_min);
306  ESP_LOGV(TAG, "Autotune: Phase Min: %f", phase_min);
307  }
308  // reset phase values for next phase
309  this->phase_min = error;
310  this->phase_max = error;
311  }
312  this->last_relay_state = relay_state;
313 
314  this->phase_min = std::min(this->phase_min, error);
315  this->phase_max = std::max(this->phase_max, error);
316 
317  // Check arrays sizes, we keep at most 7 items (6 datapoints is enough, and data at beginning might not
318  // have been stabilized)
319  if (this->phase_maxs.size() > 7)
320  this->phase_maxs.erase(this->phase_maxs.begin());
321  if (this->phase_mins.size() > 7)
322  this->phase_mins.erase(this->phase_mins.begin());
323 }
325  // Return if we have enough data to generate PID parameters
326  // The first phase is not very useful if the setpoint is not set to the starting process value
327  // So discard first phase. Otherwise we need at least two phases.
328  return std::min(phase_mins.size(), phase_maxs.size()) >= 3;
329 }
331  float total_amplitudes = 0;
332  size_t total_amplitudes_n = 0;
333  for (size_t i = 1; i < std::min(phase_mins.size(), phase_maxs.size()) - 1; i++) {
334  total_amplitudes += std::abs(phase_maxs[i] - phase_mins[i + 1]);
335  total_amplitudes_n++;
336  }
337  float mean_amplitude = total_amplitudes / total_amplitudes_n;
338  // Amplitude is measured from center, divide by 2
339  return mean_amplitude / 2.0f;
340 }
342  // Check if oscillation amplitude is convergent
343  // We implement this by checking global extrema against average amplitude
344  if (this->phase_mins.empty() || this->phase_maxs.empty())
345  return false;
346 
347  float global_max = phase_maxs[0], global_min = phase_mins[0];
348  for (auto v : this->phase_mins)
349  global_min = std::min(global_min, v);
350  for (auto v : this->phase_maxs)
351  global_max = std::min(global_max, v);
352  float global_amplitude = (global_max - global_min) / 2.0f;
353  float mean_amplitude = this->get_mean_oscillation_amplitude();
354  return (mean_amplitude - global_amplitude) / (global_amplitude) < 0.05f;
355 }
356 
357 } // namespace pid
358 } // namespace esphome
struct esphome::pid::PIDAutotuner::OscillationAmplitudeDetector amplitude_detector_
const char * name
Definition: stm32flash.h:78
struct esphome::pid::PIDAutotuner::OscillationFrequencyDetector frequency_detector_
enum esphome::pid::PIDAutotuner::RelayFunction::RelayFunctionState state
struct esphome::pid::PIDAutotuner::RelayFunction relay_function_
PIDResult calculate_pid_(float kp_factor, float ki_factor, float kd_factor)
uint32_t IRAM_ATTR HOT millis()
Definition: core.cpp:26
PIDAutotuneResult update(float setpoint, float process_variable)
void update(float error, RelayFunction::RelayFunctionState relay_state)
void print_rule_(const char *name, float kp_factor, float ki_factor, float kd_factor)
Definition: a4988.cpp:4
enum esphome::pid::PIDAutotuner::State state_
bool state
Definition: fan.h:34
PIDResult get_ziegler_nichols_pid_()
Definition: pid_autotuner.h:97