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