Optimal Experimental Design for Uncertain Systems Based on Coupled Differential Equations

We consider the optimal experimental design problem for an uncertain Kuramoto model, which consists of N interacting oscillators described by coupled ordinary differential equations. The objective is to design experiments that can effectively reduce the uncertainty present in the coupling strengths between the oscillators, thereby minimizing the cost of robust control of the uncertain Kuramoto model. We demonstrate the importance of quantifying the operational impact of the potential experiments in designing optimal experiments.


Introduction
In many real-world applications, we often have to deal with complex systems, for which we do not have complete knowledge. While collecting more data may lead to better system modeling, there exist many scientific applications in which gathering sufficient data for accurate system identification is practically impossible due to the enormous complexity of the system, prohibitively high cost of data acquisition, or both. Relevant examples abound across various domains, including multi-scale climate modeling for long-term prediction, inference of genome-scale regulatory network for predicting effective intervention strategies, characterization of a material system for optimization of targeted functional properties, just to name a few. In such cases, experimental design should target improving one's knowledge of the uncertain system on aspects that critically affect one's operational goals, be they related to control, classification, filtering, or others.
In this work, we consider the problem of optimal experimental design (OED) for an uncertain complex system based on coupled ordinary differential equations (ODEs), called the Kuramoto oscillator model [1,2]. The primary goal is to identify the optimal experiment that is expected to effectively reduce the model uncertainty in such a way that minimizes the cost of controlling the uncertain system. The mean objective cost of uncertainty (MOCU) [3] can be used to quantify the objective-based uncertainty, which then can be used to predict the optimal experiment that maximally reduces the uncertainty pertaining the cost of control [4]. Recently, MOCU-based OED has been applied to a number of applications such as gene regulatory network intervention [5,4], adaptive sequential sampling [6], robust filtering [7], and autonomous materials discovery [8]. For the first time, we tackle the experimental design problem for an uncertain coupled-ODE system. Unlike previous studies [4,5,6,8,7], we do not assume experiments can directly measure the unknown parameters. Instead, we consider realistic experiments, whose outcomes may be used to narrow down the range of unknown parameters.

Kuramoto model of interacting oscillators
We consider the Kuramoto model: where θ i = θ i (t) is the phase of the i-th oscillator with the intrinsic natural frequency ω i , a i,j (= a j,i ) represents the coupling strength constant between the i-th and the j-th oscillators, and N is the total number of oscillators in the model. The system (1) has been first introduced by Yoshiki Kuramoto in [1,2] to describe the phenomenon of collective synchronization observed in the systems of chemical and biological oscillators, in which an ensemble of oscillators spontaneously locks to a common frequency, despite the differences in the natural frequencies of the individual oscillators.
Collective synchronization phenomena have been observed in various biological and chemical systems. These include circadian rhythms, intestinal muscles, menstrual cycles, and fireflies. There are also examples in engineering and physics, for instance, electrical generators, arrays of lasers and superconducting Josephson junction arrays. There have been extensive studies for the Kuramoto model in (1). We refer to [9] and the references therein for derivation of the model and recent developments.

Uncertainty class of Kuramoto models
Now we consider a set of N Kuramoto oscillators, where the natural frequency ω i is known for all oscillators. However the interaction strength a i,j between the i-th and the j-th oscillators is not known with certainty. Instead, we assume that only a lower bound a ℓ i,j and an upper bound a u i,j is known for a i,j . This gives rise to an uncertainty class A of N interacting Kuramoto oscillators, where a = {a i,j }, 1 ≤ i < j ≤ N is an uncertain parameter vector. We assume that a is uniformly distributed in A following the prior distribution below: where

Quantifying the objective cost of uncertainty
Suppose we want to ensure the frequency synchronization of the N interacting Kuramoto oscillators with uncertain interaction strength a ∈ A by adding an additional oscillator which interacts with the N oscillators in the original system. Here the Kuramoto oscillator ensemble ϑ(t) := (θ 1 (t), · · · , θ N (t)) is said to achieve the frequency synchronization asymptotically if it locks to a common frequency such that We assume that the (N + 1)-th oscillator added to the original system for control (i.e., to achieve frequency synchronization) has a known natural frequency ω N +1 and that its interaction strength with the i-th oscillator (part of the original system) is uniform a i,N +1 = a N +1 , ∀i = 1, · · · , N . By selecting a sufficiently large interaction strength a N +1 , we can enforce all N oscillators in the original system to be synchronized with each other in terms of their oscillation frequency (i.e., angular speed). With the introduction of the additional oscillator, now we havė for i = 1, · · · , N . Although a N +1 → ∞ would guarantee synchronization, our goal is to minimize the interaction strength a N +1 as it affects the cost of control, a larger a N +1 resulting in a higher cost for synchronizing the oscillators in the system at hand.
For a given a, we define ξ(a) as the minimum value of the interaction strength a N +1 that guarantees synchronization of all oscillators. As a N +1 = ξ(a) would be optimal for a specific a, we call it ξ(a) the optimal interaction strength.
In the presence of uncertainty, we are unable to identify ξ(a) since a is unknown. Instead, we desire an optimal robust interaction strength ξ * (A) such that Note that it is robust because a N +1 = ξ * (A) guarantees synchronization for any a ∈ A. It is optimal because it is the smallest such value. As we can see from (3), a N +1 increases due to the uncertainty, which forces us to choose a larger interaction strength than might be actually needed for synchronization. The expected increase of this differential cost can be measured by computing the expected value of the cost increase based on p(a), which governs the distribution of a within the uncertainty class A. This average differential cost M (A) is referred to as the mean objective cost of uncertainty (MOCU) [3], and it quantifies the impact that the model uncertainty has on the operational objective.

Optimal Experimental Design
Suppose we want to perform additional experiments to reduce the uncertainty class. In general, the outcome of an experiment may reduce the uncertainty class, which may help us predict a better robust controller-in this case, the (N + 1)-th oscillator with a smaller interaction strength a N +1 = ξ * (A) that ensures the synchronization of all oscillators despite the uncertainty in a. A practical question arises naturally: among the possible experiments, how can we select the optimal experiment? We address this question in what follows.

Experimental design space
We restrict our experimental design space to experiments that test pairwise synchronization between oscillators. Suppose we choose the oscillator pair (i, j) for our experiment, where we initialize the angles to θ i (0) = θ j (0) and observe whether the two oscillators will become eventually synchronized in the absence of any influence from all other oscillators. We define a binary random variable B i,j for the experimental outcome, where B i,j = 1 corresponds to the eventual synchronization of the oscillator pair, while B i,j = 0 corresponds to the opposite.
If the interaction strength a i,j between the oscillators i and j were known, the outcome B i,j would be known with certainty. In fact, Theorem 1 below shows that the oscillator pair will be synchronized if and only if |ω i − ω j | ≤ K, where K = 2a i,j in this case.
However, the uncertainty regarding a i,j renders B i,j a random variable whose outcome is unknown before performing the pairwise synchronization experiment described above. Suppose our experiment results in the eventual synchronization of the two oscillators (B i,j = 1). Based on Theorem 1, the following inequality must hold: This experimental outcome allows us to update the lower bound of a i,j from a ℓ i,j toã ℓ i,j = max( 1 2 |ω i − ω j |, a ℓ i,j ). On the other hand, if the oscillators do not get synchronized (B i,j = 0), we have which allows us to update the upper bound of a i,j from a u i,j toã u i,j = min( 1 2 |ω i −ω j |, a u i,j ). In either case, the pairwise experiment can potentially reduce the uncertainty regarding a i,j , thereby shrinking the uncertainty class A.

Selecting the optimal experiment
Knowing that the aforementioned pairwise experiments can potentially reduce the uncertainty class, how should we prioritize the experiments to select the optimal one? The MOCU framework can be used to predict the optimal experiment that is expected to maximally reduce the uncertainty [4,5,6,7] in such a way that minimizes the cost of uncertainty, namely, the expected cost increase for controlling (i.e., synchronizing) the N Kuramoto oscillators due to the uncertain interaction strength. More specifically, for every experiment in the experimental design space, we first compute the expected remaining MOCU after performing the given experiment. Based on these results, we can prioritize the experiments and select the one that is expected to minimize the MOCU that remains after carrying out the experiment.
For convenience, let us denote the uncertainty class A reduced based on the experimental outcome B i,j as A|B i,j . Then the expected remaining MOCU for the synchronization experiment of the oscillator pair (i, j) can be computed by Based on the prior p(a) in (2), we can compute the probabilities for the possible experimental outcomes as follows to ensure thatã i,j ∈ [a ℓ i,j , a u i,j ]. The optimal oscillator pair, the outcome of whose pairwise synchronization experiment is expected to most effectively improve the control performance among the N 2 pairs can be predicted by (i * , j * ) = arg min 1≤i<j≤N R(i, j). (10)

Simulation Results
In this section, we present numerical experiments to demonstrate our proposed OED method described in Section 3.2.
A classical fourth-order Runge-Kutta method is used to compute the Kuramoto model in (1) for 0 ≤ t ≤ T for T = 5 with time discretization ∆t = 1/160. For the sake of simplicity, the initial conditions are set to θ i = 0, 1 ≤ i ≤ N . As the system is regular enough, the Runge-Kutta solvers provide reasonably accurate numerical solutions. For instance, the relative L 2 -error of a five-oscillator model (N = 5) is 10 −10 ∼ 10 −9 at t = T . To examine synchronization of the model, difference in θ i is measured by ∆θ i (t) := θ i (t + ∆t) − θ i (t). If the model is synchronized, the maximum difference is sufficiently small for t > T /2 such that with tol = 10 −4 . Due to the immense uncertainty in parameters, the sample size K for computing the expectation in (4) should be sufficiently large (e.g. K ≥ 20, 000) to obtain reliable experimental results. In this regard, massive computing power is highly desired, and we adopt GPU parallel computing, PyCUDA [10], using Nvidia GTX 1080-Ti.
For the additional oscillator used for control, we simply chose w 6 = mean 1≤i≤5 w i . The upper and lower bounds of interaction strength were chosen by  where d i,j is a correction constant. If d i,j ≡ 1 for all i, j, the system is already synchronized in general as all entries of the interaction strength a i,j are large enough. Hence, we introduced the correction parameter d i,j such that where 0.3 ≤ d * i,j ≤ 0.5, and I 1 and I 2 are partition of the set of indices I = {(i, j) ∈ N : 1 ≤ i, j ≤ N, i < j}. Here, the set I 1 , I 2 and the corresponding quantity d * i,j were empirically determined. To compare numerical results, we carried out a sequence of experiments based on three different experiment selection strategies: MOCU-based, random selection, and entropy-based. In the MOCU-based selection strategy, the pairwise experiment with the smallest expected remaining MOCU in (8) was chosen, and the corresponding entry of a u i,j or a ℓ i,j was updated based on the experimental outcome. More precisely, from the result of the pairwise experiment between the i-th and the j-th oscillators, if the oscillators were synchronized, then the lower bound a ℓ i,j was updated toã i,j defined in (9). Otherwise, the upper bound a u i,j was updated toã i,j . In the entropy-based approach, we selected the pairwise experiment with the largest value of a u i,j − a ℓ i,j , and the corresponding entry of a u i,j or a ℓ i,j was updated based on the outcome in the same fashion. In the random approach, we randomly chose one of the N 2 possible experiments and updated the corresponding entry of a u i,j or a ℓ i,j based on the experimental outcome. Figure 1 shows that the MOCU-based experiment selection strategy leads to the most efficient updates among the three strategies. Especially, it identified effective experiments early on, nearly reaching the minimum attainable uncertainty level in just 3 updates.
We set the natural frequency of the additional oscillator to  To start with a non-synchronized model, we selected sufficiently large natural frequencies for w 8 and w 9 . The upper and lower bounds of the interaction strengths were set to a u i,j = 1.03d i,j where d i,j is the correction parameter and that was empirically determined as in the previous example. Figure 2 shows the simulation results, which clearly demonstrate that the MOCU-based experimental design results in the most efficient updates among all three methods. Here we considered up to 18 updates (out of 36 experiments in total). As shown in Figure 2, the MOCU-based OED was able to drastically reduce uncertainty with a single update, and it reached the minimum uncertainty level just in 9 updates.

Concluding Remarks
As shown in our results, designing effective experiments for complex uncertain systems requires quantifying the state of our current knowledge of the system and measuring the impact of the remaining uncertainty on the operator performance. It is clear that we cannot expect system-agnostic black-box optimization schemes to perform well. Furthermore, experiments that aim to enhance our knowledge regarding parameters with the largest uncertainties do not necessarily help, as they may not be pertinent to the operator performance. Our work shows that the MOCUbased OED framework can effectively prioritize the experiments in the design space by quantifying the impact of their potential outcomes on the operational goal using scientific knowledge as in Theorem 1.