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 3bf682cb3bbb16bc5cc572669bbe9d55607c0f96..f7bc4e7089ba95af3fd784411fcc495439fba8b0 100644 (file)
--- 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 (file)
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 (file)
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<float>(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 a5287f6f2e816f7233943f23db470bf362ffe0e5..f8a857b0600c340a02e99423377470d4b98f0881 100644 (file)
--- 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 c33f3614b16426d71280e70ec5b1853931c30e6c..2f721e84dea9457f865627857df8bd03c6b86f20 100644 (file)
--- 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);
        }
 }