Bifurcation-diagram reconstruction estimates various attractors of a system without observing all of them but only from observing several attractors with different parameter values. Therefore, the bifurcation-diagram reconstruction can be used to investigate how attractors change with the parameter values, especially for real-world engineering and physical systems for which only a limited number of attractors can be observed. Although bifurcation diagrams of various systems have been reconstructed from time-series data generated in numerical experiments, the systems that have been targeted for reconstructing bifurcation diagrams from time series measured from physical phenomena so far have only been continuous-time dynamical systems. In this paper, we reconstruct bifurcation diagrams only from time-series data generated by electronic circuits in discrete-time dynamical systems with different parameter values. The generated time-series datasets are perturbed by dynamical noise and contaminated by observational noise. To reconstruct the bifurcation diagrams only from the time-series datasets, we use an extreme learning machine as a time-series predictor because it has a good generalization property. Hereby, we expect that the bifurcation-diagram reconstruction with the extreme learning machine is robust against dynamical noise and observational noise. For quantitatively verifying the robustness, the Lyapunov exponents of the reconstructed bifurcation diagrams are compared with those of the bifurcation diagrams generated in numerical experiments and by the electronic circuits.

Many real-world systems exhibit attractors with different parameter values of the system. We hope to estimate such attractors, although the system itself is usually unknown. By using bifurcation-diagram reconstruction, we can estimate these attractors if several time-series data can be generated by the target system at different parameter values. The purpose of this paper is to show that the bifurcation-diagram reconstruction method that Itoh *et al.* have already proposed^{1} is robust against dynamical and observational noise because time-series data measured from physical phenomena are always influenced by dynamical and observational noise. We therefore demonstrate the results of reconstructing bifurcation diagrams only from time-series data generated by electronic circuits in discrete-time dynamical systems. Herein, we expect that a filtering effect is obtained by using an extreme learning machine as a time-series predictor, which has a good generalization property.

## I. INTRODUCTION

Attractors in real-world systems can be measured as time-series data. If the system is exhibiting deterministic chaos, for example, then the data can be predicted in the short term by using a predictor that is trained to model the measured time-series data, although a target dynamical system itself is usually unknown. If the target system has bifurcation parameters and we wish to estimate its bifurcation diagram (BD), then doing so requires a large number of data that are measured while changing the parameter values that can be controlled.^{2} However, if cost and time are limited, then a sufficient amount of data often cannot be obtained. To address this problem, in 1994, Tokunaga *et al.*^{3} proposed a method of BD reconstruction, whereby the BD of the target system is estimated only from a few time-series datasets measured at different parameter values; by numerical experiments, Tokunaga *et al.* reconstructed the BDs of the Hénon map and coupled logistic/delayed-logistic maps only from time-series data without dynamical and observational noise. Tokuda *et al.*^{4} then reconstructed the BD of the Rössler equations only from time-series data with observational noise, and in 2000, Bagarinao *et al.*^{5} reconstructed the BDs of a cubic map and the FitzHugh–Nagumo equations only from time-series data with dynamical noise. In 2001, Small *et al.*^{6} estimated the BD of the Rössler equations only from a time-series dataset with observational noise. Here, the difference between their method and other methods is the time-series data used for modeling, because they used a time-series dataset generated by a system with a changing parameter. After that, Small *et al.*^{7} estimated the BD during the onset of human ventricular fibrillation.

Several research groups reconstructed BDs of several nonlinear systems from time-series data by numerical experiments. However, behavior of a real-world system cannot be faithfully reproduced on digital computers, especially when a chaotic system is targeted. For example, Shi reported that attractors of an extended logistic map system are changed depending on computational precision.^{8} Therefore, the effectiveness of the BD reconstruction against a real-world system needs to be evaluated by real data. Obviously, it is necessary to consider dynamical and observational noise in real-world problems. For these reasons, analog electronic circuits that have high reproducibility of mathematical experiments and can easily control experimental conditions are often used as experimental systems that generate time-series data based on analog computation. Moreover, analog electronic circuits operated by using physical characteristics like transistor characteristics are incompletely known physical systems, even if these are artificially designed. The characteristics cannot be fully controlled because an impurity distribution in the semiconductor causes fluctuations in physical properties like the threshold voltage of transistors. On the other hand, digital computers as well as application-specific integrated digital circuits are complete known systems, because they are not affected by physical properties due to binarization and error correction. Therefore, the analog electronic circuits are suitable to verify effectiveness of a method against unknown physical systems. In 2004, Langer and Parlitz^{9} reconstructed the BD of a Colpitts oscillator from time-series data generated by electronic circuit, thereby showing that BDs can be reconstructed for real-world systems subjected to intrinsic noise. Langer and Parlitz used not only time-series datasets but also the parameter values at which those datasets were generated. However, it is not always possible to know the parameter values at which time-series data were generated. Nevertheless, in 2018, Itoh and Adachi reconstructed the BD of the Rössler equations only from time-series data generated by electronic circuit.^{10}

In this paper, we verify that the BD reconstruction method that Itoh *et al.* have already proposed^{1} is robust against dynamical and observational noise using time-series data generated by electronic circuits that have high reproducibility of dynamical systems. Therefore, we use the electronic circuits that can realize nonlinear functions similar to nonlinear voltage waveforms inputted by an arbitrary waveform generator,^{11} unlike general nonlinear analog electronic circuits that generate chaos using the nonlinear characteristics of devices. However, dynamical systems that can be realized with this electronic circuits are discrete-time dynamical systems. In addition, as far as we know, the systems that have been targeted for reconstructing BDs from time-series data generated by electronic circuits have only been continuous-time dynamical systems. For these reasons, we target discrete-time dynamical systems in this paper. For clearly verifying the BD reconstructed from time series influenced by dynamical and observational noise, we use the logistic and sine maps as target systems that are simple discrete-time dynamical systems. This is caused by the BD reconstruction that is not always successful even when the target system is simple as shown by Itoh and Adachi.^{12} In the paper, they have shown success rates for reconstructing the BDs of logistic and sine maps.

In general, time-series data generated by electronic circuits are contaminated by observational noise resulting from analog-to-digital conversion and perturbed by dynamical noise peculiar to the circuits. For qualitatively evaluating the influence of dynamical noise, we compare between the BD of logistic map generated by real data from the electronic circuit and that with dynamical noise generated by a numerical experiment. In addition, the Lyapunov exponents are compared between the BDs generated by data from the electronic circuits and those by numerical experiments. We use two methods to estimate the Lyapunov exponents for the reconstructed BD, namely, with a high degree of accuracy from the Jacobian matrix^{1} and from time-series data^{13} for the BD generated by the electronic circuits.

Herein, we use an extreme learning machine (ELM)^{14} as a time-series predictor for reconstructing BDs. Tokunaga *et al.*^{3} used a neural network that was trained via the gradient decent learning of synaptic weights and biases, but a neural network with iterative learning has shortcomings such as sensitivity to dynamical and observational noise. Because our generated time-series data are influenced by dynamical and observational noise, we use an ELM, which has a good generalization property. In addition, Itoh and Adachi have shown previously that an ELM is suitable for the BD reconstruction.^{1,12,15–17}

The rest of this paper is organized as follows. In Sec. II, we explain the algorithm for reconstructing BDs with Lyapunov exponents using the ELM.^{1} In Sec. III, we describe the electronic-circuit realization in discrete-time dynamical systems, and in Secs. IV and V, we present our experimental results for the logistic map and the sine map, respectively. Finally, we draw conclusions in Sec. VI.

## II. ALGORITHM FOR RECONSTRUCTING BIFURCATION DIAGRAMS WITH LYAPUNOV EXPONENTS USING AN EXTREME LEARNING MACHINE

Because we use an ELM as a time-series predictor for reconstructing BDs, we begin by explaining the ELM. We then explain BD reconstruction briefly. The following explanations assume that the ELM has one input neuron and one output neuron because the target systems herein are the logistic and sine maps, which are both one-dimensional maps.

### A. Extreme learning machine

In 2004, Huang *et al.*^{14} proposed an ELM, which is a feed-forward neural network consisting of three layers. Figure 1 shows the structure of the ELM used in the present study. The training targets in the ELM are the synaptic weights of only the output neuron, which are trained by linear regression. The training of the ELM, therefore, avoids local minima and is extremely fast. We use an ELM in this study because of its good generalization property, which is required in BD reconstruction to estimate various attractors.

The output $hi\u2208R$ of hidden neuron $i$ is described by

where $wi\u2208R$ and $bi\u2208R$ are the synaptic weight and bias, respectively, of hidden neuron $i$, $y(t)\u2208R$ is the input value, $t$ is the discrete-time, and $s(\u22c5)$ is a standard sigmoidal function. The value of the output neuron is given by

where $\beta \u2208RV$ is the synaptic weight vector of the output neuron and $h(t)=[h1(t),h2(t),\u2026,hV(t)]T$ is the output vector of the hidden neurons.

The synaptic weights of the output neuron are calculated by

where $H\u2020$ is the Moore–Penrose generalized inverse of the output matrix $H=[h(1),h(2),\u2026,h(L)]T$ of the hidden neurons, $d=[d(1),d(2),\u2026,d(L)]T$ is the desired output vector, and $L$ is the number of training patterns.

### B. Bifurcation-diagram reconstruction based on the method proposed by Tokunaga *et al.*

In 1994, Tokunaga *et al.*^{3} proposed the reconstruction of BDs that estimate attractors in parameter space when parameter values change. The BD reconstruction requires several time-series datasets generated from a system with different parameter values. There are three steps to our BD reconstruction process. First, we train the synaptic weights of the output neuron in the ELM to model the given time-series datasets. Next, we apply principal component analysis (PCA) to the trained synaptic weights to estimate parameter space corresponding to the original parameter space. Finally, we plot the BD using the PCA results.

Herein, we define a nonlinear map for the input–output relations of the target system as

where $x(t)\u2208R$ and $x(t+1)\u2208R$ are the state values of the target system of the one-dimensional map at $t$ and $t+1$, respectively, $pn\u2208RC$ is a parameter vector, and $f(\u22c5,\u22c5)$ is a nonlinear map. Here, $C$ is the dimension of the parameter space. The time-series datasets generated by the system with parameter vectors $p1,p2,\u2026,pP$ are defined as $S1,S2,\u2026,SP$, respectively.

The first step of the BD reconstruction is to generate time-series predictors

to model the time-series datasets $Sn(n=1,2,\u2026,P)$, where $g(\u22c5,\u22c5)$ is a nonlinear function. Here, the trained synaptic weight vectors $\beta 1,\beta 2,\u2026,\beta P$ are trained to model the time-series datasets $S1,S2,\u2026,SP$, respectively.

The second step of the BD reconstruction is to apply PCA to the trained synaptic weight vectors. To do so, deviation vectors of the synaptic weight vectors are calculated by

where $\beta 0\u2208RV$ is the mean vector of the synaptic weight vectors. We then apply PCA to the deviation vectors $\delta \beta 1,\delta \beta 2,\u2026,\delta \beta P$ and obtain their eigenvalues and eigenvectors. Sorted in the descending order, the eigenvalues are $\varphi 1\u2265\varphi 2\u2265\cdots \u2265\varphi V$, and the eigenvectors $vi\u2208RV(i=1,2,\u2026,V)$ correspond to $\varphi i(i=1,2,\u2026,V)$. The deviation vectors are represented by

where $\gamma \xafn\u2208RV$ is the estimated vector. Herein, as in our previous study,^{1} we assume that the eigenvalues $\varphi 1,\varphi 2,\u2026,\varphi E$ contain enough information when the $E$th cumulative contribution ratio exceeds $80%$.

The third step of the BD reconstruction is to plot the BD. The deviation vectors can be approximated as

where $\gamma n\u2208RE(n=1,2,\u2026,P)$ are the estimated parameter vectors. The estimated parameter vector corresponding to the time-series dataset $Sn$ is calculated by

The sequence of estimated parameter vectors $\gamma 1\u2192\gamma 2\u2192\cdots \u2192\gamma P$ corresponds to the sequence of parameter vectors $p1\u2192p2\u2192\cdots \u2192pP$ of the target system. Here, the sequence of estimated parameter vectors is referred to as a bifurcation locus, and the sequence of parameter vectors of the target system is referred to as a bifurcation path. It is via the relationship between the bifurcation path and locus that the estimated parameter space is viewed as corresponding to the parameter space of the target system. Therefore, we reconstruct the BD using

while changing the estimated parameter vector $\gamma $ in the estimated parameter space. Here, $\delta \beta ^$ is calculated by Eq. (8) with $\gamma $. We obtain the reconstructed BD by gradually changing $\gamma $ from $\gamma n$ to $\gamma n+1$; this range corresponds to a parameter range $[pn,pn+1]$ of the original BD. It must be noted that this method can be extended to extrapolation of the reconstruction of BDs.^{1}

### C. Estimation of Lyapunov exponents for bifurcation-diagram reconstruction

Herein, we use two methods to estimate Lyapunov exponents. We begin by explaining how to estimate Lyapunov exponents from time-series data. Although this method is relatively imprecise, it can be used to estimate Lyapunov exponents for BDs generated by electronic circuits. Herein, we use the method proposed in 2012 by Yao *et al.*^{13} for estimating the largest Lyapunov exponent. In particular, this estimation method can be used with time-series data that are influenced by dynamical and observational noise. The other method involves estimating Lyapunov exponents from the reconstructed BD^{1} using the Jacobian matrix of the time-series predictor.^{18–20} Consequently, the second method estimates Lyapunov exponents more precisely. However, because the present target systems are one-dimensional maps, instead of the Jacobian matrix we use the derivative of nonlinear function of the predictor.

We begin by explaining the method for estimating Lyapunov exponents from time-series data. This method measures directly the extension rates of the distances between points on pairs of $2Np$ orbits, where $2Np$ is the number of points falling into an $\epsilon $ neighborhood of a point on the target orbit. Herein, we take $\epsilon =0.025$ and use the orbit of embedded time-series data with an embedding dimension of one. The algorithm of this estimation method is as follows:

Select an initial point $p(target)(t)\u2208R$.

Take all points $pi(\epsilon )(t)$, $(i=1,2,\u2026,2Np)$ falling into the $\epsilon $ neighborhood of $p(target)(t)$.

Form the pairs of $pi(\epsilon )(t)$, $(i=1,2,\u2026,Np)$ and $pi+Np(\epsilon )(t)$.

- Calculate the mean vector $d(\epsilon )(t)$ of the distance between the pairs bywhere $|\u22c5|$ is the absolute value.(11)$d(\epsilon )(t)=1Np\u2211i=1Np|pi(\epsilon )(t)\u2212pi+Np(\epsilon )(t)|,$
Calculate the mean vector $d(\epsilon )(t+1)$ of the distance between the pairs from $pi(\epsilon )(t+1),(i=1,2,\u2026,2Np)$.

Repeat steps 2–5 until $t=\psi +1$.

- Estimate the Lyapunov exponent by(12)$\mu =1\psi \u2211t=1\psi ln\u2061d(\epsilon )(t+1)d(\epsilon )(t).$

Herein, we use time-series data with no transients; therefore, this method cannot estimate Lyapunov exponents from periodic time series because $d(\epsilon )(t)=0$. Consequently, we set the Lyapunov exponents of periodic time series estimated by this method to zero.

Next, we explain the method for estimating Lyapunov exponents by using the derivative of nonlinear function $g(\u22c5,\u22c5)$ in Eq. (10). In this method, the Lyapunov exponent is obtained by

## III. ELECTRONIC-CIRCUIT REALIZATION IN DISCRETE-TIME DYNAMICAL SYSTEMS

We generate time-series data for discrete-time dynamical systems by using a CMOS (complementary metal-oxide-semiconductor) integrated circuit with pulse-width modulation (PWM).^{11,21} The circuit achieves arbitrary discrete-time nonlinear dynamics by voltage waveform sampling for nonlinear transformation.^{22,23} Figure 2 shows the principle of the voltage waveform sampling, where $\tau $, $Vx(t)$, $Vnon(\tau )$, and $Vrmp(\tau )$ are the continuous-time, the state-variable voltage, the nonlinear voltage waveform, and the ramped reference voltage for voltage to PWM conversion, respectively. Here, $Vx(t)(\u221dx(t))$ is transformed into a PWM signal with a pulse width $Tx(t)(\u221dVx(t))$ by comparing $Vx(t)$ with $Vrmp(\tau )$, as shown in Fig. 2(b). The nonlinear voltage waveform $Vnon(\tau )$ is sampled with this PWM signal, which is given by

where $Vc$ is the voltage that is sampled to capacitor $C$. By considering $Vc$ as $Vx(t+1)$ at the next time step $t+1$, arbitrary discrete-time nonlinear dynamics can be achieved.

In this study, the target systems are the logistic map and the sine map. To realize the dynamics of these maps, we generate the following nonlinear voltage waveforms $Vnon(l)(\tau )$ and $Vnon(s)(\tau )$, respectively, by using an arbitrary waveform generator:

Here, $p(l)$ is a bifurcation parameter of the logistic map; $p1(s)$ and $p2(s)$ are bifurcation parameters of the sine map; $a1$, $a2$, and $b1$ are coefficients for converting the numerical model into nonlinear voltage waveforms; $c1$ and $c2$ are coefficients for adjusting the maximum value of the sine map; and $Vnonbt0$ is a bias voltage for the nonlinear voltage waveforms. We adjust the maximum value of the sine map to be $0.8$ in Eq. (16) because $Vnon(s)(t)$ must be less than unity in this electronic circuit; i.e., $c1=0.8$ and $c2=1.25$. We set $a1=2.88\xd7106$, $a2=5\xd7105$, $b1=0.38\xd7106$, $Vnonbt0=1.2$, and $Vrmp(\tau )=0.414\xd7106\tau +0.6[V]$.

## IV. EXPERIMENTAL RESULTS: LOGISTIC MAP

In this section, we present the results of reconstructing BDs from time-series data generated by the electronic-circuit realization of the logistic map. We begin by describing the experimental conditions, then we show the BDs generated by MATLAB$\xae$, the electronic circuit, and the reconstructed BD. We then compare the return plots between the time-series datasets generated by MATLAB and those generated from the electronic circuit and the time-series predictor. Finally, we show the Lyapunov exponents for these BDs.

### A. Experimental conditions

The logistic map is given by

where $p(l)$ is the bifurcation parameter value of the logistic map. We use $p(l)(n)$ as the value of $p(l)$ to generate the time-series dataset $Sn$. Here, we set $p(l)(n)$ as

The length of each time-series dataset is 1000, and we use an ELM that has four neurons in its hidden layer.

### B. Results of bifurcation-diagram reconstruction

First, we estimate the dimension of the parameter space from the contribution ratio. From the PCA results, the estimated dimension is unity because the contribution ratio of the first principal component is around $100%$. Next, we show the bifurcation path and locus in Figs. 3(a) and 3(b), respectively. These figures correspond to each other because the parameter-value relations of the bifurcation path and locus are the same. Therefore, we consider that the estimated bifurcation parameter space corresponds to the parameter space of the target system.

Figures 4(a) and 4(b) show the BDs without and with dynamical noise generated by MATLAB, respectively. Figure 4(c) shows the BD generated by the electronic circuit, and Fig. 4(d) shows the reconstructed BD. The difference between Figs. 4(a) and 4(c) is that the latter is influenced by dynamical and observational noise and its fine bifurcation structure cannot be seen; in particular, the window around $p(l)=3.85$ disappears because of dynamical noise, and the fluctuations of time series are observed because of observational noise. Furthermore, the maximum state values of each parameter value in Fig. 4(c) are smaller than those in Fig. 4(a); for example, the maximum state values with $p(l)=4.0$ in Figs. 4(a) and 4(c) are around 1.0 and 0.9, respectively. Here, Fig. 4(b) shows the BD generated by Eq. (17) with added dynamical noise; that is, $x(l)(t+1)=p(l)x(l)(t)(1\u2212x(l)(t))+\xi $, where $\xi $ is Gaussian noise whose mean and standard deviation are zero and $10\u22123$, respectively. Because comparing Figs. 4(c) and 4(b) shows that these BDs have similarity, we see that Fig. 4(c) is the BD of the logistic map perturbed by dynamical noise in the electronic circuit.

Figure 4(d) shows the results of reconstructing the BD from time-series datasets generated by the electronic circuit. The bifurcation structure in Fig. 4(d) is quite clear, including the windows. We see that Fig. 4(d) is more similar to Fig. 4(a) than it is to Fig. 4(c) even though time-series datasets used to generate Fig. 4(c) are also used to produce Fig. 4(d). However, the maximum state values for each parameter value in Fig. 4(d) are approximately the same as the ones in Fig. 4(c). We therefore obtain the filtered BD by using time-series data influenced by dynamical and observational noise to reconstruct the BD.

### C. Comparison of return plots

We consider why the filtered BD is obtained only from time-series data influenced by dynamical and observational noise. In this section, we compare the return plots between the time-series datasets generated by MATLAB, the electronic circuit, and the time-series predictor, and we show the return plots when the path and locus numbers are one (i.e., $p(l)=3.55$ and $\gamma 1\u2243\u2212779$) and five (i.e., $p(l)=3.85$ and $\gamma 1\u2243909$) to confirm the apparent difference in the return plots of these numbers.

Figures 5(a) and 5(b) show the return plots when the path and locus numbers are one and five, respectively. Here, “$+$” and “$\u22c5$” are points of time series generated by MATLAB and the electronic circuit, respectively, and “$\xd7$” are points of time-series predictors. First, we compare “$+$” with “$\u22c5$” for both cases. Although these dynamics and parameter values are the same, the attractors and the dynamical ranges are different. For example, in Fig. 5(a), “$+$” is period-$8$ but “$\u22c5$” is chaotic. Note that as estimated in Sec. IV D, the Lyapunov exponent of “$\u22c5$” with this parameter value is positive. Next, we compare “$\xd7$” with “$\u22c5$” as well. In Fig. 5(a), we see that the output of the time-series predictor is period-4 even though the predictor was trained to model the “$\u22c5$” time series, which is chaotic. For this parameter value, we consider that the periodic time series is correct because the dynamical pattern of “$\u22c5$” is changed by dynamical noise.^{24} The time-series predictor has a filtering effect and models period-$4$ according to the good generalization of ELM. Comparing “$\xd7$” with “$\u22c5$” in Fig. 5(b), we see that “$\xd7$” and “$\u22c5$” almost overlap. Figure 6 shows the enlargement of a portion in Fig. 5(b). From the enlarged figure, we see that “$\u22c5$” is influenced by dynamical and observational noise, whereas “$\xd7$” might be filtered by the good generalization of the ELM. Therefore, we see that using a time-series predictor that has good generalization produces a filtered BD.

### D. Results of estimating Lyapunov exponents

We compare the Lyapunov exponents estimated using the following two methods:

the estimation method with time series for the BD generated by the electronic circuit and the reconstructed BD and

the estimation method with the derivative for the BD generated by MATLAB and the reconstructed BD.

Figures 7(a) and 7(b) show the Lyapunov exponents for the BD generated by the electronic circuit and the reconstructed BD from time-series data, respectively. Comparing Figs. 7(a) and 7(b), we see that the Lyapunov exponents for the reconstructed BD are close to those for the BD generated by the electronic circuit, except for the windows of the BD. Therefore, reconstructing the BD from time series influenced by dynamical and observational noise is successful to some extent.

Figures 8(a) and 8(b) show the Lyapunov exponents for the BD generated by MATLAB and the reconstructed BD using the derivative, respectively. We obtain more-precise estimates of the Lyapunov exponents in Fig. 8 than we do in Fig. 7. Comparing Figs. 8(a) and 8(b), the numbers of windows and the maximum values of the Lyapunov exponents agree closely. From this result, we see that the BD reconstructed from time series influenced by dynamical and observational noise is also quantitatively close to the BD generated by MATLAB.

## V. EXPERIMENTAL RESULTS: SINE MAP

In this section, as with the results for the logistic map, we present the BDs reconstructed from time-series data generated by the electronic-circuit realization of the sine map, along with the associated Lyapunov exponents.

We begin by describing the experimental conditions, then we show the BD generated by MATLAB, the one produced by the electronic circuit, and the reconstructed BD. Finally, we show the Lyapunov exponents for these BDs.

### A. Experimental conditions

The sine map in this study is given by

which corresponds to the arbitrary waveform generator in Eq. (16), where $p1(s)$ and $p2(s)$ are the bifurcation parameter values of the sine map and $p2(s)$ is fixed as zero in these numerical experiments. We use $p1(s)(n)$ as the parameter value $p1(s)$ to generate the time-series dataset $Sn$. Here, we set $p1(s)(n)$ as

The length of each time-series dataset is 1000, and we use an ELM that has four neurons in its hidden layer.

### B. Results of bifurcation-diagram reconstruction

First, we estimate the dimension of the parameter space from the contribution ratio. From the PCA results, we estimate that the dimension is one because the contribution ratio of the first principal component is around $100%$. Next, we show the bifurcation path and locus in Figs. 9(a) and 9(b), respectively. We see that these figures correspond to each other because the parameter-value relations of the bifurcation path and locus are the same. Therefore, the estimated bifurcation parameter space corresponds to the parameter space of the target system.

Figures 10(a) and 10(b) show the BDs generated by MATLAB and the electronic circuit, respectively, and Fig. 10(c) shows the reconstructed BD. Comparing Figs. 10(a) and 10(b), the latter is clearly influenced by dynamical and observational noise that is obscuring its bifurcation structure; in particular, the large window around $p(l)=2.96$ disappears. In addition, the state values of each parameter value in Fig. 10(b) are generally smaller than those in Fig. 10(a).

Figure 10(c) shows the results of reconstructing the BD from time-series data sets generated by the electronic circuit. The bifurcation structure in Fig. 10(c) is quite clear, including the windows. We see that Fig. 10(c) is more similar to Fig. 10(a) than it is to Fig. 10(b) even though time-series datasets were used to generate Fig. 10(b) are also used to produce Fig. 10(c). However, the maximum state values for each parameter value in Fig. 10(c) are approximately the same as the ones in Fig. 10(b). We, therefore, again obtain the filtered BD by using time-series data influenced by dynamical and observational noise to reconstruct the BD.

### C. Results of estimating Lyapunov exponents

As in Sec. IV D, we compare the Lyapunov exponents obtained using two estimation methods.

Figures 11(a) and 11(b) show the Lyapunov exponents for the BD generated by the electronic circuit and the reconstructed BD from time-series data, respectively. Comparing Figs. 11(a) and 11(b), we see that the Lyapunov exponents for the reconstructed BD are close to the ones for the BD generated by the electronic circuit, except for the windows of the BD. Therefore, reconstructing the BD from time series influenced by dynamical and observational noise is again successful.

Figures 12(a) and 12(b) show the Lyapunov exponents for the BD generated by MATLAB and the reconstructed BD using the derivative, respectively. Here, we obtain more-precise estimates of the Lyapunov exponents in Fig. 12 than in Fig. 11. Comparing Figs. 12(a) and 12(b), the numbers of windows and the maximum values of the Lyapunov exponents agree well. From this result, we see that the BD reconstructed from time series influenced by dynamical and observational noise is also quantitatively close to the BD generated by MATLAB.

### D. Results of reconstructing a two-dimensional bifurcation diagram

In this section, we reconstruct the two-dimensional BD of the sine map. We use $p1(s)(n)$ and $p2(s)(n)$ as the parameter values $p1(s)$ and $p2(s)$, respectively, to generate the time-series dataset $Sn$. Here, we set $p1(s)(n)$ and $p2(s)(n)$ as given in Table I. The other experimental conditions are the same as those in the numerical experiments for the one-dimensional BD of the sine map.

. | n
. | |||||
---|---|---|---|---|---|---|

. | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . |

$p1(s)(n)$ | 2.700 | 2.700 | 2.700 | 2.800 | 2.800 | 2.800 |

$p2(s)(n)$ | 0.000 | 0.050 | 0.095 | 0.000 | 0.050 | 0.095 |

. | n
. | |||||
---|---|---|---|---|---|---|

. | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . |

$p1(s)(n)$ | 2.700 | 2.700 | 2.700 | 2.800 | 2.800 | 2.800 |

$p2(s)(n)$ | 0.000 | 0.050 | 0.095 | 0.000 | 0.050 | 0.095 |

We show the results of reconstructing the two-dimensional BD of the sine map. First, we estimate the dimension of the parameter space from the contribution ratio. From the PCA results, we estimate that the dimension is two because the cumulative contribution ratio of the first and second principal components is around $100%$. Figures 13(a) and 13(b) show the BD generated by MATLAB and the reconstructed BD, respectively. In the figures, periodicities in each time-series data are indicated by a color bar. Here, the yellow points represent period-20 or greater, including chaos. We see that the estimated parameter space is partially expanded and contracted, but these figures have similarity nevertheless. In particular, we obtain the bifurcation structure in the global area of the two-dimensional parameter space, even though we used only time-series datasets generated with Table I.

## VI. CONCLUSION

Whereas the target systems of previous studies were limited to continuous-time dynamical systems, we have now reconstructed BDs from time-series data generated by electronic circuits in discrete-time dynamical systems. We obtained a clear bifurcation structure from the reconstructed BD. Comparing the reconstructed BD with the BDs generated by MATLAB and the electronic circuit, the reconstructed BD is more similar to the BD generated by MATLAB than it is to the BD generated by the electronic circuit. From these results, we obtained the filtered BD by reconstructing the BD using time-series data influenced by dynamical and observational noise. In addition, we also reconstructed a two-dimensional BD that was not reconstructed in previous studies, only from time-series data generated by the electronic circuit.

Moreover, we estimated the Lyapunov exponents of the BDs using two methods. We compared the Lyapunov exponents of the reconstructed BD with those of the BD generated by the electronic circuit and using the estimation method from time-series data. We also compared the Lyapunov exponents of the reconstructed BD with those of the BD generated by MATLAB using the estimation method with the derivative. By comparing the Lyapunov exponents, we think that reconstructing the BD from time series influenced by dynamical and observational noise is successful because the Lyapunov exponents of the reconstructed BD are approximately the same as those of the other BDs well.

From these results, we confirmed that our BD reconstruction is robust against dynamical and observational noise. Here, the results may be obtained with the gradient decent learning such as the method proposed by Tokunaga *et al.*^{3} However, we think that the ELM is suitable for the BD reconstruction since the ELM has better generalization property, faster learning speed, and fewer parameters than networks trained with the gradient decent learning. This is because the BD reconstruction is not always successful even when the target system is simple, and the parameter tuning for the time-series predictor is required.^{12}

Consequently, the method used in this study is useful to analyze measured time-series data influenced by dynamical and observational noise obtained from unknown dynamical systems. Here, it is significant that we use the noisy time-series data because time-series data measured from physical phenomena including operation of electronic circuits are always influenced by dynamical and observational noise.

In the future work, we will reconstruct the BDs from time-series generated by electronic circuits of non-autonomous systems that exhibit deterministic chaos.

## ACKNOWLEDGMENTS

This research is partially supported by JSPS KAKENHI (No. 15H05707), AMED (No. JP19dm0307009), WPI, MEXT JAPAN, and NEC Corporation.

## REFERENCES

*2018 International Symposium on Nonlinear Theory and Its Applications (NOLTA)*(IEICE, 2018), p. 439.

*International Joint Conference on Neural Networks (IJCNN)*(IEEE, 2017), p. 1809.

*Proceedings of ELM-2017*, edited by J. Cao, C. Vong, Y. Miche, and A. Lendasse (Springer, 2019), p. 58.

*European Conference on Circuit Theory and Design (ECCTD)*(IEEE, 2011), p. 126.