## 1 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.

## 2 Uncertain Kuramoto Model

### 2.1 Kuramoto model of interacting oscillators

We consider the Kuramoto model:

(1) |

where is the phase of the -th oscillator with the intrinsic natural frequency , represents the coupling strength constant between the -th and the -th oscillators, and 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.

### 2.2 Uncertainty class of Kuramoto models

Now we consider a set of Kuramoto oscillators, where the natural frequency is known for all oscillators. However the interaction strength between the -th and the -th oscillators is not known with certainty. Instead, we assume that only a lower bound and an upper bound is known for . This gives rise to an uncertainty class of interacting Kuramoto oscillators, where

is an uncertain parameter vector. We assume that

is uniformly distributed in

following the prior distribution below:(2) |

where .

### 2.3 Quantifying the objective cost of uncertainty

Suppose we want to ensure the frequency synchronization of the interacting Kuramoto oscillators with uncertain interaction strength by adding an additional oscillator which interacts with the oscillators in the original system.
Here the Kuramoto oscillator ensemble
is said to achieve *the frequency synchronization asymptotically* if it locks to a common frequency such that

We assume that the -th oscillator added to the original system for control (i.e., to achieve frequency synchronization) has a known natural frequency and that its interaction strength with the -th oscillator (part of the original system) is uniform , . By selecting a sufficiently large interaction strength , we can enforce all 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 have

for . Although would guarantee synchronization, our goal is to minimize the interaction strength as it affects the cost of control, a larger resulting in a higher cost for synchronizing the oscillators in the system at hand.

For a given , we define as the minimum value of the interaction strength that guarantees synchronization of all oscillators. As would be optimal for a specific , we call it the optimal interaction strength. In the presence of uncertainty, we are unable to identify since is unknown. Instead, we desire an optimal robust interaction strength such that

(3) |

Note that it is robust because guarantees synchronization for any . It is optimal because it is the smallest such value. As we can see from (3), 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

(4) |

based on , which governs the distribution of within the uncertainty class . This average differential cost 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.

## 3 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 -th oscillator with a smaller interaction strength that ensures the synchronization of all oscillators despite the uncertainty in . A practical question arises naturally: among the possible experiments, how can we select the optimal experiment? We address this question in what follows.

### 3.1 Experimental design space

We restrict our experimental design space to experiments that test pairwise synchronization between oscillators. Suppose we choose the oscillator pair for our experiment, where we initialize the angles to

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

for the experimental outcome, where corresponds to the eventual synchronization of the oscillator pair, while corresponds to the opposite.If the interaction strength between the oscillators and were known, the outcome would be known with certainty. In fact, Theorem 1 below shows that the oscillator pair will be synchronized if and only if , where in this case.

###### Theorem 1.

Consider the Kuramoto model of two-oscillators:

(5) |

with the initial angles . Then, for any solutions and to (5), there holds as if and only if .

###### Proof.

First assume that . By symmetry, we may assume without loss of generality, that and that . By subtracting the equations of (5), one has

(6) |

where . Since , there is a number such that . We can easily check that , indicating that is a unique (mod ) stable critical point of (6). This implies that (mod ) as unless , in turn . When , which is an unstable critical point, we have for all . These prove our assertion. On the other hand, if , one has from (6) that

(7) | |||||

This completes the proof. ∎

However, the uncertainty regarding renders 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 (). Based on Theorem 1, the following inequality must hold:

This experimental outcome allows us to update the lower bound of from to . On the other hand, if the oscillators do not get synchronized (), we have

which allows us to update the upper bound of from to . In either case, the pairwise experiment can potentially reduce the uncertainty regarding , thereby shrinking the uncertainty class .

### 3.2 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 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 reduced based on the experimental outcome as . Then the expected remaining MOCU for the synchronization experiment of the oscillator pair can be computed by

(8) | |||||

Based on the prior in (2

), we can compute the probabilities for the possible experimental outcomes as follows

where we define

(9) |

to ensure that . The optimal oscillator pair, the outcome of whose pairwise synchronization experiment is expected to most effectively improve the control performance among the pairs can be predicted by

(10) |

## 4 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 for with time discretization . For the sake of simplicity, the initial conditions are set to , . As the system is regular enough, the Runge-Kutta solvers provide reasonably accurate numerical solutions. For instance, the relative -error of a five-oscillator model () is at . To examine synchronization of the model, difference in is measured by

If the model is synchronized, the maximum difference is sufficiently small for such that

with . Due to the immense uncertainty in parameters, the sample size for computing the expectation in (4) should be sufficiently large (e.g. ) 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.

As a paradigm example, we first implemented a 5-oscillator Kuramoto model. The following natural frequencies were used for the experiment in Figure 1:

For the additional oscillator used for control, we simply chose

The upper and lower bounds of interaction strength were chosen by

where is a correction constant. If for all , the system is already synchronized in general as all entries of the interaction strength are large enough. Hence, we introduced the correction parameter such that

where , and and are partition of the set of indices Here, the set and the corresponding quantity 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 or was updated based on the experimental outcome. More precisely, from the result of the pairwise experiment between the -th and the -th oscillators, if the oscillators were synchronized, then the lower bound was updated to defined in (9). Otherwise, the upper bound was updated to . In the entropy-based approach, we selected the pairwise experiment with the largest value of , and the corresponding entry of or was updated based on the outcome in the same fashion. In the random approach, we randomly chose one of the possible experiments and updated the corresponding entry of or 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.

Next, we conducted numerical simulations with a larger number of oscillators, , in which case the number of possible experiments increases to . We used the following natural frequencies:

We set the natural frequency of the additional oscillator to

To start with a non-synchronized model, we selected sufficiently large natural frequencies for and . The upper and lower bounds of the interaction strengths were set to

where 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.

## 5 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 MOCU-based 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.

## References

- [1] Yoshiki Kuramoto. Self-entrainment of a population of coupled non-linear oscillators. In International symposium on mathematical problems in theoretical physics, pages 420–422. Springer, 1975.
- [2] Yoshiki Kuramoto. Chemical Oscillations, Waves, and Turbulence. Courier Corporation, 2003.
- [3] Byung-Jun Yoon, Xiaoning Qian, and Edward R. Dougherty. Quantifying the objective cost of uncertainty in complex dynamical systems. IEEE Transactions on Signal Processing, 61(9):2256–2266, 2013.
- [4] Roozbeh Dehghannasiri, Byung-Jun Yoon, and Edward R Dougherty. Optimal experimental design for gene regulatory networks in the presence of uncertainty. IEEE/ACM Transactions on Computational Biology and Bioinformatics (TCBB), 12(4):938–950, 2015.
- [5] Roozbeh Dehghannasiri, Byung-Jun Yoon, and Edward R. Dougherty. Efficient experimental design for uncertainty reduction in gene regulatory networks. BMC Bioinformatics, 16(13):S2, 2015.
- [6] Ariana Broumand, Mohammad Shahrokh Esfahani, Byung-Jun Yoon, and Edward R Dougherty. Discrete optimal bayesian classification with error-conditioned sequential sampling. Pattern Recognition, 48(11):3766–3782, 2015.
- [7] G. Zhao, X. Qian, B.-J. Yoon, F. J. Alexander, and E. R. Dougherty. Model-based robust filtering and experimental design for stochastic differential equation systems. IEEE Transactions on Signal Processing, 68:3849–3859, 2020.
- [8] Anjana Talapatra, Shahin Boluki, Thien Duong, Xiaoning Qian, Edward Dougherty, and Raymundo Arróyave. Autonomous efficient experiment design for materials discovery with Bayesian model averaging. Physical Review Materials, 2(11):113803, 2018.
- [9] Steven Strogatz. From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators. Physica D, 143(1):1–20, 2000.
- [10] Andreas Klöckner, Nicolas Pinto, Yunsup Lee, B. Catanzaro, Paul Ivanov, and Ahmed Fasih. PyCUDA and PyOpenCL: A Scripting-Based Approach to GPU Run-Time Code Generation. Parallel Computing, 38(3):157–174, 2012.

Comments

There are no comments yet.