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>&);
|
