Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <algorithm>
7 : #include <memory>
8 : #include <string>
9 : #include <type_traits>
10 : #include <unordered_map>
11 : #include <unordered_set>
12 :
13 : #include "ControlSystem/CalculateMeasurementTimescales.hpp"
14 : #include "ControlSystem/CombinedName.hpp"
15 : #include "ControlSystem/Metafunctions.hpp"
16 : #include "ControlSystem/Tags/SystemTags.hpp"
17 : #include "DataStructures/DataVector.hpp"
18 : #include "Domain/Creators/DomainCreator.hpp"
19 : #include "Utilities/TMPL.hpp"
20 :
21 : namespace control_system {
22 : /*!
23 : * \ingroup ControlSystemGroup
24 : * \brief Calculate the next expiration time for the FunctionsOfTime.
25 : *
26 : * If \p delay_update is false, this returns
27 : *
28 : * \f{equation}
29 : * T_\mathrm{expr}^\mathrm{FoT} = t + N * \tau_\mathrm{m}^\mathrm{new},
30 : * \f}
31 : *
32 : * and if it is true,
33 : *
34 : * \f{equation}
35 : * T_\mathrm{expr}^\mathrm{FoT} = t + \tau_\mathrm{m}^\mathrm{old}
36 : * + N * \tau_\mathrm{m}^\mathrm{new},
37 : * \f}
38 : *
39 : * where \f$T_\mathrm{expr}^\mathrm{FoT}\f$ is the expiration time for the
40 : * FunctionsOfTime, \f$t\f$ is the update time,
41 : * \f$\tau_\mathrm{m}^\mathrm{old/new}\f$ is the measurement timescale, and
42 : * \f$N\f$ is the number of measurements per update.
43 : *
44 : * If \p delay_update is true, we update the functions of time
45 : * one (old) measurement before they actually expire.
46 : *
47 : * The choice of having the functions of time expire exactly one old measurement
48 : * after they are updated is arbitrary. They could expire any time between the
49 : * update time and one old measurement after the update. This decision was made
50 : * to minimize time spent waiting for the functions of time to be valid.
51 : *
52 : * Since functions of time are valid at their expiration time, we are actually
53 : * able to do the next measurement if the expiration time is at that
54 : * measurement. Thus we delay any potential waiting that may happen until the
55 : * subsequent measurement (and by that time, most, if not all, functions of time
56 : * should have been updated because an entire horizon find happened in the
57 : * meantime). If the expiration time was earlier than the next measurement, we'd
58 : * have to pause the Algorithm on the DG elements and wait until all the
59 : * functions of time have been updated.
60 : *
61 : * If \p delay_update is false, we update the functions of time
62 : * immediately, which allows for less parallelization but potentially
63 : * makes the control system more responsive.
64 : */
65 1 : double function_of_time_expiration_time(
66 : double time, const DataVector& old_measurement_timescales,
67 : const DataVector& new_measurement_timescales, int measurements_per_update,
68 : bool delay_update);
69 :
70 : /*!
71 : * \ingroup ControlSystemGroup
72 : * \brief Calculate the first expiration time for the FunctionsOfTime.
73 : *
74 : * The first expiration time is
75 : *
76 : * \f{equation}
77 : * T_\mathrm{expr}^\mathrm{FoT} = t + N * \tau_\mathrm{m},
78 : * \f}
79 : *
80 : * where $t$ is the initial time, $\tau_\mathrm{m}$ is the measurement
81 : * timescale, and \f$N\f$ is the number of measurements per update if
82 : * \p delay_update is true, and one less if it is false.
83 : */
84 1 : double function_of_time_initial_expiration_time(
85 : double time, const DataVector& measurement_timescales,
86 : int measurements_per_update, bool delay_update);
87 :
88 : /*!
89 : * \ingroup ControlSystemGroup
90 : * \brief Calculate the next expiration time for the MeasurementTimescales.
91 : *
92 : * \f{align}
93 : * T_\mathrm{expr}^\mathrm{m} &= T_\mathrm{expr}^\mathrm{FoT} - \frac{1}{2}
94 : * \tau_\mathrm{m}^\mathrm{new} \\
95 : * \f}
96 : *
97 : * where \f$T_\mathrm{expr}^\mathrm{m}\f$ is the measurement expiration time,
98 : * \f$T_\mathrm{expr}^\mathrm{FoT}\f$ is the function of time expiration time as
99 : * calculate from `function_of_time_expiration_time()`, and
100 : * \f$\tau_\mathrm{m}^\mathrm{new}\f$ is the new measurement timescale. The
101 : * reason for the factor of a half is as follows:
102 : *
103 : * We update the functions of time one (old) measurement before the expiration
104 : * time. Based on how dense triggers are set up, which control_system::Trigger
105 : * is a dense trigger, you calculate the next trigger (measurement) time at the
106 : * current measurement time. However, at the function of time expiration time we
107 : * need updated damping timescales from all control systems in order to
108 : * calculate when the next measurement is going to be (and in turn, the next
109 : * measurement expiration time). Thus, at the measurement that occurs at the
110 : * function of time expiration time, our measurement timescales can't be valid
111 : * and we must wait for updated ones. To achieve this, we set the measurement
112 : * expiration time *before* the function of time expiration time, but *after*
113 : * the previous measurement (the update measurement). The factor of one half is
114 : * just to guarantee we are more than epsilon before the function of time
115 : * expiration time and more than epsilon after the update measurement.
116 : */
117 1 : double measurement_expiration_time(double time,
118 : const DataVector& old_measurement_timescales,
119 : const DataVector& new_measurement_timescales,
120 : int measurements_per_update,
121 : bool delay_update);
122 :
123 : /*!
124 : * \ingroup ControlSystemGroup
125 : * \brief Calculate the first expiration time for the MeasurementTimescales.
126 : *
127 : * This is related to `function_of_time_initial_expiration_time()` in
128 : * the same way that `measurement_expiration_time()` is related to
129 : * `function_of_time_expiration_time()`.
130 : */
131 1 : double measurement_initial_expiration_time(
132 : double time, const DataVector& measurement_timescales,
133 : int measurements_per_update, bool delay_update);
134 :
135 : /*!
136 : * \ingroup ControlSystemGroup
137 : * \brief Construct the initial expiration times for functions of time that are
138 : * controlled by a control system
139 : *
140 : * The expiration times are constructed using inputs from control system
141 : * OptionHolders as an unordered map from the name of the function of time being
142 : * controlled to the expiration time.
143 : *
144 : * The expiration time for each individual function of time is computed as
145 : * $\tau_\mathrm{exp} = \alpha_\mathrm{update} \tau_\mathrm{damp}$ where
146 : * $\alpha_\mathrm{update}$ is the update fraction supplied as input to the
147 : * Controller and $\tau_\mathrm{damp}$ is/are the damping timescales
148 : * supplied from the TimescaleTuner ($\tau_\mathrm{damp}$ is a DataVector
149 : * with as many components as the corresponding function of time, thus
150 : * $\tau_\mathrm{exp}$ will also be a DataVector of the same length).
151 : *
152 : * However, this expiration time calculated above is not necessarily the
153 : * expiration that is returned by this function. We group functions of time by
154 : * the `control_system::protocols::Measurement` that their corresponding control
155 : * systems use. This is because one measurement may be used to update many
156 : * functions of time. So the actual expiration time that is used for all the
157 : * functions of time in this group is the *minimum* of the $\tau_\mathrm{exp}$
158 : * of each function of time in the group.
159 : *
160 : * These groups are calculated using `control_system::system_to_combined_names`
161 : * where the list of control systems comes from the passed in `OptionHolder`s.
162 : *
163 : * If the control system isn't active then expiration time is
164 : * `std::numeric_limits<double>::infinity()`, regardless of what the groups'
165 : * expiration time is.
166 : */
167 : template <size_t Dim, typename... OptionHolders>
168 1 : std::unordered_map<std::string, double> initial_expiration_times(
169 : const double initial_time, const int measurements_per_update,
170 : const bool delay_update,
171 : const std::unique_ptr<::DomainCreator<Dim>>& domain_creator,
172 : const std::optional<OptionHolders>&... option_holders) {
173 : std::unordered_map<std::string, double> initial_expiration_times{};
174 :
175 : using control_systems = tmpl::list<typename OptionHolders::control_system...>;
176 :
177 : // First string is name of each control system. Second string is combination
178 : // of control systems with same measurement
179 : std::unordered_map<std::string, std::string> map_of_names =
180 : system_to_combined_names<control_systems>();
181 :
182 : std::unordered_map<std::string, double> combined_expiration_times{};
183 : std::unordered_set<std::string> infinite_expiration_times{};
184 : for (const auto& [control_system_name, combined_name] : map_of_names) {
185 : (void)control_system_name;
186 : if (combined_expiration_times.count(combined_name) != 1) {
187 : combined_expiration_times[combined_name] =
188 : std::numeric_limits<double>::infinity();
189 : }
190 : }
191 :
192 : [[maybe_unused]] const auto combine_expiration_times =
193 : [&initial_time, &measurements_per_update, &delay_update, &domain_creator,
194 : &map_of_names, &combined_expiration_times,
195 : &infinite_expiration_times](const auto& option_holder) {
196 : const std::string& control_system_name = std::decay_t<
197 : decltype(option_holder)>::value_type::control_system::name();
198 : const std::string& combined_name = map_of_names[control_system_name];
199 :
200 : // If not active, leave the expiration time as infinity
201 : if (not option_holder.has_value()) {
202 : infinite_expiration_times.insert(control_system_name);
203 : return;
204 : }
205 :
206 : auto tuner = option_holder->tuner;
207 : Tags::detail::initialize_tuner(make_not_null(&tuner), domain_creator,
208 : initial_time, control_system_name);
209 :
210 : const auto& controller = option_holder->controller;
211 : const DataVector measurement_timescales =
212 : calculate_measurement_timescales(controller, tuner,
213 : measurements_per_update);
214 : const double min_measurement_timescale = min(measurement_timescales);
215 :
216 : const double initial_expiration_time =
217 : function_of_time_initial_expiration_time(
218 : initial_time, DataVector{1, min_measurement_timescale},
219 : measurements_per_update, delay_update);
220 :
221 : combined_expiration_times[combined_name] = std::min(
222 : combined_expiration_times[combined_name], initial_expiration_time);
223 : };
224 :
225 : EXPAND_PACK_LEFT_TO_RIGHT(combine_expiration_times(option_holders));
226 :
227 : // Set all functions of time for a given measurement to the same expiration
228 : // time
229 : for (const auto& [control_system_name, combined_name] : map_of_names) {
230 : initial_expiration_times[control_system_name] =
231 : infinite_expiration_times.contains(control_system_name)
232 : ? std::numeric_limits<double>::infinity()
233 : : combined_expiration_times[combined_name];
234 : }
235 :
236 : return initial_expiration_times;
237 : }
238 : } // namespace control_system
|
