From: Lukáš Jiřiště Date: Mon, 11 Nov 2024 16:00:11 +0000 (+0100) Subject: Create Measurement to track pressure over time X-Git-Url: https://git.ljiriste.work/?a=commitdiff_plain;h=b9517f93f1e69d34b2aebb1907ed3e696f4a53f1;p=Servomatic.git Create Measurement to track pressure over time This is an untested implementation. The new class should provide the functionality necessary to conveniently control the cell. Minor changes for compatibility and eventual testing in other files. --- diff --git a/Makefile b/Makefile index 3bf682c..f7bc4e7 100644 --- a/Makefile +++ b/Makefile @@ -15,6 +15,7 @@ endif SRCS := Servomatic.ino \ PressureSensor.ino \ Valve.ino \ + Measurement.ino \ NAME := Servomatic diff --git a/Measurement.h b/Measurement.h new file mode 100644 index 0000000..d5b27d4 --- /dev/null +++ b/Measurement.h @@ -0,0 +1,39 @@ +#ifndef MEASUREMENT_H +# define MEASUREMENT_H + +#include "PressureSensor.h" + +class Measurement +{ + private: + static const size_t MAX_ELEMENTS = 100; + static const unsigned long START_INTERVAL_MS = 5000; + + const PressureSensor &m_sensor; + unsigned long m_interval_ms; + unsigned long m_start_time; + float m_pressure[MAX_ELEMENTS]; + unsigned long m_time[MAX_ELEMENTS]; + size_t m_size; + + float calc_equilibrium_dist() const; + float calc_decay_speed(float equilibrium_pressure) const; + float get_error(float equilibrium_pressure) const; + float get_error_diff(float equilibrium_pressure, float step = 0.001) const; + float get_equilibrium_pressure(size_t max_iter = 10) const; + + public: + Measurement() = delete; + Measurement(const Measurement &other) = delete; + Measurement(const PressureSensor &sensor); + ~Measurement(); + + Measurement &operator=(const Measurement &other) = delete; + + void get_measurement(); + void thin_out(); + void clear(); + void measure_until(float threshold = 0.01); +}; + +#endif // MEASUREMENT_H diff --git a/Measurement.ino b/Measurement.ino new file mode 100644 index 0000000..3fe75c2 --- /dev/null +++ b/Measurement.ino @@ -0,0 +1,155 @@ +#include "Measurement.h" + +Measurement::Measurement(const PressureSensor &sensor) + : m_size{0} + , m_sensor{sensor} + , m_interval_ms{START_INTERVAL_MS} +{ +} + +Measurement::~Measurement() +{ +} + +void Measurement::get_measurement() +{ + if (m_size == MAX_ELEMENTS) + thin_out(); + m_pressure[m_size] = m_sensor.get_pressure(); + m_time[m_size] = millis() - m_start_time; +} + +void Measurement::thin_out() +{ + m_interval_ms *= 2; + m_size /= 2; + for (size_t i{1}; i < m_size; ++i) + { + m_pressure[i] = m_pressure[2 * i]; + m_time[i] = m_time[2 * i]; + } +} + +void Measurement::clear() +{ + m_size = 0; +} + +// The function calls get_measurement repeatedly. It uses these measurements to +// estimate how much of the process (sorption) is done. +// +// The loop last as long as the following condition holds: +// +// p_c - p_e +// --------- >= threshold +// p_0 - p_e +// +// where p_c is the current pressure, p_e is the expected equilibrium pressure +// and p_0 is the start pressure. The fraction could be interpreted as the +// normalized distance from equilibrium. +// +// I decided to use an exponential decay model for the pressure dependence on time +// +// p(t) = (p_0 - p_e) * exp(-a * t) + p_e +// +// as I find it quite reasonable. The "a" constant encodes the decay speed. +// Using this model the above condition reduces to +// +// exp(-a * t) >= threshold +// +// The model has 2 parameters, and it is difficult to get them analytically +// using the least squares method. But getting the optimal "a" parameter for a +// given p_e is quite easy using least squares, and it gives this: +// +// a(p_e) = sum(ln((p_i - p_e) / (p_0 - p_e)) * t_i) / sum((t_i)^2) +// +// Using this we can then numerically optimize p_e. Interval halving is used +// because it is foolproof, the precision needed is low and the evaluation is +// not computationally demanding. +void Measurement::measure_until(float threshold = 0.01) +{ + while (m_size < 20 || threshold < calc_equilibrium_dist()) + { + get_measurement(); + delay(m_interval_ms); + } +} + +float Measurement::calc_equilibrium_dist() const +{ + float decay_speed{calc_decay_speed(get_equilibrium_pressure())}; + + return (exp(-decay_speed * m_time[m_size - 1])); +} + +float Measurement::calc_decay_speed(float equilibrium_pressure) const +{ + float time_squares_sum{0}; + float time_logp_sum{0}; + + for (size_t i{0}; i < m_size; ++i) + { + time_squares_sum += static_cast(m_time[i]) * m_time[i]; + time_logp_sum += m_time[i] * log((m_pressure[i] - equilibrium_pressure) + / (m_pressure[0] - equilibrium_pressure)); + } + return (- time_logp_sum / time_squares_sum); +} + +float Measurement::get_error(float eq_pressure) const +{ + float decay_speed{calc_decay_speed(eq_pressure)}; + float residual{}; + float error{0}; + + for (size_t i{0}; i < m_size; ++i) + { + residual = m_pressure[i] - (m_pressure[0] - eq_pressure) * exp(- decay_speed * m_time[i]) - eq_pressure; + error += residual * residual; + } + return (error); +} + +float Measurement::get_error_diff(float eq_pressure, float step) const +{ + float error_minus{}; + float error_plus{}; + + error_plus = get_error(eq_pressure + step); + if (eq_pressure < step) + { + error_minus = get_error(0); + return ((error_plus - error_minus) / (step + eq_pressure)); + } + error_minus = get_error(eq_pressure - step); + return ((error_plus - error_minus) / (2 * step)); +} + +float Measurement::get_equilibrium_pressure(size_t max_iter) const +{ + float min_arg{0}; + float min_val{get_error_diff(min_arg)}; + float max_arg{1.1 * m_sensor.m_max_pressure}; + float max_val{get_error_diff(max_arg)}; + float test{}; + float test_val{}; + + for (; max_iter > 0; --max_iter) + { + test = (max_arg + min_arg) / 2; + test_val = get_error_diff(test); + if (test_val == 0) + return (test); + else if ((test_val < 0 && min_val < 0) || (test_val > 0 && min_val > 0)) + { + min_arg = test; + min_val = test_val; + } + else + { + max_arg = test; + max_val = test_val; + } + } + return ((max_arg + min_arg) / 2); +} diff --git a/PressureSensor.h b/PressureSensor.h index a5287f6..f8a857b 100644 --- a/PressureSensor.h +++ b/PressureSensor.h @@ -7,7 +7,6 @@ class PressureSensor typedef int pin; const pin m_analog_pin; - const float m_max_pressure; const float m_resistance; // in Ohms static const float MIN_CURRENT = 4e-3; // in A @@ -18,6 +17,8 @@ class PressureSensor float signal_to_pressure(int ad_signal); public: + const float m_max_pressure; + PressureSensor() = delete; PressureSensor(pin analog_pin, float max_pressure = 10, float resistance = 240); PressureSensor(const PressureSensor &other) = delete; diff --git a/Servomatic.ino b/Servomatic.ino index c33f361..2f721e8 100644 --- a/Servomatic.ino +++ b/Servomatic.ino @@ -22,15 +22,16 @@ int main() void main_loop() { - PressureSensor test{A5, 10}; + PressureSensor sensor_1{PRESSURE_SENSOR1_PIN, 10}; Valve test_valve{6, 760, 2400}; int usec; + Measurement measurement{sensor_1}; while (1) { - delay(5000); + measurement.measure_until(); test_valve.open(10); - delay(5000); + measurement.measure_until(); test_valve.close(100); } }