Piecewise-Linear Identification of Nonlinear ... - Semantic Scholar

1542

IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 54, NO. 7, JULY 2007

Piecewise-Linear Identification of Nonlinear Dynamical Systems in View of Their Circuit Implementations Oscar De Feo and Marco Storace, Member, IEEE

Abstract—We address here an aspect of the problem concerning circuit implementations of nonlinear dynamical systems that depend on control parameters. In particular, the problem of the identification of such systems is addressed in two steps. The first step uses the state-space reconstruction (through time delay reconstruction associated with principal component analysis) on the basis of scalar time series measured in the systems to be identified. The second step deals with the approximation of the flow in the reconstructed space (through a piecewise-linear approximation technique). The proposed method is first validated with two examples concerning known systems and then applied successfully to two realistic cases. Index Terms—Bifurcation, function approximation, identification, modeling, nonlinear systems, piecewise-linear (PWL) approximation techniques.

I. INTRODUCTION

T

HIS paper is concerned with modeling and circuit synthesis for a class of coupled nonlinear dynamical systems, which have attracted attention in the last few years. Such systems are described as networks consisting of interacting elementary units [1], which in turn are often modeled by low-dimensional ordinary differential equations (ODEs). Examples can be found in many fields: from thermodynamics [2], to artificial intelligence (neural networks) [3], neuronal and gene networks [4], computer simulation (Internet models) [5], psychology (cognitive systems) [6], sociology and anthropology (social organizations, group dynamics) [7], economics (heterogeneous agents models) [8], ecology and virology (population dynamics) [9], biology (cellular and metabolic networks) [10], electrical power grids (blackouts modeling) [11], technology [12], and traffic [13]. The degree of complexity of such models increases with the size of the network, i.e., the number of elementary units and, consequently, the number of ODEs. As a result, computer simuManuscript received August 2, 2004; revised September 7, 2005, April 9, 2006, and December 21, 2006. This work was supported by the MIUR, within the PRIN framework under Project 2004092944 and Project 2006093814, by the University of Genoa, and by the European Project APEREST: IST-200134893 and OFES-01.0456. This paper was recommended by Associate Editor L. Trajkovic. O. De Feo is with the Department of Microelectronic Engineering, University College Cork and Tyndall National Institute, Cork, Ireland (e-mail: Oscar. [email protected]). M. Storace is with the Biophysical and Electronic Engineering Department, University of Genoa, I-16145 Genoa, Italy (e-mail: [email protected]). Digital Object Identifier 10.1109/TCSI.2007.899613

lations for analyzing these models or investigating possible new behaviors require increasingly more: • computing resources, which is the limiting factor when we deal with large numbers of elementary units [2]; • computation time, which is the limiting factor when we consider dynamical systems used to process fast time-varying inputs. Such limitations may be overcome by using circuits governed by the original ODEs (or by good approximations of them) and having the same intrinsically parallel structure as the complex dynamical systems to be simulated. In other words, such circuits consist of a very large number of interconnected elementary units that are able to process large amounts of data in parallel, thus overcoming the bottleneck imposed by serial digital processing. The resulting circuit networks may have either regular or irregular structure. In the first case, each state component interacts directly only with a prescribed limited number of neighbor state components (defining a regular template) in the network of elementary dynamical systems. Chua’s cellular nonlinear networks [14] belong to this class. In the second case, the interactions are still limited, but they do not define a locally regular structure, as in the case of the small worlds and food webs [1]. This general (circuit) modeling problem requires an algorithmic synthesis procedure. Regarding this requirement, the piecewise-linear (PWL) approximation technique proposed in [15]–[18] provides a synthesis technique to find direct circuit implementations (either analog [19] or mixed-signal [20]) of PWL models. By assuming the existence of a reliable circuit synthesis technique for the PWL models considered, the problem reduces to the approximation/identification of such models. Such an approximation/identification problem can be tackled with various levels of knowledge about the dynamical systems to be emulated, thus resulting in either a direct “white-box” (approximation) or an inverse “black-box” (identification) modeling problem. In the direct problem [21], [22], a model of the elementary unit is known and a PWL approximation is derived from it. Examples of direct problems include circuit approximations of dynamical systems used to process sequences of images [18], [23] or bio-inspired neural networks based on biological neuron models, such as the Hodgkin–Huxley model of the axon membrane [17], [24], [25]. By resorting to continuation methods [26] applied to smoothed versions of the PWL approximations, we showed in [21] that the PWL technique can provide approximations that are both qualitatively and quantitatively similar to the

1549-8328/$25.00 © 2007 IEEE

DE FEO AND STORACE: PWL IDENTIFICATION OF NONLINEAR DYNAMICAL SYSTEMS

1543

Fig. 1. Block scheme of the identification/approximation method. (a) State-space reconstruction. (b) PWL approximation.

original systems. In [27], also a more elaborate multi-resolution algorithm has been proposed. In the case of the inverse problem, tackled in [28] with applications to three known systems, we deal with the identification of black-box dynamical systems from measured time-series. For instance, this is one of the main problems to be solved in the synthesis of networks for pattern classification based on chaotic dynamics proposed in [29]–[32]. Indeed, in such networks one needs to identify properly a set of ODEs that mimic the behavior of a given black-box system by using a class of models that is suitable for circuit implementation. Although the problem of identifying black-box systems from time series was introduced in [33], the main results in the field of nonlinear dynamical systems have been obtained quite recently, starting from the early 1980s, when works on state-space reconstruction were developed by Packard et al. [34], Ruelle [35], and Takens [36]. Time series analysis and state-space reconstruction were then developed [37]–[41] and applied to many fields to identify nonlinear dynamical systems from measured time series that are derived from the state of the system. In this paper, we use well-known techniques for the identification of nonlinear dynamical systems and combine these techniques with the PWL approximation method proposed in [15]–[18]. As a first step, we tackle the inverse problem of statespace reconstruction from observed data through “time delay reconstruction” exploiting principal component analysis (PCA); this technique is well known for the identification of linear and nonlinear systems [37], [42]. Once the state space has been reconstructed, we need a modeling technique to extract from the reconstructed data a model that approximates the original black-box system. We use the PWL approximation technique [15]–[18] to fit PWL models to the results of the state-space reconstruction. Using this combination of well-known techniques gives us two advantages. • The definition of a step by step method, which can be interpreted as the essential pre-requisite towards the design of

circuits that are able to emulate the behaviors of complex systems consisting of a very large number of elementary units. • The possibility of also identifying the dependence of the dynamical systems on control parameters. Concerning this last point, to validate the proposed method, we first test it by identifying known systems with multiple (even coexisting) invariants. The identification of such systems is performed by exploiting the available control parameters to explore all the parameter space regions characterized by qualitatively different dynamical behaviors. By properly modulating the parameters, we can build a time series containing information about all system invariants. Finally, the proposed procedure is applied to two cases where models of the dynamical systems to be identified from measured time series are not available. The rest of the paper is organized as follows. We introduce our technique for PWL identification of nonlinear dynamical systems in Section II. In Section III, we apply our method to four examples. Finally, in Section IV, the results are discussed and conclusions are drawn. II. TWO-STEPS IDENTIFICATION METHOD In this section, we introduce the basic elements of the identification/approximation problem considered. Our method consists of two steps, corresponding to the two sequences of operations shown in Fig. 1(a) and (b), i.e., state-space reconstruction and PWL approximation applied to the reconstructed state space, respectively. The goal of the method is to find a black-box model of a dynamical system of the form (1) where the state vector , the parameter vector and are column vectors the (continuous) flow ( and are the dimensions of the state space and parameter

1544

IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 54, NO. 7, JULY 2007

space, respectively), is a limited compact domain, is a scalar is a scalar function. (output) variable, and which have We assume that we have samples been obtained by sampling, with a constant step , a signal that contains additive measurement noise . We wish to find a candidate black-box system of the form . As is unknown, in (1) which could have generated general, we cannot reconstruct directly so we aim, instead, to reconstruct a (output predictive) topologically equivalent model of the state space, of the form

(2) where the reconstructed state vector and the (continare column vectors ( is the uous) flow dimension of the reconstructed state space), is the image of the compact domain in the reconstructed state space, and is the reconstructed output variable, obtained as a linear combination of the components. The model (2) should preserve the most significant invariants of the system (1), such as its state-space dimension, geometrical invariants, and Lyapunov exponents. A. Step 1: State-Space Reconstruction is Gaussian noise, the state-space reBy assuming that construction can be performed effectively using a combination of PCA [39], [41] and time delay reconstruction [36], [37], [43], [44]. Without loss of generality, and to keep the identification independent of the average and absolute energies, the time series is unbiased and normalized, thus obtaining a zero-mean . Through the time delay reand unit-variance time series -sample trajectory, from the construction we construct an time series , in an over-dimensioned ( -dimensional) state space. Such an -point trajectory is time delay reconstructed . The choice from the original time series with a lag time (embedding dimension) can be critical. Nevof both and ertheless, such parameters can be estimated through statistical techniques such as the mutual information [37], [44]. The main idea of the time delay reconstruction process [36], [37], [43], [44] is that we do not need the time derivatives in order to form a coordinate system to capture the structure of the orbits in the state space. We could, instead, use the lagged time series directly. In a nonlinear system, each element of such a collection of time lags is some (unknown) nonlinear combination of the actual physical variables that comprise the source of the measurements. Although we do not know the relationships between physical variables and measurements, it is irrelevant for the purpose of creating a topologically equivalent state space [36], [43]. In pattern recognition, PCA is widely used to extract features from a cluster of data. The attraction of PCA is its ability to project a collection of random vectors onto a linear subspace while preserving their variances. In the state-space reconstruction framework, PCA can be used to down-project an oversized state space along its principal components. The eigen decomposition of stationary discrete-time signals is also known as Hotelling transform [45]. For continuous-time stationary

random processes, the eigen decomposition method was formulated by Karhunen [46] and subsequently generalized by Loeve [47]. A summary of the Karhunen–Loeve expansion theory can be found in [48], [49] and references therein. The combination of time-delay reconstruction and PCA onto the projects a large delay embedding [34], [35] of principal eigendirections of the autocovariance matrix of the time series, simultaneously providing a denoising effect. This noise reduction is important for the subsequent flow estimation step, which is based on numerical differentiation of the reconstructed trajectory (see Section II-B-3). After performing this first step of the identification/approximation process, which is summarized in the block scheme shown in Fig. 1(a), we find the reconstructed state vector and a coefficient vector that minimizes the distance between and . By applying PCA to time delay reconstruction, in this first (colstep, we aim to extract principal components lected in the vector ). The number is chosen to be equal to the embedding dimension we would estimate through a method such as the Cao’s averaged false nearest neighbor estimator [50]. The projection of the -dimensional trajectory, provided by the time delay reconstruction, onto the -dimensional subspace defined by the PCA provides the reconstructed state trajectory. Of course, if the original trajectory is obtained for a fixed parameter . On the contrary, if the configuration, we expect to find original trajectory is obtained by varying control parameters, an increased dimension can be found, since state variables and parameters play the same role, from an identification point of view. Nevertheless, if the parameters are varied randomly and/or slowly with respect to the system dynamics, the risk of finding a number of principal components larger than is very low. B. Step 2: PWL Approximation Method In the second step [summarized in Fig. 1(b)], the samples of the reconstructed state-space vector , or a nonlinear transformation of them (see Section II-B-2), are used to estimate, by . numerical differentiation, samples of the flow This operation justifies the denoising requirement stated above, as the numerical differentiation acts as a high-pass filter on the samples of . This produces a PWL approximation for each component of the flow as a linear combination of a set of PWL basis functions

The -dimensional compact domain of and , i.e., , where defined as ) of the form hyperrectangle (rectangle, if

(3) is is a

(4) Each component of the domain can be subdivided into subintervals of amplitude , leading to a boundary configuration , and can be easily partitioned into simplices [15]–[18].

DE FEO AND STORACE: PWL IDENTIFICATION OF NONLINEAR DYNAMICAL SYSTEMS

The critical elements in the PWL approximation method are: (3); 2) 1) the choice of samples for determining the weights the selection of a domain partition which results in an approximation error that is lower than a prescribed or desired threshold; 3) the computation of the “optimal” (in some sense) vectors . The automatic definition of a suitable domain partition is an open problem, and will not be treated in this paper. In the examples considered, the domain is chosen and partitioned on the basis of the available data, with the aim of obtaining reasonably of accurate approximations by using the minimum number basis functions. are used for the The samples of each component of nonparametric regression process to estimate the weights (the projections of the reconstructed flow components on the PWL basis ) that minimize the distance between the PWL flow and , as described in the next subsections. In summary, the state-space is reconstructed (Step 1) using standard techniques [37]. The novelty of this paper is to apply the PWL technique to the approximation of the reconstructed flow (Step 2) and to show that there are no particular restrictions with respect to other more diffused (smooth) approximation methods. The main limit of the PWL technique lies in the required to produce relatively large number of coefficients an accurate model of a given function . This drawback is counterbalanced by the straightforward circuit implementations of such models. 1) Sample Generation: A significant aspect of the identification procedure concerns the spatial distribution of the fitting samples. Some simplices of the domain will be densely covered by the reconstructed trajectory samples, whereas other simplices will contain few or no samples at all. This also depends on the hyperrectangular shape of the domain in the PWL technique considered, since it is unlikely that the state part of is . uniformly covered by the reconstructed state trajectories In order to find a function that provides a good approximation to the flow according to a chosen domain partition, should be almost uniformly distributed over the samples of the state part of . Therefore, we preprocess the reconstructed state vector. Generally speaking, in order to have a good (almost uniform) distribution of samples over the simplicial domain, we need transients starting from uniformly distributed initial conditions. Consequently, if we aim to identify systems with invariants whose number, shape, and stability depend on control parameters (contained in the vector ), we need a strategy to vary such parameters during the evolution of the observed time series. Such a strategy should provide the maximum information on the system dynamics with the minimum number of samples of the reconstructed state-space trajectory. Moreover, the parameters should be varied on a time scale which is greater than that of the system dynamics, in order to obtain transient trajectories . in the observed signal The solution we adopted in this work is a “brute force” random variation of the control parameters over a significant region of the parameter space. The PWL model is obtained in the control space by using as domain points not only the

1545

Fig. 2. Bautin system. (a) Bifurcation diagram of the Bautin normal form (7); the curves T and H are a fold of cycles and Hopf bifurcation curves, respectively; the black points along the circular line indicate the pairs (p ; p ) reached during the driving of the control parameter ; the  coordinates of the three black points a, b, and c are:  = =8,  = =2, and  = 5=4 + =15. (b) Detail of the driving signal (t).

state trajectory samples, but also the corresponding parameter values. In so doing, the dependence of the model on the parameters is also identified. 2) Sample Distribution: The only case where the state vector eventually explores per se quite a large part of the state space is when it evolves on a chaotic attractor. However, in this case too the trajectory wanders through a limited part of the state space, visiting only a subset of the simplices, with a nonuniform distribution of samples. Therefore, we need to adapt the available samples to the sharp hyper-rectangular shape of the domain. In other words, we have to apply a nonlinear transformation to the in order to (i) adapt the smooth shape of the state trajectory chaotic attractor explored by the trajectory to the sharp shape of the domain and (ii) have a better (almost uniform) distribution of the samples over the domain used for the PWL approxi-dimensional domain in the mation. Hence, we map the control space visited by the vector to an hyperand distribute the available samples as cube

1546

IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 54, NO. 7, JULY 2007

Fig. 3. Bautin model identification. (a) Detail of the observed variable y~(t) = x (t) + n(t). (b) Invariant manifolds in the control space (; x ; x ) together with the trajectories obtained through the random driving of the control parameter . (c) Identified invariant manifolds in the control space (;  ;  ) and sample pointed out in Fig. 2. trajectories at the control parameter values 

uniformly as possible over it (according to their statistical properties). This should ensure a more reliable PWL approximation, as the regions of the domain without samples (without information) are compressed in a thin shell at the boundaries of the hypercube and hopefully we have samples in most of the simplices. Of course, since the parameters are directly controllable, it is sufficient to apply a nonlinear transformation only to the , the mapping for the parameters reconstructed state vector partiremaining linear. We remark that the lengths of the tions along the th state coordinate are equal in the hypercube and nonequal (according to the nonlinear inverse mapping) in the domain . With this in mind, in the second example proposed in Section III (concerning the Colpitts oscillator), we use the nonlinear invertible transformation

the hyperbolic tangent rather than a generic function (such as those provided by the inverse distribution method) lies in the simple circuit implementation of its inverse (widely used in conventional neural networks). This aspect is significant (from a circuit implementation point of view) whenever the observable of the identified system must be compared with the observable of the original system as, for instance, in synchronization, prediction, or control frameworks. As an alternative, we could make the comparison by transforming the observed signal. In such a case, the performance measures should be redefined according to the new metric. 3) Flow Estimation: From the reconstructed (and “sat, we urated,” if necessary) trajectory (domain) estimate samples of the flow vector (codomain) by numerical differentiation

(5)

(6)

where is estimated using the inverse cumulative distribution function method in order to distribute the samples almost uni. We note that, for peformly over the hypercube riodic-like strange attractors, the hyperbolic tangent is qualitatively similar to the functions provided by the inverse distribution (quantile) function method. The main advantage of using

implicitly depends on the control paramThe trajectory fitting samples are then used to find a flow eters . The that approximates . Of course, the availability of rich information is always desirable. Nevertheless, large amount of data may lead to memory problems during the computation of the weights of the PWL

DE FEO AND STORACE: PWL IDENTIFICATION OF NONLINEAR DYNAMICAL SYSTEMS

Fig. 4. Simulations of the identified Bautin model. Phase portraits at the control parameter values 

1547

chosen in Fig. 2.

model. In such a case, one can reduce the quantity of data by imposing a fitting factor of a fixed number of points per subdivision per dimension and by interpolating the available data to obtain fake samples over a more regular grid. The opposite problem (the presence of empty simplices) may be solved by linearly interpolating samples belonging to the nearest simplices. 4) Model Regression: Once a proper set of samples and a domain partition are obtained, we can derive the weights for the PWL approximation through optimization techniques. In the examples presented in Section III, we obtain the coefficients through the most direct solution, i.e., least-squares optimization [17]. III. EXAMPLES In this section, we apply the PWL identification method in order to verify that systems admitting different (even coexisting) attractors can be also identified [28]. The first two examples concern well-known systems (Bautin normal form and the Col, with pitts oscillator). We assume that , i.e., that the observable of the system to be identified is one of the components of the state vector. We assume, in addition, that the measurements include additive zero-mean Gaussian noise. In the examples we consider, the Gaussian noise has a standard deviation equal to 5% of the standard deviation of (corresponding to a signal-to-noise ratio of 13 dB). The last two examples are of a more applied nature. The third case concerns the identification of an autonomous chaotic dynamical system that is capable of modeling time series belonging to a given class (spoken vowels), for pattern recognition purposes [30]. In [30], this case was considered by using an identification technique based on a Lur’e type model. After successful identification, a generalized synchronization phenomenon is exploited to define a pattern recognition system [51], [52]. The fourth example concerns the identification from real patch clamp data [53] of a neuron model that depends on a control parameter. The advantage of using the PWL approximation technique in the second step of the proposed methodology lies in the corresponding direct circuit synthesis procedure, since, from a com-

Fig. 5. Colpitts system. (a) Detail of the random driving of the control parameter g (t). (b) Detail of the observed variable y~(t) = x (t) + n(t).

putational point of view, other methods are more effective and efficient. In all cases, the state-space reconstruction parameters are estimated using the numerical criteria given in Section II-A, whereas the numbers of subdivisions per dimension have been estimated heuristically, as a trade-off between model complexity and accuracy. Moreover, in the last two examples, are composed of several independent the time series observations. The measurements have been mean detrended and normalized in amplitude in the range [ 1,1].

1548

IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 54, NO. 7, JULY 2007

Fig. 6. Sample distribution for the identification of the Colpitts model. Projection in the reconstructed state space of: (a) the reconstructed trajectory  (t), (b) the saturated trajectory,  (t), and (c,d) Parallelepipeds in the state space containing points of the trajectories (a,b), respectively.

A. Bautin System The Bautin system is described by (7) The bifurcation diagram of system (7) (see [21], [26] for details) is shown in Fig. 2(a) for the sake of clarity. The trajectory of the system is forced by the parameters and , which step randomly amongst values that are uniformly distributed along the circle of unit radius shown in Fig. 2(a). In other words, the actual forcing parameter is the angle [Fig. 2(b)], which remains held at each value for a random time in the range [5, 15]. The projection of the trajectory in the parandomly visits the three regions , , rameter plane and [Fig. 2(a)], characterized by the presence of different invariants: a stable equilibrium point in , an unstable equilibrium point and a stable cycle in , and a stable equilibrium, a stable cycle, and an unstable cycle in . A detail of the time series (containing approximately points sampled with a time step ) of the observed variable is shown in Fig. 3(a). The corresponding trajectory in the control space is shown in Fig. 3(b) (it is easy to see the transients). The state-space reconstruction has been obtained by

, , and . The control parameter setting is added to the PCA space as a third dimension before the PWL approximation step. of the flow is obtained over The PWL approximation the domain , with , , , , , and . The domain is partitioned by and subdivisions along each performing dimension. Fig. 3(c) shows the result of the identification in the . By comparing Fig. 3(b) and (c), we see control space the relative effectiveness of the identification performed by our method. In particular, the state-space trajectories corresponding to the points , , and in Fig. 2(a) are clearly visible in Fig. 3(c) and reported in Fig. 4. B. Colpitts Oscillator The Colpitts oscillator, whose dynamics were analyzed in detail in [54], [55], is described by

(8)

DE FEO AND STORACE: PWL IDENTIFICATION OF NONLINEAR DYNAMICAL SYSTEMS

1549

Fig. 7. Colpitts model identification results. (a) Bifurcation diagram with respect to the parameter g for the original Colpitts system. (b) Bifurcation diagram with respect to the normalized parameter p for a full-knowledge PWL-approximated system (de-saturated). (c) Bifurcation diagram with respect to the normalized parameter p for the identified system (de-saturated). The p coordinates of the three dashed lines a, b, and c are: p = 0:2, p = 0:5, and p = 0:8. (d) Bifurcation diagram with respect to the normalized parameter p for the system obtained by PWL approximating the numerically estimated vector field (de-saturated).

In order to apply an identification procedure, the parameters have been set to the following values: , , and . The parameter is varied randomly in the range [0.37,0.47] to force the system to explore the control space. is shown in Fig. 5(a). The signal A detail of the signal remains held at each value for a random time in the range . A section of the time series (containing about samples) obtained by sampling with a fixed step is shown in Fig. 5(b). The state-space reconstruction has been performed by , , and . The normalized control choosing is added to the PCA space parameter as a fourth dimension before the PWL identification step. shown in The reconstructed state-space trajectory Fig. 6(a) is “saturated” according to the technique described in in transformation (5), yielding Section II-B-2, with shown in Fig. 6(b). From the latter, we the trajectory estimate the flow vector by numerical differentiation, thus fitting samples needed to determine . obtaining the Comparison between Fig. 6(c) and (d) indicates the improvement in the covering of the approximation space. All figures show projections in the reconstructed state space . In terms of the parallelepipeds deriving from the domain partition, Fig. 6(c) and (d) show that the coverage of the space is increased, as a result of the transformation, from about 60% to about 72%.

of the flow is obtained over The PWL approximation the normalized saturated domain , partitioned with subdivisions along each state coordinate and with subdivisions along the parameter axis. In the bifurcation diagrams shown in Fig. 7, the maxima of the original output [(a)] and the maxima of the of the three systems reconstructed output variables [(b)–(d)] are plotted. The three PWL approximations represented in Fig. 7(b)–(d) are obtained with different levels of knowledge about the original system. Fig. 7(a) and (c) show the bifurcation diagrams (with respect to and , respectively) of the original system and the system , respectively. In this identified from the noisy time series case we reconstruct the state space and we have measurement noise. The two diagrams exhibit a good degree of qualitative similarity, even if there are several quantitative differences. In particular, Fig. 8 shows a comparison between the asymptotic state trajectories corresponding to the three values of the con, , and trol parameter labelled as in the original system (first column), the saturated system (second column), the PWL-identified system (third column), and the de-saturated PWL system (fourth column). The identification process changes some details of the bifurcation diagram, preserving its qualitative shape and the qualitative features of the flow. The Feigenbaum route to chaos is not as evident in

1550

IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 54, NO. 7, JULY 2007

Fig. 8. Simulations of the identified Colpitts model: phase portraits at the control parameter values p pointed out in Fig. 7. First column: original trajectories projected in the reconstructed state space  (t). Second column: saturated trajectories  (t). Third column: trajectories of the PWL approximation of the saturated flow. Fourth column: de-saturated trajectories.

Fig. 7(c) as in Fig. 7(a), and chaos appears for a value of that is lower than the corresponding value of that marks the appearance of chaotic dynamics in Fig. 7(a). Also, the merging of the two main chaotic intervals into a single one occurs at a lower value of the bifurcation parameter . There are three potentially critical phases to which one can ascribe the quantitative differences in the dynamical behaviors of the original and the approximate system, i.e., the state-space reconstruction, estimation of the vector field from samples of a trajectory, and the PWL approximation. One might wonder if there is a phase which affects the final result more than the others. Fig. 7(b) and (d) shows the bifurcation diagrams of two further PWL approximations of the original system. The numerical parameters related to scaling of the data, the saturation procedure, and the simplicial partition are the same as for the proposed identification result in Fig. 7(c). However, in Figs. 7(b) and (d) the output variables concern dynamical are PWL approximations of systems whose vector fields the original noiseless vector field projected along the principal components of the original state space , with full knowledge of the derivatives (b) or with the derivatives estimated numerically

(d). In both cases we do not reconstruct the state space and we do not have measurement noise. In Fig. 7(b), the Feigenbaum route to chaos is quite evident. However, the two main chaotic intervals do not merge into a single one. As in Fig. 7(c), chaos appears for a lower value of than the corresponding value of which marks the appearance of chaotic dynamics in Fig. 7(a). The main differences are in the left part of the plots, where the approximation induces the presence of a second stable limit cycle coexisting with the one homologous to the limit cycle present in the original system and shown in Fig. 8 (first row and first column). The abrupt in Fig. 7(b) is evidence of transition occurring at about a change in the basin of attraction for the initial condition chosen to generate the bifurcation diagrams. It is not easy to decide which of the results in Fig. 7(b) and (c) is more representative of the real system. This is a sign of the strong influence of the PWL approximation on the final result. Finally, in Fig. 7(d) a further early appearance of chaotic dynamics and the coexistence of stable limit cycles for low values of the bifurcation parameter are evident. In conclusion, the main differences between the dynamical behaviors of the original and the approximate systems arise

DE FEO AND STORACE: PWL IDENTIFICATION OF NONLINEAR DYNAMICAL SYSTEMS

1551

during the PWL approximation phase. Of course, the larger the partition numbers, the more accurate the PWL approximations and the closer dynamical behaviors of the approximate systems. of partitions, As shown in [21], by increasing the numbers the approximation level shifts from qualitative to quantitative.

C. Vowel Signal The time series considered consists of measured vowels [a]’s spoken by a single speaker. The data have been obtained by recording 50 sustained (2 s each one) pronunciations of the vowel [a], with a normal microphone, at a sampling frequency of 44.1 kHz, and at a resolution of 16 bits/sample (see also [30]). Such signals are the same as those considered in [30]; similarly to there, they are processed here with the very same objective of identifying chaotic systems from pseudo-periodic signals, and subsequently exploiting chaotic synchronization to realize pattern recognizers. The signals have been pre-processed by applying a second-order Butterworth low-pass filter at 20 kHz. (containing A window of the dimensionless time series samples, with a sampling step approximately ms) is shown in Fig. 9(a). According to the estimation criteria recalled in Section II-A, after the time series is unbiased, , the state-space reconstruction is obtained by choosing , and . The PWL approximation of the flow is obtained over the domain , with and . The domain is partitioned by performing , and subdivisions along each state-space component. Finally, owing to the multiple initial points provided by the 50 available recordings, the reconstructed state trajectory explores the reconstructed state space quite well. No saturation is required for a successful identification in this case. Fig. 9(c) shows a projection in a three-dimensional subspace of a typical trajectory of the identified system. A fragment of the corresponding output time series is shown, for comparison, in Fig. 9(b). The maximum Lyapunov exponent of the identiand the box counting fractal fied system is (both are numeridimension of the chaotic attractor is cally estimated using the TISEAN tool [38]). Fig. 9(d) shows the peak-to-peak plot [56] corresponding to this time series and confirms that the identification leads to a chaotic dynamical system in spite of the unknown nature of the measured time series, as also shown in [30]. In other words, chaos self-emerges to model diversity. Finally, the identified model could be implemented in circuitry using the methods described in [19], [20], to obtain realtime pattern recognizers. However, this further step is beyond the scope of this paper.

D. Patch Clamp Signals The signals considered in this example are action potentials measured from a single neuron with a patch clamp recording

Fig. 9. Identification of a PWL model for a spoken vowel [a]. (a) Window of y g. (b) Window of the output time series f^ y g from the the time series f~ identified system. (c) Projection in a three-dimensional subspace of a typical trajectory of the identified system. (d) Peak-to-peak plot of the reconstructed y g. time series f^

technique. They are part of the data set considered in [53] (courtesy of Markram and Toledo1). The recording was performed as described in [57]–[59]. In brief, somatic whole-cell recordings ) were made. Signals were (with a pipette resistance of 13 sampled at 20 kHz, Butterworth second-order bandpass filtered and kHz), digitized using an ITC-18 in( terface (Instrutech, Great Neck, NY) and stored on the computer hard disk for off-line analysis (Igor Wavemetrics, Lake Oswego, OR). The data have been further pre-processed through one-dimensional wavelet de-noising. The measurements depend on the membrane input current (the control parameter). We identify a system depending on one parameter, as for the first two proposed examples. considered contains about samThe time series ples, with a sampling step ms. Three windows of the time series are shown in the second row of Fig. 10, where the output time series are labeled (instead of ) for consistency 1http://bmi.epfl.ch/LNMC.html

1552

IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 54, NO. 7, JULY 2007

Fig. 10. Identification of a PWL model for a collection of patch clamp signals measured from a single neuron. First row: input signals i(t). Second row: measured output signals. Third row: output signals fu ^ g (simply denoted as fu ^g) provided by the identified PWL system.

with the symbol adopted for the membrane potential in electrophysiology. The state-space reconstruction from the normalized , , and unbiased time series is obtained by choosing . The PWL approximation of the flow is and obtained over the domain , with and . The domain is partitioned by performing , and subdivisions along each domain component. Finally, as in the case of the Bautin system identification example (Section III-A), no saturation of the state space has been used. Fig. 10 shows the results of the identification and PWL approximation process. The first row of figures shows the corresponding to three measurements. The input signals second and third rows show the output signals measured from the neuron and provided by the obtained PWL system, respectively. The figure shows results corresponding to a suprathreshold input (column a), a subthreshold input (column c), and a time-varying ramp input (column b). The comparison between the second and third rows of the figures is evidence of a successful identification result. The measurements concerning subthreshold inputs were quite poor compared to the suprathreshold measurements. Nevertheless, the result reported in the bottom right figure shows good agreement with the real data. In this example, circuit synthesis of the identified PWL model could yield to a basic building block of an artificial brain, which currently can be simulated only by resorting to very complex computer simulations (see, for instance, the Blue Brain project [60]). IV. CONCLUDING REMARKS This paper can be viewed as a first step towards the definition of a method for PWL modeling/identification and circuit

synthesis of nonlinear black-box systems. In summary, we have combined the state-space reconstruction and the PWL approximation method to identify both ideal and realistic systems. In the first case, even though the considered systems were known, no information was used, except for the observed time series. The reported results show successful identification of all systems, both from state and parameter dependence points of view. In particular, the results obtained in the two realistic cases considered illustrate that the proposed methodology may be successfully applied to the identification of realistic systems (also depending on input control parameters, as in the case of the neuron model), whenever a circuit implementation is required. The flow approximation step could be taken by resorting to many different modeling techniques. In particular, a class of linear models such as radial basis functions (RBFs) [61], [62] can be used as an alternative to the PWL functions considered. RBFs are mainly used for classification problems and their circuit realizations [63], [64] can be quite cumbersome if large sets of RBFs are needed. From an approximation theory point of view, the (single-resolution) PWL basis functions considered could be less efficient than RBFs. Indeed, in the case of functions that are wrinkled only over limited regions of the domain, the regular (simplicial) partition of the domain, on which they rely, may result in a futilely over-dimensioned network. Nevertheless, the method proposed in this paper can also be applied directly to RBFs. A multi-resolution version of the PWL approximation technique has been proposed in [27]. Finally, the PWL approximation is obtained by exploiting only information about the state derivative over the domain . More refined results in terms of accuracy of the dynamical behaviors of the approximate system may be on the basis of obtained by deriving the “optimal” weights variational criteria that take into account a priori the expected dynamical behaviors [65]–[67].

DE FEO AND STORACE: PWL IDENTIFICATION OF NONLINEAR DYNAMICAL SYSTEMS

ACKNOWLEDGMENT The authors wish to thank H. Markram, M. Toledo and co-workers (see footnote 1) for providing the patch clamp signals data, and P. Kennedy and M. Parodi for useful comments and discussions. A part of this work was performed while the first author was with the Laboratory of Nonlinear Systems, Swiss Federal Institute of Technology Lausanne (EPFL), Switzerland.

REFERENCES [1] S. Strogatz, “Exploring complex networks,” Nature, vol. 410, pp. 268–276, Mar. 2001. [2] G. D’anna, P. Mayor, A. Barrat, V. Loreto, and F. Nori, “Observing brownian motion in vibration-fluidized granular matter,” Nature, vol. 424, pp. 909–912, Aug. 2003. [3] W. Maass, T. Natschläger, and H. Markram, “Real-time computing without stable states: A new framework for neural computation based on perturbations,” Neural Comput., vol. 14, no. 11, pp. 2531–2560, Nov. 2002. [4] T. Mestl, C. Lemay, and L. Glass, “Chaos in high dimensional neural and gene networks,” Phys. D, vol. 98, no. 1, pp. 33–52, Nov. 1996. [5] J. Eckmann and E. Moses, “Curvature of co-links uncovers hidden thematic layers in the world-wide web,” Proc. Nat. Acad. Sci., vol. 99, no. 9, pp. 5825–5829, Apr. 2002. [6] J. D. W. Tschacher, The Dynamical Systems Approach to Cognition—Concepts and Empirical Paradigms Based on Self-Organization, Embodiment, and Coordination Dynamics. Singapore: World Scientific, 2003. [7] R. Ferrer i Cancho and R. Solé, “The small world of human language,” Proc. Royal Soc. London B, vol. 268, no. 1482, pp. 2261–2265, Nov. 2001. [8] W. B. Arthur, “Complexity and the economy,” Science, vol. 284, no. 5411, pp. 107–109, Apr. 1999. [9] J. Montoya and R. Solé, “Small world patterns in food webs,” J. Theor. Biol., vol. 214, no. 3, pp. 405–412, Feb. 2002. [10] H. Jeong, B. Tombor, R. Albert, Z. Oltavi, and A. Barabási, “The largescale organization of metabolic networks,” Nature, vol. 407, no. 6804, pp. 651–654, Oct. 2000. [11] B. A. Carreras, V. E. Lynch, I. Dobson, and D. E. Newman, “Critical points and transitions in an electric power transmission model for cascading failure blackouts,” Chaos, vol. 12, no. 4, pp. 985–994, Dec. 2002. [12] R. F. i Cancho, C. Janssen, and R. Solé, “Topology of technology graphs: Small world patterns in electronic circuits,” Phys. Rev. E, vol. 64, no. 4, p. 046 119(1-5), Oct. 2001. [13] A. Visser and K. Nagel, “Traffic simulation,” Future Generation Comput. Syst., vol. 17, no. 5, pp. 625–626, Mar. 2001. [14] L. Chua, CNN: A Paradigm for Complexity, ser. A. Singapore: World Scientific , 1998, vol. 31. [15] P. Julián, A. Desages, and O. Agamennoni, “High-level canonical piecewise-linear representation using a simplicial partition,” IEEE Trans. Circuits Syst. I, Fundam. Theory Appl., vol. 46, no. 4, pp. 463–480, Apr. 1999. [16] P. Julián, A. Desages, and B. D’Amico, “Orthonormal high level canonical PWL functions with applications to model reduction,” IEEE Trans. Circuits Syst. I, Fundam. Theory Appl., vol. 47, no. 5, pp. 702–712, May 2000. [17] M. Storace, P. Julián, and M. Parodi, “Synthesis of nonlinear multiport resistors: A PWL approach,” IEEE Trans. Circuits Syst. I, Fundam. Theory Appl., vol. 49, no. 8, pp. 1138–1149, Aug. 2002. [18] M. Storace, L. Repetto, and M. Parodi, “A method for the approximate synthesis of cellular nonlinear networks—Part 1: Circuit definition,” Int. J. Circuit Theory Appl., vol. 31, no. 3, pp. 277–297, May–Jun. 2003. [19] M. Storace and M. Parodi, “Towards analog implementations of PWL two-dimensional nonlinear functions,” Int. J. Circuit Theory Appl., vol. 33, no. 2, pp. 147–160, Mar.–Apr. 2005. [20] M. Parodi, M. Storace, and P. Julián, “Synthesis of multiport resistors with piecewise-linear characteristics: A mixed-signal architecture,” Int. J. Circuit Theory Appl., vol. 33, no. 4, pp. 307–319, Jul.–Aug. 2005.

1553

[21] M. Storace and O. De Feo, “Piecewise-linear approximation of nonlinear dynamical systems,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 51, no. 4, pp. 830–842, Apr. 2004. [22] M. Storace and O. De Feo, “PWL approximation of nonlinear dynamical systems—Part I: Structural stability,” J. Phys., vol. 22, pp. 208–221, 2005. [23] L. Repetto, M. Storace, and M. Parodi, “A method for the approximate synthesis of cellular nonlinear networks—Part 2: Circuit reduction,” Int. J. Circuit Theory Appl., vol. 31, no. 3, pp. 299–313, May–June 2003. [24] A. Hodgkin and A. Huxley, “A quantitative description of membrane current and its applications to conduction in nerve,” J. Physiol., vol. 117, pp. 500–544, 1952. [25] M. Parodi and M. Storace, “On a circuit representation of the hodgkin and huxley nerve axon membrane equations,” Int. J. Circuit Theory Appl., vol. 25, no. 2, pp. 115–124, Mar.–Apr. 1997. [26] Y. Kuznetsov, Elements of Applied Bifurcation Theory, 2nd ed. New York: Springer-Verlag, 1998. [27] M. Gaggero, M. Parodi, and M. Storace, “Multiresolution PWL approximations,” presented at the Proc. Eur. Conf. Circuit Theory and Design (ECCTD), Cork, Ireland, Aug. 29–Sep. 2 2005, paper 051, unpublished. [28] O. De Feo and M. Storace, “PWL approximation of nonlinear dynamical systems Part—II: Identification issues,” J. Phys., vol. 22, pp. 30–42, 2005. [29] O. De Feo and M. Hasler, Modelling Diversity by Chaos and Classification by Synchronization. Dordrecht, Germany: Kluwer, 2003, pp. 25–40. [30] O. De Feo, “Self-emergence of chaos in the identification of irregular periodic behavior,” Chaos, vol. 13, no. 4, Dec. 2003. [31] O. De Feo, “Qualitative resonance of Shil’nikov-like strange attractors—Part I: Experimental evidence,” Int. J. Bifurc. Chaos, vol. 14, no. 3, pp. 873–891, Mar. 2004. [32] O. De Feo, “Qualitative resonance of Shil’nikov-like strange attractors—Part II: Mathematical analysis,” Int. J. Bifurc. Chaos, vol. 14, no. 3, pp. 893–912, Mar. 2004. [33] G. Yule, “On a method for investigating periodicites in disturbed series with special reference to Wolfe’s sunspot number,” Phil. Trans. Roy. Soc., vol. 226, pp. 267–298, 1927. [34] N. Packard, J. Crutchfield, J. Farmer, and R. Shaw, “Geometry from a time series,” Phys. Rev. Lett., vol. 45, no. 9, pp. 712–716, Sep. 1980. [35] D. Ruelle, Chaotic Evolution and Strange Attractors: The Statistical Analysis of Time Series for Deterministic Nonlinear Systems. Cambridge, MA: Cambridge Univ. Press, 1989. [36] F. Takens, Detecting Strange Attractors in Fluid Turbolence, ser. Lecture notes in mathematics. Berlin, Germany: Springer, 1981, pp. 366–381. [37] H. Kantz and T. Schreiber, Nonlinear Time Series Analysis, 2nd ed. Cambridge, MA: Cambridge University Press, 2003. [38] P. Hegger, H. Kantz, and T. Schreiber, “Practical implementation of nonlinear time series methods: The TISEAN package,” Chaos, vol. 9, no. 2, pp. 413–435, Jun. 1999. [39] H. Abarbanel, R. Brown, J. Sidorowich, and L. Tsimring, “The analysis of observed chaotic data in physical systems,” Rev. Modern Phys., vol. 65, no. 4, pp. 1331–1392, Oct. 1993. [40] J. Gibson, J. Farmer, M. Casdagli, and S. Eubank, “An analytic approach to practical state-space reconstruction,” Phys. D, vol. 57, no. 1–2, pp. 1–30, Jun. 1992. [41] M. Casdagli, S. Eubank, J. Farmer, and J. Gibson, “State-space reconstruction in the presence of noise,” Phys. D, vol. 51, no. 1–3, pp. 52–98, Aug. 1991. [42] P. van Overschee and B. De Moor, Subspace Identification of Linear Systems: Theory, Implementation, Applications. Dordrecht, The Netherlands: Kluwer Academic, 1996. [43] R. Mañé, On the Dimension of the Compact Invariant Set of Certain Non-Linear Maps, ser. Lecture notes in mathematics. Berlin, Germany: Springer, 1981, pp. 230–242. [44] M. Small, Applied Nonlinear Time Series Analysis: Applications in Physics, Physiology and Finance. Singapore: World Scientific, 2005. [45] H. Hotelling, “Analysis of a complex of statistical variables into principal components,” J. Educ. Psychol., vol. 24, pp. 417–441, 1933, 498–520. [46] K. Karhunen, “Uber lineare methoden in der wahrscheinlichkeitsrechnung,” Amer. Acad. Sci., Fennicade, Ser. A, I, vol. 37, pp. 3–79, 1947, translation: RAND Corporation, Santa Monica, California, Rep. T-131, Aug. 1960.

1554

IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 54, NO. 7, JULY 2007

[47] M. Loeve, Fonctions aléatories de second ordre. Paris: Herman, 1948. [48] A. Newman, Model Reduction via the Karhunen-Loeve Expansion Part I: An Exposition Institute for System Research, University of Maryland, Tech. Rep. 96-32, 1996. [49] A. Newman, Model Reduction via the Karhunen-Loeve Expansion Part I: Some Elementary Examples Institute for System Research, University of Maryland, Tech. Rep. 96-32, 1996. [50] L. Cao, “Practical method for determining the minimum embedding dimension of a scalar time series,” Phys. D, vol. 110, no. 1–2, pp. 43–50, Dec. 1997. [51] O. De Feo, “Transforming subharmonic chaos to homoclinic chaos suitable for pattern recognition,” Int. J. Bifurc. Chaos, vol. 15, no. 10, pp. 3345–3357, Oct. 2005. [52] O. De Feo, “Tuning chaos synchronization and anti-synchronization for applications in temporal pattern recognition,” Int. J. Bifurcation Chaos, vol. 15, no. 12, pp. 3905–3921, Dec. 2005. [53] M. Toledo-Rodriguez, B. Blumenfeld, C. Wu, J. Luo, B. Attali, P. Goodman, and H. Markram, “Correlation maps allow neuronal electrical properties to be predicted from single-cell gene expression profiles in rat neocortex,” Cerebral Cortex, vol. 14, no. 12, pp. 1310–1327, Dec. 2004. [54] G. Maggio, O. De Feo, and M. Kennedy, “Nonlinear analysis of the Colpitts oscillator and applications to design,” IEEE Trans. Circuits Syst. I, Fundam. Theory Appl., vol. 46, no. 9, pp. 1118–1130, Sep. 1999. [55] O. De Feo, G. Maggio, and M. Kennedy, “The Colpitts oscillator: Families of periodic solutions and their bifurcations,” Int. J. Bifurc. Chaos, vol. 10, no. 5, pp. 935–958, May 2000. [56] M. Candaten and S. Rinaldi, “Peak-to-peak dynamics: A critical survey,” Int. J. Bifurc. Chaos, vol. 10, no. 8, pp. 1805–1819, Aug. 2000. [57] H. Markram, J. Lubke, M. Frotscher, A. Roth, and B. Sakmann, “Physiology and anatomy of synaptic connections between thick tufted pyramidal neurones in the developing rat neocortex,” J. Physiol., vol. 500, no. 2, pp. 409–440, Apr. 1997. [58] A. Gupta, Y. Wang, and H. Markram, “Organizing principles for a diversity of GABAergic interneurons and synapses in the neocortex,” Science, vol. 287, no. 5451, pp. 273–278, Jan. 2000. [59] Y. Wang, A. Gupta, M. Toledo-Rodriguez, C. Wu, and H. Markram, “Anatomical, physiological, molecular and circuit properties of nest basket cells in the developing somatosensory cortex,” Cerebral Cortex, vol. 12, no. 4, pp. 395–410, Apr. 2002. [60] [Online]. Available: http://bluebrainproject.epfl.ch/ [61] M. Orr, Introduction to Radial Basis Function Networks intro.ps [Online]. Available: http://www.anc.ed.ac.uk/mjo/ [62] M. Orr, Recent Advances in Radial Basis Function Networks recad.ps [Online]. Available: http://www.anc.ed.ac.uk/mjo/ [63] L. Theogarajan and L. Akers, “A scalable low-voltage analog Gaussian radial basis circuit,” IEEE Trans. Circuits Syst. II, Analog Digit. Signal Process., vol. 44, no. 11, pp. 977–979, Nov. 1997. [64] J. Madrenas, M. Verleysen, P. Thissen, and J. Voz, “A CMOS analog circuit for Gaussian functions,” IEEE Trans. Circuits Syst. II, Analog Digit. Signal Process., vol. 43, no. 1, pp. 70–74, Jan. 1996. [65] L. Repetto, “A Piecewise-Linear Approach to Approximate Circuit Syntheses of Nonlinear Systems” Ph.D. dissertation, Biophysical and Electronic Engineering Department, University of Genoa, Genova, Italy, 2004 [Online]. Available: http://www.dibe-mail.dibe.unige.it/department/ncas/REPETTO/Thesis_C.pdf

[66] L. Repetto, M. Parodi, and M. Storace, “A procedure for the computation of accurate PWL approximations of nonlinear dynamical systems,” Int. J. Circuit Theory Appl., vol. 34, no. 2, pp. 237–248, Mar.–Apr. 2006. [67] M. Storace and F. Bizzarri, “Towards accurate PWL approximations of parameter-dependent nonlinear dynamical systems with equilibria and limit cycles,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 54, no. 3, pp. 620–631, Mar. 2007.

Oscar De Feo received the B.Sc. degree (summa cum laude) in industrial electronics from Maxwell Technical High School, Milan, Italy, the M.Sc. degree (summa cum laude) in computer science engineering from Politecnico di Milano, Milan, Ital, and the Ph.D. degree in Technical Sciences from the Swiss Federal Institute of Technology Lausanne (EPFL), Lausanne Switzerland, in 1990, 1995, and 2001, respectively. He is currently a Lecturer of the University College Cork , Cork, Ireland. In the past, he held the position of Research Associate at the EPFL and of Post-Doctoral Researcher at the Ecole Normale Supérieure Paris, Paris, France. He also held several Visiting Researcher positions throughout Europe: International Institute for Applied Systems Analysis (Laxenburg, Austria-IIASA); Fondazione Eni Enrico Mattei (Milan, Italy); Research Institute for Mathematics and Informatics (Amsterdam, The Netherlands). His research interests are at the interdisciplinary interface of applied nonlinear dynamical systems, machine learning, and life sciences. During his research career, he has combined linear and nonlinear modelling techniques, control theory, mixed analytical-numerical tools and experimental data analysis to address concrete problems spanning the fields of: population biology; evolutionary ecology; theoretical, experimental, and clinical neuroscience; machine learning; and industrial electronics. Dr. De Feo received the Mikhalevich Award from the IIASA, the Best Paper Award at the European Conference on Circuit Theory and Design (ECCTD’99), and the prize from the Chorafas foundation for the work accomplished during his Ph.D. thesis.

Marco Storace was born in Genoa, Italy, in 1969. He received the “Laurea” (M.Sc.) five-year degree (summa cum laude) in Electronic Engineering in March 1994 and the Ph.D. degree in Electrical Engineering in April 1998, both from the University of Genoa, Italy. Since 2004 he has been an Associate Professor in the Department of Biophysical and Electronic Engineering, University of Genoa. He was a visitor to EPFL, Lausanne, Switzerland, in 1998 and 2002. His main research interests are in the area of nonlinear circuit theory and applications, with emphasis on circuit models of nonlinear systems (e.g., hysteresis, biological neurons), methods for the piecewiselinear approximation of nonlinear systems and for the consequent circuit synthesis, methods of analysis and synthesis of cellular circuits, and nonlinear dynamics. He is the author or coauthor of about seventy scientific papers, more than an half of which have been published in international journals.