➠
➡ PWL IDENTIFICATION OF DYNAMICAL SYSTEMS: SOME EXAMPLES Oscar DE FEO ∗
Marco STORACE †
Laboratory of Nonlinear Systems Swiss Federal Institute of Technology Lausanne EPFL-IC-LANOS, CH-1015 Lausanne, Switzerland E-mail:
[email protected] Biophysical and Electronic Engineering Department, University of Genoa Via Opera Pia 11a, I-16145, Genova, Italy E-mail:
[email protected] ABSTRACT The problem of the identification of nonlinear dynamical systems with a piecewise-linear technique for the purpose of their circuit implementations is addressed. The identification is achieved on the basis of time series measured on the systems. Two academic examples are given, where the systems to be identified are completely known and the trajectories are properly obtained by controlling bifurcation parameters. 1. INTRODUCTION Nonlinear black-box modeling/identification is one of the two main steps of a very recently developed chaos-based temporal pattern recognition/classification network [1, 2]. In this paper, we shall take a step forward in this direction for the purpose of finding an effective circuit synthesis technique to implement such a kind of networks. Toward this end, we shall resort to the piecewiselinear (PWL) approximation technique proposed by Juli´an et al. in [3, 4, 5, 6]. In [7], we faced the problem of the PWL approximation of dynamical systems of the kind
N˙ = B (N; F)
(1)
where the (continuous) flow B : S ⊂ 4(n+q) → 4n , the state vector N ∈ 4n , and the parameter vector F ∈ 4q are column vectors (n and q are positive integer numbers) and S is a limited domain. The approximation of the flow B by a continuous PWL function BPWL was obtained with minimization techniques starting from a set of NS samples of B properly distributed over the domain S. In particular, by resorting to continuation methods [8] applied to a smoothed version of BPWL , we showed that the PWL approximation is structurally stable to parameter variations. The main advantage of the PWL approximation method is that it provides a synthesis technique to find direct circuit implementations of PWL models. This can be useful in many contexts, in particular, in the framework defined in [1, 2]. In this paper, we aim to exploit the PWL approximation technique to identify black-box systems, only some outputs of which are known. Toward this end, we shall consider two examples. In the first, the identification will be derived from a single time series, obtained by sampling a component of the trajectory of a ∗ This work was supported by the European Projects APEREST: IST2001-34893 and OFES-01.0456. † This work was supported by the MIUR, within the PRIN and FIRB frameworks, and by the University of Genoa.
0-7803-8251-X/04/$17.00 ©2004 IEEE
known system (the Colpitts oscillator). Such a trajectory corresponds to a proper set of parameters, and visits a large part of the allowed state space. The main drawback of this approach is that it is difficult to identify systems with multiple coexisting attractors, as each trajectory evolves within the basin of attraction of just one of them, whereas the others remain ”invisible.” The second example concerns a system with two coexisting attractors, and the identification is performed by forcing the trajectory to evolve in the basins of attraction of both of them. Toward this end, we need one or more control parameters to drive the system in all the regions of the parameter space, characterized by qualitatively different dynamical behaviors. By properly modulating such parameters, we can build up a trajectory containing information about the two system attractors. 2. PWL IDENTIFICATION METHOD Given a trajectory component y(t) of a dynamical system described by an (unknown) flow B , we can properly sample it, thus obtaining a time series {˜ yk }(k = 1, . . . , M ). The first step of the PWL identification method is to obtain a reconstruction of the trajectory from which the time series was extracted. Among the possible reconstruction techniques, in this paper we shall apply the standard technique known as Principal Component Analysis (PCA) to a properly scaled version {yk } of {˜ yk }. The time series {yk } has zero mean and unitary variance. Through the PCA, we aim to find the P principal components ξ1 , . . . , ξP (collected in the vector ) corresponding to a flow BPCA able to justify at least 80% of the time series variance. Of course, if the original trajectory is obtained for a fixed parameter configuration (i.e., it is a state trajectory), we expect to find P = n principal components. On the contrary, if the original trajectory is obtained by varying some control parameters, an increased number P of principal components can be found, as the state variables and parameters play the same role from an identification point of view. The PCA-reconstructed trajectory (made up of M × P points obtained by projecting {yk } along the principal components) is then used for the identification of a PWL function BPWL that approximates the reconstructed flow BPCA over a domain SP CA that is scaled with respect to S due to the normalization {˜ yk } −→ {yk }. A significant aspect of the identification procedure concerns the spatial distribution of the fitting samples. Some subregions of the scaled domain SP CA will be densely covered by the reconstructed trajectory samples, whereas other subregions will contain few or no samples at all, depending on the shape of the attractor. In order to find out the function BPWL that identifies best the
IV - 665
ISCAS 2004
➡
➡ 55 6
y
x
3
3.5
40
1 25
−1.5 −4 −81
x
2
10
x
−59
1
−37 −15
55
40
25
10
7 −5
t
−5 0
62.5
125
187.5
250
(b)
(a)
Figure 1: The Colpitts Shil’nikov-like strange attractor (see text for the parameter values) used for the identification: (a) trajectory in the state space; (b) time series of the observed variable y(t) = x2 (t).
ξ3
ξ3
2.6
2.6
ξ
1.4
2
−1.4 0.2
−1 ξ
−0.35
1
−2.2 −1.8
ξ
1.4
2
0.2
−0.6
0.6
1.75 1.8
−0.35
1
−2.2 −1.8
0.7
−1.4
−1 ξ
0.7 −0.6
1.75
0.6
1.8
3 2.8
3 2.8
(a)
(b)
Figure 2: Results of the PWL identification method for the Colpitts oscillator with fixed parameter values: (a) three-dimensional PCA reconstruction of the observed variable y(t) = x2 (t); (b) free trajectory of the PWL identified system. flow BPCA all over SP CA (according to a chosen simplicial partition of it), the samples should be uniformly distributed in the PCA space. Therefore, we shall preprocess the reconstructed state vector by properly interpolating the original M × P samples and by imposing that the flow BPWL at the boundary of SP CA be directed inside to avoid divergency. By doing so, we finally obtain a regular grid of NS samples uniformly distributed in the PCA space. As a final step, such samples are used for the identification of a PWL model of the system by assuming that each dimensional component ξi (i = 1, . . . , P ) is partitioned into m subintervals. Of course, this approach could be greatly improved by using multiple trajectories starting from initial points uniformly distributed over the domain.
3. SYSTEM IDENTIFICATION WITH FIXED PARAMETERS: THE COLPITTS OSCILLATOR The model given by Eqs. (2) describes the Colpitts oscillator, which is widely used in electronic applications [9]:
8> i h >><x˙ = Q(1g− k) −α (e − 1) + x i g h x˙ = −1+x (1 − α )e >> Qk >:x˙ =− Qk(1 − k) hx + x i − 1 x ∗
1
F
∗
2 3
F
g∗
−x2
−x2
1
3
3
2
Q
(2)
3
The model’s dynamics was analyzed in detail in [10]. To apply the first step of the identification procedure described in the previous section, the parameters have been kept at the following values in order to obtain a trajectory wandering on a Shil’nikov-like chaotic attractor: αF = 0.996, k = 0.5, g ∗ = 4.1687, and Q = 1.6596. Figure 1(a) shows the trajectory in the state space, and Fig. 1(b) shows a window of the time series {yk }(k = 1, . . . , M = 105 ) obtained by sampling the second component of the state trajectory by a fixed sampling step ∆t = 2π/100. The PCA provides P = 3 principal components, as expected. Figure 2(a) shows the state trajectory reconstructed through the PCA. The PWL approximation BPWL of the flow BPCA is obtained
IV - 666
➡
➡ 1
p
0.22
B a bc
2
Tc
0.5
p
1
0.11 +
H
A
0
C
DH
0
−0.5
−0.11 −
H −1 −1
p
1
−0.5
0
0.5
1
t
−0.22 0
300
600
900
1200
(b)
(a)
Figure 3: (a) Bifurcation diagram of the Bautin normal form (3) and (b) the parameter forcing signal used for the identification. The coordinates of the three black points a, b, and c are: a = (−0.18, 0.8), b = (−0.04, 0.8), and c = (0.05, 0.8). The curve Tc denotes a fold bifurcation of cycles, whereas the half-lines H + and H − mark Hopf bifurcations, and the point DH corresponds to a degenerate Hopf bifurcation. The dashed line in diagram (a) indicates the parameter region visited by the driving signal shown in plot (b). The dashed line in (b) is the average value p¯1 = −0.02 of the driving signal. over the scaled domain SP CA = { ∈ 43 : a1 ≤ ξ1 ≤ b1 , a2 ≤ ξ2 ≤ b2 , a3 ≤ ξ3 ≤ b3 }, with a1 = −1.25, b1 = 2.61, a2 = −1.81, b2 = 3.03, a3 = −2.25, and b3 = 2.58. The domain SP CA is partitioned by performing m = 6 subdivisions along each principal component. The NS fitting samples needed to find out BPWL are obtained from the M × P points of the reconstructed state vector by using the strategy described in the previous section. The regular grid of NS points is defined by assuming a fitting factor of 8 points per subdivision per dimension (i.e., NS = (8m)3 ). Figure 2(b) shows the state trajectory in the PWL identified system. A comparison of Fig. 1(a) with Fig. 2(a) points out the qualitative similarities among the original, PCA-reconstructed and PWLidentified trajectories. 4. SYSTEM IDENTIFICATION WITH PARAMETER MODULATION: THE BAUTIN SYSTEM In this second example, we aim to apply the PWL identification method in order to identify a system (the Bautin system) admitting coexisting attractors. The Bautin system is described by the following equations: x˙ 1 =p1 x1 − x2 − x1 (x21 + x22 )(x21 + x22 − p2 ) x˙ 2 =p1 x2 + x1 − x2 (x21 + x22 )(x21 + x22 − p2 )
(3)
The bifurcation diagram of system (3) (see [8, 7] for details) is shown in Fig. 3(a) for the sake of clarity. The parameter chosen to drive the trajectory of the system is p1 . In other words, we shall construct our time series by replacing the constant parameter p1 in Eqs. (3) with a time-varying parameter defined as follows:
4 (t mod T ) − p1 (t) = p¯1 + A
T
T 2
+ν
order to let the projection of the trajectory in the parameter plane (p1 , p2 ) oscillate over the three regions A, B, and C (see Fig. 3(a)), characterized by the presence of different attractors (an equilibrium in A, an equilibrium and a cycle in B, and a cycle in C). The trajectory in the control space (p1 , x1 , x2 ) is shown in Fig. 4(a), and the time series {yk } (containing about M = 125000 samples) of the observed variable y(t) = x1 (t) is presented in Fig. 4(b). A comparison of this figure with Fig. 3(b) points out that the frequency of the driving signal p1 (t) is slow, as compared with the trajectory dynamics. This fact, together with the relatively small number (100) of correlation samples we used for the PCA, leads the reconstruction process to find out only P = 2 principal components, corresponding to the state variables. The control parameter p1 (t) is added to the PCA space as a third dimension before the PWL identification step. The PWL approximation BPWL of the flow BPCA is obtained over the scaled domain SP CA = { ∈ 43 : a1 ≤ ξ1 ≤ b1 , a2 ≤ ξ2 ≤ b2 , a3 ≤ ξ3 ≤ b3 }, with a1 = −2.1, b1 = 2.1, a2 = −2.1, b2 = 2.1, a3 = −0.2, and b3 = 0.16. The domain SP CA is partitioned by performing m = 5 subdivisions along each principal component. The NS fitting samples needed to find out BPWL are obtained from the PCA-reconstructed trajectory by assuming a fitting factor of 9 points per subdivision per dimension (i.e., NS = (9m)3 ). Figure 5 shows the three phase portraits of the identified PWL system for parameters fixed at the points a, b, and c in the (p1 , p2 ) plane (see Fig. 3(a)). The results evidence the good identification performed by the described method: in the first figure (point a), there is one stable equilibrium point only, the second figure (point b) points out the coexistence of a stable equilibrium point and two cycles (one stable and one unstable), whereas, in the third figure (point c), only one stable cycle is present. 5. CONCLUDING REMARKS
(4)
2000 where p¯1 = −0.02, A = 0.18, T = and ν is white Gaus3 sian noise with standard deviation σ = 2.5 10−3 (see Fig. 3(b)). The second parameter of the Bautin system is fixed at p2 = 0.8 in
This paper can be viewed as a first step towards the definition of a method for PWL modeling/identification of nonlinear black-box systems. To sum up, we have applied the PWL approximation method proposed by Juli´an et al. in [3, 4, 5, 6] to two academic identification problems. Even though the two considered systems
IV - 667
➡
➠ 1.2
x
2
y
0.6
1.2 0.6
0 0
p 1
−1.2 −0.2
1
1.2 0.6
−0.6 x
−0.1
0
0 −0.6 0.2 −1.2
0.1
−0.6
t
−1.2 0
300
600
(a)
900
1200
(b)
Figure 4: The Bautin normal form used for the identification (see text for the description of the forcing parameter): (a) trajectory in the control space; (b) time series of the observed variable y(t) = x1 (t).
2 ξ
ξ2
2
ξ2
1
0
−1
ξ1
−2 −2
−1
0
1
2
ξ1 −2
−1
0
1
2 −2
ξ1 −1
0
1
2
Figure 5: Result of the identification procedure for the modulated Bautin normal form: phase portraits of the identified PWL system corresponding to the parameter values at points a (first figure), b (second figure), and c (third figure) in Fig. 3(a), respectively. reduction,” IEEE Transactions on Circuits and Systems—I, vol. 47, pp. 702–712, 2000.
were perfectly known, no information about them was used, except for the observed time series. The reported results point out a good identification of both systems.
[5] M. Storace, P. Juli´an, and M. Parodi, “Synthesis of nonlinear multiport resistors: a PWL approach,” IEEE Transactions on Circuits and Systems—I, vol. 49, pp. 1138–1149, 2002.
6. REFERENCES
[6] M. Storace, L. Repetto, and M. Parodi, “A method for the approximate synthesis of cellular nonlinear networks - Part 1: Circuit definition,” International Journal of Circuit Theory and Applications, vol. 31, pp. 277–297, 2003.
[1] O. De Feo and M. Hasler, “Modeling and calssification of approximately periodic signals using chaotic systems,” in International Conference on Nonlinear Theory and its Applications N OLT A, (Xi’an, PRC), pp. 11–15, October 7-11 2002.
[7] M. Storace and O. De Feo, “PWL approximation of dynamical systems: An example,” in International Symposium on Circuits and Systems ISCAS, (Bangkok, Thailand), pp. 654–657, May 2003.
[2] O. De Feo, Modeling Diversity by Strange Attractors with Application to Temporal Pattern Recognition. PhD thesis, Swiss Federal Institute of Technology Lausanne, Lausanne, Switzerland, 2001. [3] P. Juli´an, A. Desages, and O. Agamennoni, “High-level canonical piecewise linear representation using a simplicial partition,” IEEE Transactions on Circuits and Systems—I, vol. 46, pp. 463–480, 1999. [4] P. Juli´an, A. Desages, and B. D’Amico, “Orthonormal high level canonical PWL functions with applications to model
[8] Y. Kuznetsov, Elements of Applied Bifurcation Theory. New York, NY: Springer-Verlag, 2nd ed., 1998. [9] S. Sedra and K. Smith, Microelectronic Circuits. New York, NY: Oxford University Press, 4th ed., 1998. [10] O. De Feo, G. Maggio, and M. Kennedy, “The Colpitts oscillator: Families of periodic solutions and their bifurcations,” Int. J. Bif. & Chaos, vol. 10, pp. 935–958, 2000.
IV - 668