Optimal Sensor Scheduling and Selection for Continuous-Discrete Kalman Filtering with Auxiliary Dynamics

Q1: Besides the auxiliary dynamics, our method considers the continuous-discrete setup and the flexibility of not requiring an upper bound on the number of measurements. The method by [Le Ny et al., 2009] is for a continuous-continuous setup, while the method in [Marelli et al., 2019] is for a discrete-discrete setup. We will make sure to clarify this point further in the paper.

Q2: Scalability: For each new sensor we consider, we will have an intensity rate, which we will need to optimize for. This is typical for optimization-based scheduling procedures. Considering efficient optimization schemes for large-scale problems will require investigating techniques such as distributed optimization, parralization, and sparsity considerations. Another challenge for scalability is the dimension of the states of the SSM we want to estimate. If we have $n_x$ states, then the covariance matrix of the CD-KF will be in the order of $n_x^2$. However, for large dimensional SSM, there exist research for efficient approximations and methods to deal with the covariance matrix problem. One solution, for example, is to consider a diagonal approximation for the covariance matrix or a low-rank approximation (see [4]). These approximations can be easily integrated into our framework. We will include a remark discussing this in the final paper.

Concavity Assumption: We emphasize that even under the concavity assumption, our approach spans a wide range of dynamical systems (it covers all linear parameter varying systems). This work is the first to explore sensor scheduling with auxiliary dynamics. Additionally, our method can be implemented within a Receding Horizon (RH) framework (refer to Reviewer nw2r's response for a detailed description), where each optimization iteration is addressed by either linearizing the non-convex/non-concave auxiliary state or employing a convex/concave approximation. If we use nonconvex/nonconcave dynamics for $\xi_p$, then we will lose the theoretical guarantees of the method. However, we may still obtain satisfactory results depending on how good the approximation $\frac{d{\mathbb{E}[\xi_p]}}{dt}\approx f_p(\mathbb{E}[\xi],u,t)+\sum^{N_s}_{s=1} \lambda_s(t)g_s(\mathbb{E}[\xi],u,t))$ (or an upper bound approximate). This approximation is similar to the approximation used in the EKF. We conducted an experiment with non-concave dynamics for the sensor degradation states $\zeta_1$ and $\zeta_2$ of the example in the paper. The experiment demonstrated that our method remains effective (link: postimg.cc/jWcFJ4f3) (as allowed by ICML). We will include this discussion with the experimental results in the paper.

Q3: The KF is usually chosen in the literature of sensor scheduling instead of the PF for the fact that the covariance matrix dynamics are independent of the actual measurements and because it scales better with large dynamical systems. However, we can still extend our approach with the EKF by utilizing an RH setup where, for each short horizon, we use the linearized dynamics around the current estimate. We will include a remark about this point in the paper.

Q4: We apologize to the reviewer as we did not understand what the reviewer meant by "high-dimensionality of the problem" in connection with our paper. If this refers to scalability, then we have addressed it above. The nonlinear dynamics of the problem were also addressed above.

Comparisons: We acknowledge and agree with the reviewer's suggestion to provide comparisons. We have conducted comparisons with heuristic approaches for scheduling measurements (greedy approach and random sampling of measurement times) that we will provide in the paper (Table link: postimg.cc/qz2Q3GVc). The results suggest that our method ("Optimized") outperforms the greedy and random scheduling approaches for our example. To assess the deterministic scheduling computed based on the optimized measurement rates (denoted "Optimized" ), we compared it with a method (denoted as "M-Optimized" in the table) that is based on sampling $M_c$ realizations of measurement times according to the corresponding Poisson process with the optimized rates. Afterwards, we pick the measurement times corresponding to the realization with the minimal cost. The results show that our deterministic quantization provides similar results to "M-Optimized" without having to sample multiple realizations, which can be computationally expensive and unrealistic since we do not have the real measurements to compute the cost for each realization. To the authors' best knowledge, no reinforcement learning methods applicable to this setup have been found for comparison (i.e., multiple sensor scheduling in a continuous-discrete setting that does not require training on pre-obtained data with uniform sampling).

[4] Chang, Peter G., et al. "Low-rank extended Kalman filtering for online learning of neural networks from streaming data." Conference on Lifelong Learning Agents. PMLR, 2023.

Unknown dynamics:
Our formulation as an Optimal Control Problem (OCP) with a differentiable cost function and constraints opens avenues for extension to uncertain dynamics using robust/stochastic OCP methods [2,3]. For completely unknown dynamics, planning is challenging since we must schedule measurements based on dynamics we do not yet know. We believe our approach paves the way for future research in this direction. One approach to handling unknown parameters is a Receding Horizon (RH) setup. We solve a finite-horizon OCP, assuming fixed system dynamics and parameters over the prediction horizon based on current estimates, then apply the control until the first measurement. The subsequent measurement updates the parameters, and the OCP is resolved for the next horizon. This iterative process adapts as more information becomes available. In the revised version, we will include successful results using RH for a moving target (link: postimg.cc/S2t1gQDW (as allowed by ICML)).

Performance characteristics:
Performance, particularly regarding discretization steps, depends on the auxiliary dynamics (e.g., stiffness, fast/slow dynamics, stability) and the KF’s covariance dynamics. Different OCP methods offer trade-offs between accuracy and computational performance (see Appendix D). Our work employed Euler discretization with variable steps while keeping a fixed number of discretization points. In the revision, we will include a figure demonstrating the trade-off between computation time and discretization points of a specific example focused on the KF dynamics (link: postimg.cc/ykL46hPv).

Penalizing auxiliary states:
Often, the auxiliary state carries physical meaning, making its penalization more intuitive. For instance, in our examples, the energy state depends on measurement rates, velocities (actions), and the robot’s position. Penalizing measurement rates is less straightforward than penalizing energy consumption when energy is the limiting factor.

Constraints handling:
Both direct and indirect methods for constrained OCP lack exact guarantees of constraint satisfaction due to numerical errors inherent in the OCP solution and ODE integration. Some methods (e.g., direct collocation) can achieve high accuracy but at increased computational cost (See Appendix D). Ultimately, the choice of method and integration scheme depends on the specific dynamics, much like the selection of ODE solvers. We will include this important remark in the paper.

For gradient computations and optimization in the example, we used an interior-point method with JuMP—a Julia package that employs automatic differentiation to compute the gradient and Hessian of the Lagrangian.

Regarding $\xi$:
We acknowledge the reviewer's concern. To clarify, $\xi^*$ from Proposition 5.1 can be any curve (ensuring (10) is well-defined, though not stated explicitly), so $\hat{\xi}$ and $\xi^*$ need not be close. Rather, $\xi^*$ serves as a placeholder, allowing substitution with $\hat{\xi}$. Proposition 5.1 states that for any curve $\xi^*$, the mean covariance $\bar{\Sigma}(\xi^*)$ is bounded by $\hat{\Sigma}(\xi^*)$. This result is applied by substituting $\xi^*$ with the curve $\hat{\xi}$ obtained from (12) and (6b) with initial condition $\xi_0$, which, per Proposition 5.2, bounds the mean curve $\bar{\xi}$ of $\xi=(\xi_p,\xi_u)$ from the SDE (10) and ODE (6b). The deterministic quantities $\bar{\xi}$ (and $\bar{\Sigma}$) serve as computationally tangible approximations of the stochastic $\xi$ (and $\Sigma$) (approximating mean behaviour). Alternatively, sampling methods can be used to compute statistical representations for the trajectories in the OCP. However, this will introduce non-differentiability and will be computationally intensive. We will note this important remark in the revision. We apologize for the initial lack of clarity.

Experimental design questions:

1. As mentioned, the RH approach can handle moving targets (we have implemented an example for it). Nonetheless, many scenarios involve a fixed process location (e.g., gas leak, specific object temperature) or applications independent of process location, such as underwater measurements with sensor biofouling. A detailed example of the latter will be provided in the paper.
2. The process assumes a SSM representation of the Matern kernel, with output $x_p$. The parameters $A$ and $\sigma$ in equation (13.b) are based on this representation. We will update the figure to a more user-friendly version. See also the reply to reviewer 6UyF.

[2] Leeman, Antoine P., et al. "Robust optimal control for nonlinear systems with parametric uncertainties via system level synthesis." the 62nd IEEE Conference on Decision and Control (CDC). IEEE, 2023.

[3] Bemporad, Alberto, and Manfred Morari. "Robust model predictive control: A survey." Robustness in identification and control. London: Springer London, 2007. 207-226.

Q1: In the last paragraph of Section 8.1, we leverage the fact that Gaussian process (GP) regression (with many common stationary covariance kernels) is equivalent to Kalman smoothing of a specific linear state-space model (see ref. Sarkka & Hartikainen, 2012 in the paper). This equivalence allows us to perform GP regression efficiently on the temporal process $x_p$ via Kalman smoothing. This example demonstrates how our approach can be used for sensor scheduling in connection with GP regression. We will provide an appendix section in the paper clarifying this point.

Q2: Our work focuses on Bayesian inference in SSMs, a topic that has been explored in previous papers at ICML and other ML-related conferences (e.g., [1]). Additionally, our work can be applied to sensor scheduling in GP regression under dynamic environments. We believe that GP regression continues to be a subject of significant interest in the machine learning community.

Q3: We thank the reviewer for spotting the sign error. The right-hand side of the equation should be $c_e \exp\bigl(r_e \|p_r - p_e\|^2\bigr) - c_u v - c_u \omega - \sum_{s=1}^2 \sum_{i=1}^{N_s} c_s \,\delta_{t^s_i}$.

Q4: We see the reviewer's point. We wrote it with the intention of providing more clarity for the reader on the total measurement process. But as the reviewer points out, it is not used; it may therefore contribute to more confusion for the reader than clarity. We will remove it from the paper.

Q5: We apologize for the ambiguity of the figure. This point has also been mentioned by the second reviewer. We will fix the figure to be more reader-friendly for the final submission. The symbol $x_p$ represents the process we want to measure. The filter estimate is now denoted as $\hat{x}_p$, and the RTS smoother estimate is denoted as $\hat{x}^s_p$ (representing the GP regression output). We have adjusted the legends and description of the figure with the new notation for the filter estimate and the smoother estimate (link: postimg.cc/hJZ0ByTn (as allowed by ICML)).

[1]: Duran-Martin, Gerardo, et al. "Outlier-robust Kalman Filtering through Generalised Bayes." International Conference on Machine Learning. PMLR, 2024.