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