From: Lukáš Jiřiště Date: Wed, 13 Nov 2024 11:01:51 +0000 (+0100) Subject: Reverse the role of p_e and a in optimization X-Git-Url: https://git.ljiriste.work/?a=commitdiff_plain;h=b9b00f24c837ec90efd4ce78fefe0398c6c74cf9;p=Servomatic.git Reverse the role of p_e and a in optimization The p_e can be calculated for a fixed a without logarithmization, so the minimum should be the true minimum. Also some other minor changes to the function around were made. --- diff --git a/Measurement.h b/Measurement.h index 85a32ea..4b5c596 100644 --- a/Measurement.h +++ b/Measurement.h @@ -17,10 +17,11 @@ class Measurement 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; + float calc_equilibrium_pressure(float equilibrium_pressure) const; + float get_error(float decay_speed) const; + float get_error_diff(float decay_speed, float step) const; + float get_decay_speed(size_t max_iter = 20) const; + float screen_for_minimum(float upper_bound) const; public: Measurement() = delete; diff --git a/Measurement.ino b/Measurement.ino index 082a45a..6eadc2d 100644 --- a/Measurement.ino +++ b/Measurement.ino @@ -61,16 +61,18 @@ void Measurement::clear() // 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: +// using the least squares method. But getting the optimal "p_e" parameter for a +// given decay_speed "a" 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) +// p_e(a) = sum((p_i - p_0 * exp_i) * (1 - exp_i)) / sum((1 - exp_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. +// where exp_i = exp(-a * t_i). +// +// Using this we can then numerically optimize a. Interval halving is used +// because it is foolproof and the evaluation is not computationally demanding. void Measurement::wait_for_equilibrium(float threshold = 0.01) { + clear(); while (m_size < 20 || threshold < calc_equilibrium_dist()) { get_measurement(); @@ -80,67 +82,107 @@ void Measurement::wait_for_equilibrium(float threshold = 0.01) float Measurement::calc_equilibrium_dist() const { - float decay_speed{calc_decay_speed(get_equilibrium_pressure())}; - - return (exp(-decay_speed * m_time[m_size - 1])); + return (exp(-get_decay_speed() * m_time[m_size - 1])); } -float Measurement::calc_decay_speed(float equilibrium_pressure) const +float Measurement::calc_equilibrium_pressure(float decay_speed) const { - float time_squares_sum{0}; - float time_logp_sum{0}; + float numerator_sum{0}; + float denominator_sum{0}; + float exponential_term{}; 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)); + exponential_term = exp(-decay_speed * m_time[i]); + numerator_sum += (m_pressure[i] - m_pressure[0] * exponential_term) + * (1 - exponential_term); + denominator_sum += (1 - exponential_term) * (1 - exponential_term); } - return (- time_logp_sum / time_squares_sum); + return (numerator_sum / denominator_sum); } -float Measurement::get_error(float eq_pressure) const +float Measurement::get_error(float decay_speed) const { - float decay_speed{calc_decay_speed(eq_pressure)}; + float eq_pressure{calc_equilibrium_pressure(decay_speed)}; float residual{}; float error{0}; - for (size_t i{0}; i < m_size; ++i) + if (decay_speed > 0) { - residual = m_pressure[i] - (m_pressure[0] - eq_pressure) * exp(- decay_speed * m_time[i]) - eq_pressure; - error += residual * residual; + 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; + } + } + else + { + for (size_t i{0}; i < m_size; ++i) + { + residual = m_pressure[i] - eq_pressure; + error += residual * residual; + } } return (error); } -float Measurement::get_error_diff(float eq_pressure, float step) const +float Measurement::get_error_diff(float decay_speed, float step) const { float error_minus{}; float error_plus{}; - error_plus = get_error(eq_pressure + step); - if (eq_pressure < step) + error_plus = get_error(decay_speed + step); + if (decay_speed < step) { error_minus = get_error(0); - return ((error_plus - error_minus) / (step + eq_pressure)); + return ((error_plus - error_minus) / (step + decay_speed)); } - error_minus = get_error(eq_pressure - step); + error_minus = get_error(decay_speed - step); return ((error_plus - error_minus) / (2 * step)); } -float Measurement::get_equilibrium_pressure(size_t max_iter) const +float Measurement::screen_for_minimum(float upper_bound) 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_arg{upper_bound}; + float last_arg{}; + float test_value{}; + float last_value{}; + + test_value = get_error(test_arg); + do + { + last_value = test_value; + last_arg = test_arg; + test_arg /= 10; + test_value = get_error(test_arg); + } + while (test_value <= last_value); + return (last_arg); +} + +float Measurement::get_decay_speed(size_t max_iter) const +{ + float min_arg{}; + float min_val{}; + float max_arg{}; + float max_val{}; float test{}; float test_val{}; + min_arg = screen_for_minimum(1); + max_arg = min_arg; + do + { + min_arg /= 2; + max_arg *= 2; + min_val = get_error_diff(min_arg, min_arg / 1000); + max_val = get_error_diff(max_arg, max_arg / 1000); + } + while ((min_val > 0 && max_val > 0) || (min_val < 0 && max_val < 0)); for (; max_iter > 0; --max_iter) { test = (max_arg + min_arg) / 2; - test_val = get_error_diff(test); + test_val = get_error_diff(test, test / 1000); if (test_val == 0) return (test); else if ((test_val < 0 && min_val < 0) || (test_val > 0 && min_val > 0))