Line data Source code
1 0 : // Distributed under the MIT License. 2 : // See LICENSE.txt for details. 3 : 4 : #pragma once 5 : 6 : #include <array> 7 : #include <cstddef> 8 : #include <deque> 9 : #include <optional> 10 : #include <pup.h> 11 : #include <pup_stl.h> 12 : 13 : #include "DataStructures/DataVector.hpp" 14 : #include "Options/String.hpp" 15 : #include "Utilities/Gsl.hpp" 16 : #include "Utilities/Serialization/PupStlCpp17.hpp" 17 : 18 : /// \ingroup ControlSystemGroup 19 : /// A weighted exponential averager of \f$Q\f$ and its derivatives 20 : /// implementing Appendix A in \cite Hemberger2012jz 21 : /// 22 : /// The purpose of the averager is to provide \f$Q\f$ and 'smoothed' numerical 23 : /// derivatives of \f$Q\f$ up to the `DerivOrder`'th derivative. The 0th 24 : /// derivative of \f$Q\f$ is not typically averaged, since no differencing is 25 : /// necessary and no noise is introduced. The `average_0th_deriv_of_q` option 26 : /// allows for specifying that the 0th derivative should be averaged in addition 27 : /// to the derivatives. This may be desirable for control systems where 28 : /// \f$Q\f$ contains some 'noise', e.g. size control typically uses an averaged 29 : /// \f$Q\f$ since the raw computed \f$Q\f$ is a function of the minimum on a 30 : /// surface and may jump around discontinuously. 31 : /// 32 : /// The averager is designed to support DerivOrders 1, 2, and 3. If an 33 : /// additional DerivOrder is needed, finite differencing needs to be implemented 34 : /// to specifically handle that order (as it seems like overkill to generalize 35 : /// the differencing stencil at this time). 36 : template <size_t DerivOrder> 37 1 : class Averager { 38 : static_assert(DerivOrder != 0, "Cannot average just a function value."); 39 : static_assert( 40 : DerivOrder <= 3, 41 : "Averager only instantiated for deriv order less than or equal to 3. If " 42 : "you want a higher deriv order, you'll need to add it yourself."); 43 : 44 : public: 45 0 : struct AverageTimescaleFraction { 46 0 : using type = double; 47 0 : static constexpr Options::String help = { 48 : "Time scale of exponential averaging"}; 49 : }; 50 : 51 0 : struct Average0thDeriv { 52 0 : using type = bool; 53 0 : static constexpr Options::String help = { 54 : "Whether to average the 0th derivative"}; 55 : }; 56 : 57 0 : using options = tmpl::list<AverageTimescaleFraction, Average0thDeriv>; 58 0 : static constexpr Options::String help{ 59 : "Averager: Performs exponential averaging of the control signal at " 60 : "multiple times in order to provide smoother derivatives of the control " 61 : "signal."}; 62 : 63 : /// `avg_timescale_frac` determines the exponential averaging timescale 64 : /// through \f$\tau_\mathrm{avg} = \f$`avg_timescale_frac`\f$\times \tau\f$, 65 : /// where \f$\tau\f$ is the damping time. 66 : /// `avg_timescale_frac` must be positive. 67 : /// `average_0th_deriv_of_q` determines whether the call operator returns an 68 : /// averaged or unaveraged quantity for the 0th derivative of \f$Q\f$. 69 : /// `true` returns an averaged 0th derivative of \f$Q\f$. 70 : /// `false` returns the raw 0th derivative of \f$Q\f$. 71 : /// The derivatives are always averaged (to reduce noise due to numerical 72 : /// differentiation), so the `average_0th_deriv_of_q` option only specifies 73 : /// whether to return an averaged value for the 0th derivative piece of the 74 : /// function. 75 1 : Averager(double avg_timescale_frac, bool average_0th_deriv_of_q); 76 : 77 0 : Averager() = default; 78 : // Explicitly defined move constructor due to the fact that the std::deque 79 : // move constructor is not marked 80 0 : Averager(Averager&& rhs); 81 0 : Averager& operator=(Averager&& rhs); 82 0 : Averager(const Averager&) = default; 83 0 : Averager& operator=(const Averager&) = default; 84 0 : ~Averager() = default; 85 : 86 : template <size_t DDerivOrder> // NOLINT 87 0 : friend bool operator==(const Averager<DDerivOrder>&, // NOLINT 88 : const Averager<DDerivOrder>&); // NOLINT 89 : 90 : /// Returns \f$Q\f$ and its derivatives at \f$t=\f$`time`, provided there is 91 : /// sufficient data. The averager is limited by the need for at least 92 : /// (`DerivOrder` + 1) data points in order to provide the `DerivOrder`'th 93 : /// derivative. If sufficient data is available, it returns \f$Q\f$ and 94 : /// averaged derivatives of \f$Q\f$ up to the `DerivOrder`'th derivative, at 95 : /// \f$t=\f$`time`. If `using_average_0th_deriv_of_q()` is `true`, then the 96 : /// returned 0th derivative of \f$Q\f$ is also averaged. In the case that 97 : /// there is insufficient data, the operator returns an invalid 98 : /// `std::optional`. 99 1 : const std::optional<std::array<DataVector, DerivOrder + 1>>& operator()( 100 : double time) const; 101 : /// A function that allows for resetting the averager. 102 1 : void clear(); 103 : /// The function responsible for updating the averaged values 104 : /// at \f$t=\f$`time`. Requires `raw_q` (the raw components of \f$Q(t)\f$) 105 : /// and `timescales` (the associated damping times for each component). 106 1 : void update(double time, const DataVector& raw_q, 107 : const DataVector& timescales); 108 : /// Returns the latest time at which the averager has sufficient data to 109 : /// return \f$Q\f$ and its derivatives. 110 1 : double last_time_updated() const; 111 : /// Returns the exponentially averaged time at \f$t=\f$`time`. The time is 112 : /// averaged along side \f$Q\f$ to determine the effective time at which 113 : /// the average value is computed. The effective time is retarded, due to the 114 : /// weighting of past times. 115 1 : double average_time(double time) const; 116 : /// Returns a bool corresponding to whether `average_0th_deriv_of_q` 117 : /// is `true`/`false`. 118 1 : bool using_average_0th_deriv_of_q() const { return average_0th_deriv_of_q_; }; 119 : /// Returns the averaging timescale fraction 120 1 : double avg_timescale_frac() const { return avg_tscale_frac_; } 121 : 122 : /// Assign the minimum of the measurement timescales related to a specific 123 : /// control system to time_between_measurements_ 124 1 : void assign_time_between_measurements( 125 : const double current_time_between_measurements) { 126 : time_between_measurements_ = current_time_between_measurements; 127 : } 128 : 129 : /// Returns `true` if the averager is ready to receive a measurement 130 1 : bool is_ready(const double time) const { 131 : if (UNLIKELY(times_.empty())) { 132 : return true; 133 : } 134 : return time >= times_[0] + time_between_measurements_; 135 : } 136 : 137 0 : void pup(PUP::er& p); 138 : 139 : private: 140 : /// Returns the function and numerical derivatives up to the 141 : /// `DerivOrder`'th derivative using backward finite differencing that 142 : /// supports non-uniform time spacing. Since we use a stencil size of 143 : /// (`DerivOrder` + 1), the order of the finite difference varies for 144 : /// different derivatives. Assuming \f$M\lte\f$`DerivOrder`, the \f$M\f$'th 145 : /// derivative is of order (`DerivOrder` + 1 - \f$M\f$). 146 1 : std::array<DataVector, DerivOrder + 1> get_derivs() const; 147 : 148 0 : double avg_tscale_frac_; 149 0 : bool average_0th_deriv_of_q_; 150 0 : std::optional<std::array<DataVector, DerivOrder + 1>> averaged_values_{}; 151 0 : std::optional<std::array<DataVector, DerivOrder + 1>> nullopt_ = std::nullopt; 152 0 : std::deque<double> times_; 153 0 : std::deque<DataVector> raw_qs_; 154 0 : double weight_k_ = 0.0; 155 0 : double tau_k_ = 0.0; 156 : // If this time_between_measurements_ isn't set, the default should just be 157 : // that no measurements happen (i.e. infinity) 158 0 : double time_between_measurements_{std::numeric_limits<double>::infinity()}; 159 : }; 160 : 161 : template <size_t DDerivOrder> 162 0 : bool operator!=(const Averager<DDerivOrder>&, const Averager<DDerivOrder>&);