- Email: [email protected]
Implementation of Input Design for Structured Nonlinear System Identification via Linear IMC Control
Proceedings of the 15th IFAC Symposium on System Identification Saint-Malo, France, July 6-8, 2009
Implementation of Input Design for Structured Nonlinear System Identification via Linear IMC Control Carlo Novara ∗ Tyrone L. Vincent ∗∗ Kameshwar Poolla ∗∗∗ ∗
Dipartimento di Automatica e Informatica, Politecnico di Torino, Italy. email: [email protected] ∗∗ Division of Engineering, Colorado School of Mines, Golden, CO 80401. email: [email protected] ∗∗∗ Dept. of Mechanical Engineering, University of California, Berkeley, CA 94720. email: [email protected]
Abstract: Structured nonlinear system identification is a block oriented identification method for systems that can be described by the interconnection of known linear dynamic systems and unknown static nonlinearities. An identification algorithm for this class of systems was previously presented in Hsu et al. (2008), and the related input design problem was solved in Vincent et al. (2009) when there is no feedback from the nonlinear blocks. In the present paper, we propose a strategy, based on linear Internal Model Control (IMC), for the implementation of the input design on a general structure. The strategy relies on the design of an IMC controller which minimizes a bound on the tracking error, defined as the difference between desired and actual input to the static nonlinearities. An application related to input design for the identification of semi-active suspensions is presented. 1. INTRODUCTION
e1
Specifying complex systems via interconnections of a network of simpler subsystems is a common way of developing a hierarchical model. Since these interconnected subsystems are often unknown or hard to model, identification of the elements in the network is a necessary first step toward control, signal estimation, or fault detection. A theory of system identification of interconnected systems has been developed for the situation where the subsystems consist of known linear dynamic blocks and unknown static nonlinearities in Hsu et al. (2008). Even though the linear block is considered known, this is a very flexible model structure. Note, for example, that NARX models, which are among the most used models for nonlinear systems identification, can be described as structured systems composed of a known linear dynamic block consisting of delays along with an unknown static nonlinearity. Other examples are spray deposition systems (see Hsu et al. (2006)), which are used in the manufacture of photo-voltaics, and the example that will be discussed in Section 5, where identification of a semi-active suspension is performed. An example interconnected system is shown in Figure 1. The elements Li are known time invariant linear dynamic systems while Ni are unknown static nonlinearities to be identified. The problem of identifying the unknown static nonlinearities is called the structured identification problem and is described in Hsu et al. (2008). The data structure used in the structured identification problem is a Linear Fractional Transformation (LFT) e.g. ⋆ Supported in part by NSF under Grants ECS 03-02554 and ECS 01-34132.
978-3-902661-47-0/09/$20.00 © 2009 IFAC
168
u1
u2
L1
N1
N3
L3
L2
N2
y
e2 Fig. 1. Example of an interconnected system Packard and Doyle (1993). As shown in Figure 2, all the linear dynamics are collected into the block L. The static nonlinear functions are gathered into the block N .
Fig. 2. LFT model structure. The linear block L is partitioned comformably as Lyu Lye Lyw . L= Lzu Lze Lzw
The nonlinear block N has a block-diagonal structure that we represent as
10.3182/20090706-3-FR-2004.0320
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009
N =
N1
..
.
.
Nn This notation signifies a specific association of some components of z as inputs to nonlinearity Ni . Without loss of generality, we assume that each component Nk is singleoutput. The case when the signal z, which is the input to the nonlinear block, is either directly measurable or can be reconstructed from knowledge of u and y alone has been shown to yield a particularly tractable identification problem in Hsu et al. (2008). However, this assumes that z satisfies appropriate conditions of persistence of excitation. In order to obtain suitable data, some care must be taken in the generation of an appropriate input for the experiment. This is the subject of this paper. Throughout the development of system identification theory, the experiment design problem has remained prominent, e.g. Mehra (1974b); Goodwin and Payne (1977); Ljung (1999); Rivera et al. (2003); Gevers (2005); Jansson and Hjalmarsson (2005); Rolain et al. (2006); Wahlberg et al. (2006). A typical approach for input design for linear dynamic systems when utilizing a parameterized model structure is to consider either maximizing a suitable measure of the Fisher information matrix or minimizing a measure of its inverse, which defines the CramerRao lower bound on the parameter error covariance for unbiased estimators. This optimization is done subject to constraints such as maximum amplitude or power of the input sequence. Since the trace of the information matrix is a quadratic function of the input trajectory for linear parameterizations, the input design can be posed as an linear-quadratic optimal control problem as in Mehra (1974a). However, it is usually preferred to minimize the Cramer-Rao bound directly, which is non-convex in terms of the input trajectory. An alternative is then to utilize a two step method, by which the Cramer-Rao bound is optimized with respect to the input distribution or spectrum, and then generate an input that achieves the desired distribution or spectrum in a separate step, e.g. Mehra (1974b); Jansson and Hjalmarsson (2005). Input design for nonlinear system identification presents special challenges, as both spatial and frequency content of input signals must be considered (Leontaritis and Billings (1987)). Kalaba and Spingarn (1977); Espie and Macchietto (1989) examined extensions of the optimal control approach to nonlinear systems, but this is practical only for systems with a few parameters. In Hjalmarsson and M˚ artensson (2007), optimization of the Cramer-Rao bound and the use of two step process, similar to the linear case, is suggested. For nonlinear systems higher order moments of the input are needed to characterize the bound, however it is shown in Hjalmarsson and M˚ artensson (2007) that by assuming a Gaussian input sequence with the autocorrelation as the design parameters, or by using a periodic deterministic sequence with the sequence values as the design parameters, the design problem can be expressed via a polynomial matrix inequality, and sum of squares relaxations used to find the solution. Another line of investigation relevant to the Hammerstein/Wiener model structure is the use of so called separable inputs
169
as in Billings and Fakhouri (1982); Hunter and Korenberg (1986); Enqvist and Ljung (2005), which have properties that allow the identification of linear and nonlinear blocks to be considered separately. However, since we consider the linear block known, our problem is distinct from Hammerstein/Wiener. There are two distinct issues that need to be addressed when considering input design for the structured nonlinear system identification problem. The first is characterization of an optimal sequence for z, that is, one that is “highly exciting”. The second is the specification of u in order to achieve this z. The problem of developing an optimal sequence u for the special case of Lzw and Lze equal to zero is addressed in Vincent et al. (2009). In the present paper, we propose a solution for the general case based on linear Internal Model Control (IMC) (see e.g. Goodwin et al. (2001)). The strategy is to design the optimal input u∗ using an approximate system, and then implement the input design by means of feedback, which can help to compensate for deviations of the actual plant from the design model. We first derive a sufficient condition for the stability of the feedback system and give a bound on the tracking error, defined as the difference between the optimal signal z ∗ and the actual one. Then, we describe how to design an IMC controller that minimizes this bound. The proposed LFT-IMC implementation strategy is applied to the problem of input design for the identification of a Nonlinear AutoRegressive with eXternal input (NARX) system, representing the semi-active suspension of a vehicle. 2. NOTATION ∆ ξ −1 I tr(X) k·k k · k∞
first difference unit-delay operator the identity trace of matrix X Euclidean vector norm or ℓ2 signal norm Induced ℓ2 operator norm or H∞ norm 3. STRUCTURED ID ALGORITHM
For convenience, the structured identification algorithm of Hsu et al. (2008) and the optimal input design method of Vincent et al. (2009) are summarized. The key to the identification algorithm is a regularization term that is a type of complexity metric for N , called the dispersion function. This term can be described after establishing some notation. Define Nu to be the structured nonlinearity that contains all elements of N without repeats, and let ζ [i] and ω [i] denote input and output sequences for the ith element of Nu . Definition 1. Given length L sequence z, define Πz to be the operator that maps z to ζ = ζ [1] × · · · × ζ [q] , with the elements of ζ sorted by magnitude and channel, and repeated elements removed. The signal ζ contains all of the locations in the input space of Nu that will be visited by z. Since Πz removes
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 only redundant information, there exists another operator, ¯ z such that z = Π ¯ z ζ. which will be denoted Π With Πz available, we can find ω = ω [1] × · · · × ω [q] , the outputs of Nu sorted compatibly with ζ, as ω = Πz w. If w was truly generated by a static nonlinearity, then any elements of w that are eliminated by Πz will also be repeated, as a static function cannot be multi-valued. On the other hand, if these elements are not repeated, we immediately know that w was not generated by a static nonlinearity. The set of signals w that contain repetitions compatible with z can be succinctly expressed as those vectors satisfying the constraint ¯ z Πz w = 0. (1) I −Π
We can now define a measure of complexity for (z, w), the dispersion function, as the sum of squared lengths of the line segments making up the linear interpolant of the scatter plot of (z, w), shifted so that the dispersion of (z, 0) is zero. For all signals z, and w, with w satisfying (1), the dispersion function is
¯ z w 2 D(z, w) = QΠ (2) ¯ = ∆ × · · · × ∆, ∆ is the first difference operator, where Q and k·k is the ℓ2 norm. Since Πz and ∆ are linear, the dispersion function is a quadratic function respect to w, e.g. D(z, w) = wT QT Qw ¯ z. with Q = QΠ
The principal assumptions for application of the identification algorithm are summarized below. A.1 The input-output data (uk , yk )L k=1 is known, and the signal y − Lyu u is bounded. A.2 L is known, Lye = I, and Lyw , Lzu are stable. A.3 z is measurable, i.e. it can be determined uniquely from input/output data alone. A.4 There is no undermodeling. That is, there exists a true nonlinearity N true and true signal etrue that generated the input-output data such that βD(z, N true (z)) < 1 and L1 ketrue k2 < 1. A.5 etrue and zk are independent. k Remark 1. We now offer some commentary on the assumptions. A.1 This assumption is standard. A.2 This assumption is crucial. The assumption Lye = I can be made without loss of generality by prefiltering. A.3 This assumption is essential, and represents a significant restriction on the types of interconnected systems that can be treated. However, we believe it is exactly these systems that can be readily identified. Note that the NARX model structure satisfies this condition. A.4 This assumption is standard. Note that this assumption is expected to hold for the particular L chosen for the experiment. A.5 This assumption is standard. Since z is measurable and known, a choice of N is equivalent to a choice for the sequence w = N (z). The signal e that is consistent with the input/output data follows immediately: e = y − Lyu u − Lyw w. (3)
170
Our preference is for estimates of w that achieve a small error while at the same time achieving a small dispersion. The identification algorithm relies on the optimization of an objective function that balances these two goals: 1 min kek2 + βD(z, w) (4) w L ¯ z Πz )w = 0. s.t. (I − Π A bound on the expected estimate error due to measurement noise is derived in Vincent et al. (2009). This bound can be used to evaluate the quality of the estimate obtained by solving optimization problem (4). For a slightly more compact notation, let H = Lyw . Let N be ¯ z. an orthonormal matrix that spans the nullspace of Az Π The bound is given by noise 2 1 E k˜ ω k ≤ trP −1 , (5) L where 1 ¯′ ′ ¯ ′¯ ′ ¯ Π H H Πz + β Q Q N. P =N L z In Vincent et al. (2009), an input design method is presented, which relies on the minimization of this bound with respect to the sequence z, with the constraint that z must be able to be created by u through Lzu . A technique is also proposed for deriving a sequence u such that Lzu u is equal to this optimal sequence z. The technique can be applied in the special case Lzw = 0 and Lze = 0. In the next section, we propose a solution for the general case Lzw 6= 0 and Lze 6= 0. 4. INPUT DESIGN IMPLEMENTATION As discussed in the previous sections, input design consists in finding an asymptotically realizable sequence z ∗ which allows an optimal (in some sense) estimation of N , and a corresponding sequence u∗ such that z ∗ = Lzu u∗ . However, the problem of implementation arises, since, for a general LFT stucture, the presence of non-zero Lzw and Lze leads always to a deviation of the actual z from the desired sequence z ∗ . In this section, we propose an approach to input design implementation for a general LFT stucture, aimed at reducing these interfering effects. The strategy can be used to implement any input design method and, in particular, the input design method of Vincent et al. (2009) when Lzw 6= 0 and Lze 6= 0.
The approach is to solve the input design problem with an approximate system, and then implement the input design using feedback, which can help to compensate for deviations of the actual plant from the design model. Let KN be a linear approximation to an a-priori estimate of N , which is to be improved by further experimentation. Define ∆N (z) = N (z)−KN z. In this way, ∆N represents a nonlinear deviation from the known linear gain KN . This is illustrated in Figure 3a, where the linear block KN and the nonlinear block N share the same signal path structure. Let L¯ be the system obtained by the feedback connection between L and KN . The subsystems of L¯ are given by L¯ya = Lya + Lyw KN L¯za and L¯za = (I − Lzw KN )−1 Lza , where a stands for u, e or w. By setting L¯zw = L¯ze = 0, the approximate system in Figure 3b can be used for the
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 design of an input sequence u∗ and associated z ∗ = Lzu u∗ leading to an optimal (in some sense) estimation of ∆N (i.e. of N ).
Fig. 4. LFT-IMC input generation. The IMC input generator is the structure inside the dashed lines.
Fig. 3. Approximation for design: (a) Original system with nonlinear block split into known linear approximation KN and unknown deviation ∆N . (b) Approximate system for design after closing the loop containing KN , and setting Lzw = Lze = 0. However, in general, u∗ cannot be used as the input of the system in Figure 3a, because of the interfering effect of Lzw and Lze . The appropriate input sequence could be found by inverting the LFT nonlinear system with respect to the output z. However, the inversion of such a nonlinear system appears to be a very hard problem. A simpler solution, which uses the knowledge of the system structure, is based on linear Internal Model Control (IMC) (see e.g. Goodwin et al. (2001)). The configuration for the IMC-based input generator is shown in Figure 4. A controller Q is chosen to generate u based on u∗ and the error between the desired and measured z. Letting d = L¯zw w ¯ + L¯ze e and w ¯ = ∆N (z), the system ∗ from u to z can be represented by the block diagram of Figure 5, which is a modified version of the typical linear IMC structure. Recall that z ∗ = L¯zu u∗ and let S be the system from d to z (sensitivity transfer function). As evident from Figure 5, if the closed loop system is stable (below we explain which kind of stability) and Q is a right-inverse of L¯zu , then S = 0,
z = z∗,
i.e. the interfering effect of Lzw and Lze is null and the actual sequence z is equal to the optimal sequence z ∗ . However, the exact inversion of L¯zu cannot be performed in many cases, since it may be tall and/or strictly proper and/or with some unstable zero. We thus look for an approximate right-inverse of L¯zu , i.e. for a controller Q such that the tracking error z ∗ − z is “small”.
171
Fig. 5. LFT-IMC block diagram for output z. We now introduce the notion of ℓ2 stable system, present a sufficient condition for the ℓ2 stability of the LFT-IMC system and, under this condition, give a bound on the tracking error. We then describe how to minimize this bound with respect to Q. Definition 2. A nonlinear system Snl : (z, y) = Snl (u, e) is ℓ2 stable if sup kSnl (u, e)k < ∞. k(u,e)k≤1
Let us define the following quantities: γ∆ = sup k∆N (z)k ,
kzk≤1
I − Lzu Q Lzw , γw =
∞ γe = I − Lzu Q Lze ∞ . Theorem 1. Consider the LFT-IMC system of Figure 4. If γw γ∆ < 1, then (i) the LFT-IMC system is ℓ2 stable, and (ii) the tracking error is bounded as
kz ∗ − zk ≤ γw γ∆ Lzu ku∗ k + γe kek . ∞
(6)
Proof.
i) Considering that d = Lzw w + Lze e, the control system of Figure 5 is equivalent to the control system of Figure 6, where a linear dynamic block Lb is connected in feedback with the nonlinear static block ∆N .
The Small Gain Theorem (see e.g. Zhou et al. (1996)) asserts that the overall feedback system is ℓ2 stable if γγ∆ < 1,
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 where γ = kM k∞ and M is the system from w to z. Since this system is I − Lzu Q Lzw , we have that γ = γw and claim i) follows.
these systems allows the design of performing controllers an thus the improvement of the suspensions performances in terms of vehicle comfort and stability.
ii) The output z of the control system of Figure 6 is given by z = Lzu u∗ + I − Lzu Q Lzw w + Lze e . The tracking error is then given by z ∗ − z = − I − Lzu Q Lzw w + Lze e . It follows that
The system considered for input design is a semi-active suspension model identified from experimental data. This model is a simplified version of those obtained in Milanese et al. (2004) from the same experimental data, and is of the NARX form yk = N yk−1 , yk−2 , u1k−1 , u2k−1 , u2k−2 + ek (8)
kz ∗ − zk ≤ γw k∆N (z)k + γe kek
≤ γw γ∆ kzk + γe kek
≤ γw γ∆ Lzu ku∗ k + γe kek . ∞
where yk is the suspension force, u1k is the semi-active damper control current, u2k is the suspension differential speed, and ek is noise. The nonlinearity N is a Nonlinear Set Membership function, see Milanese and Novara (2004), estimated from experimental data. The sampling time is 1/512 s. A NARX system can always be represented in LFT form. In particular, an LFT representation of the suspension system (8) is given by: uk = [u1k u2k ]T yk = wk + ek , T zk = wk−1 + ek−1 , wk−2 + ek−2 , u1k−1 , u2k−1 , u2k−2 , L ] , Lyw = Lye =1, yu = [ 0 0 0 0 1/ξ 0 0 1/ξ 2 1/ξ 0 Lzu = , Lzw = Lze = 0 . 0 1/ξ 0 0 1/ξ 2 0
Fig. 6. LFT-IMC block diagram for output z. Inequality (6) holds for both infinite and finite length signals (a consequence of the fact that the Small Gain Theorem holds for both these classes of signals). In the experiment design setting considered in this paper, we deal with signals of length L. Inequality (6) can thus be expressed as
kz ∗ − zkpow ≤ γw γ∆ Lzu ∞ ku∗ kpow + γe kekpow √ where kλkpow = kλk / L is the power norm (semi-norm for L → ∞). The power (semi-) norm is more suitable than the ℓ2 norm when L is large. Also, when the noise is white and L tends to infinity, kekpow is the standard deviation of the noise e. As a consequence of Theorem 1, it is clear how the controller Q should be designed to minimize the effect of L¯zw and L¯ze . Specifically, the controller to implement, QF , can be found from the the following convex H∞ optimization problem: (7) QF = arg min k(I − L¯zu Q) L¯zw Fw L¯ze Fe k∞ Q
where Fw and Fe are filters that capture the a-priori knowledge of the frequency spectra of w ¯ and e. 5. EXAMPLE This example is about input design for the identification of semi-active suspensions. These suspensions are characterized by a damping coefficient which is controllable by means of a current. The availability of accurate models for
172
According to the implementation strategy of Section 4, a linear approximation KN of N has been obtained from the static characteristic of the suspension provided by the constructor. Considering the feedback connection between L and KN and setting Lzw = Lze = 0, the NARX system (8) can be represented in the approximate LFT form of Figure 3b where ∆N (z) = N (z) − KN z and
Lyu = [ 0 0 ] , Lyw = Lye = 1, Lzu = (I − Lzw KN )−1 Lzu . The following sequence u∗ has been considered for the input generator: T u∗k = [u∗1 u∗2 k = 0, 1, ..., L − 1, k k ] , ∗1 where L = 5000, u is a random current signal uniformly distributed in the interval [0,1.8] A, and u∗2 is a random differential speed signal with zero mean and standard deviation 1 m/s. We can suppose that this sequence had been obtained by means of any input design technique. A white noise sequence e with zero mean and standard deviation 0.1 kN (kilonewton) has also been used to simulate a real experiment context. The IMC-based input generator shown in Figure 4 has been implemented, where Q = QF has been designed solving the optimization problem (7). The filters Fw = (0.08897ξ + 0.08897)2/(ξ − 0.8221)2 and Fe (ξ) = 1 have been used to capture the knowledge on the frequency spectra of the signals w ¯ and e. In particular, the filter Fw has been chosen by analyzing the spectra of the signal KN z obtained from a preliminary simulation of L and assuming a similar spectra for w ¯ = ∆N (z). Note that such a simulation does not require any experiment on the true system.
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 In order to evaluate the performances of the IMC-based input generator, the following indices have been considered: T EN C = kzN C − z ∗ kpow T EIMC = kzIMC − z ∗ kpow where z ∗ = Lzu u∗ is the desired signal to use for the identification of N , zN C = Lzu u∗ + Lze e + Lzw w is the output of the LFT system without feedback control, and zIMC = Lzu (u∗ − QF d) + Lze e + Lzw w is the output of the LFT system with IMC control, see Figures 4 and 5. A simulation has been performed using the signals u∗ and e described above. The values T EN C = 0.034 kN and T EIMC = 0.012 kN have been obtained for the performance indices. These values indicate that an improvement of about 65% in the performance index has been obtained by the IMC-based input generator with respect to the uncotrolled LFT system. Such an improvement can also be observed in Figure 7, where the first component (suspension force) of three signals z ∗ , zN C and zIMC is reported for a part of the simulation. 4 3
z1 (Force) [kN]
2 1 0 −1 −2 −3
0
0.5
1
1.5
Time [s]
Fig. 7. LFT-IMC input generation. First component of z (suspension force). Bold (black): z ∗1 . Dashed (blue): 1 1 zN C . Solid (red): zIMC . REFERENCES S. A. Billings and S. Y. Fakhouri. Identification of systems containing linear dynamics and static nonlinear elements. Automatica, 18(1):15–26, 1982. M. Enqvist and L. Ljung. Linear approximations of nonlinear FIR systems for separable input processes. Automatica, 41(3):459–473, 2005. D. Espie and S. Macchietto. The optimal design of dynamic experiments. AIChE Journal, 35(2):223–229, 1989. Michel Gevers. Identification for control: From early achievements to the revival of experiment design. European Journal of Control, 11:335–352, 2005. G. C. Goodwin and R. L. Payne. Dynamic System Identification: Exeriment Design and Data Analyisis. Academic Press, New York, 1977. Graham C. Goodwin, Stefan E. Graebe, and Mario E. Salgado. Control System Design. Prentice Hall, New Jersey, 2001.
173
H˚ akan Hjalmarsson and Jonas M˚ artensson. Optimal input design for identification of non-linear systems. In Proc. American Control Conference, pages 1572–1576. 2007. Kenneth Hsu, Tyrone L. Vincent, and Kameshwar Poolla. A kernel based approach to structured nonlinear system identification part I: Algorithms. In 14th IFAC Symposium on System Identification, Newcastle, Australia, 2006. Kenneth Hsu, Kameshwar Poolla, and Tyrone L. Vincent. Identification of structured nonlinear systems. IEEE Transactions on Automatic Control, 58(11):2497–2513, 2008. IW Hunter and MJ Korenberg. The identification of nonlinear biological systems: Wiener and Hammerstein cascade models. Biological Cybernetics, 55(2):135–144, 1986. H. Jansson and H. Hjalmarsson. Input design via LMIs admitting frequency-wise model specifications in confidence regions. Automatic Control, IEEE Transactions on, 50(10):1534–1549, 2005. RE Kalaba and K. Spingarn. Optimal input system identification for nonlinear dynamic systems. Journal of Optimization Theory and Applications, 21(1):91–102, 1977. IJ Leontaritis and SA Billings. Experimental design and identifiability for non-linear systems. International Journal of Systems Science, 18(1):189–202, 1987. Lennart Ljung. System Identification: Theory for the User. Prentice Hall, New Jersey, 2nd edition, 1999. R. Mehra. Optimal inputs for linear system identification. Automatic Control, IEEE Transactions on, 19(3):192– 200, 1974a. R. Mehra. Optimal input signals for parameter estimation in dynamic systems–Survey and new results. Automatic Control, IEEE Transactions on, 19(6):753–768, 1974b. M. Milanese and C. Novara. Set membership identification of nonlinear systems. Automatica, 40/6:957–975, 2004. M. Milanese, C. Novara, F. Mastronardi, and D. Amoroso. Experimental modeling of vertical dynamics of vehicles with controlled suspensions. In SAE World Congress, Detroit, Michigan, 2004. A. Packard and J.C. Doyle. The complex structured singular value. Automatica, 29:71–110, 1993. D.E. Rivera, H. Lee, M.W. Braun, and H.D. Mittelmann. Plant-friendly system identification: a challenge for the process industries. 13th IFAC Symposium on System Identification (SYSID 2003), pages 917–922, 2003. Y. Rolain, R. Pintelon, and J. Schoukens. Generation and processing of robust periodic excitations. In 14th IFAC Symposium on System Identification, Newcastle, Australia, pages 1347–1351, 2006. T. Vincent, C. Novara, K. Hsu, and K. Poolla. Experiment design for structured nonlinear system identification. In 15th IFAC Symposium on System Identification, SaintMalo, France, 2009. Bo Wahlberg, H˚ akan Hjalmarsson, and M¨arta Barenthin. On optimal input design in system identification. In 14th IFAC Symposium on System Identification, Newcastle, Australia, pages 499–504, 2006. Kemin Zhou, John C. Doyle, and Keith Glover. Robust and Optimal Control. Prentice-Hall, Upper Saddle River, NJ, 1996.
Implementation of Input Design for Structured Nonlinear System Identification via Linear IMC Control Carlo Novara ∗ Tyrone L. Vincent ∗∗ Kameshwar Poolla ∗∗∗ ∗
Dipartimento di Automatica e Informatica, Politecnico di Torino, Italy. email: [email protected] ∗∗ Division of Engineering, Colorado School of Mines, Golden, CO 80401. email: [email protected] ∗∗∗ Dept. of Mechanical Engineering, University of California, Berkeley, CA 94720. email: [email protected]
Abstract: Structured nonlinear system identification is a block oriented identification method for systems that can be described by the interconnection of known linear dynamic systems and unknown static nonlinearities. An identification algorithm for this class of systems was previously presented in Hsu et al. (2008), and the related input design problem was solved in Vincent et al. (2009) when there is no feedback from the nonlinear blocks. In the present paper, we propose a strategy, based on linear Internal Model Control (IMC), for the implementation of the input design on a general structure. The strategy relies on the design of an IMC controller which minimizes a bound on the tracking error, defined as the difference between desired and actual input to the static nonlinearities. An application related to input design for the identification of semi-active suspensions is presented. 1. INTRODUCTION
e1
Specifying complex systems via interconnections of a network of simpler subsystems is a common way of developing a hierarchical model. Since these interconnected subsystems are often unknown or hard to model, identification of the elements in the network is a necessary first step toward control, signal estimation, or fault detection. A theory of system identification of interconnected systems has been developed for the situation where the subsystems consist of known linear dynamic blocks and unknown static nonlinearities in Hsu et al. (2008). Even though the linear block is considered known, this is a very flexible model structure. Note, for example, that NARX models, which are among the most used models for nonlinear systems identification, can be described as structured systems composed of a known linear dynamic block consisting of delays along with an unknown static nonlinearity. Other examples are spray deposition systems (see Hsu et al. (2006)), which are used in the manufacture of photo-voltaics, and the example that will be discussed in Section 5, where identification of a semi-active suspension is performed. An example interconnected system is shown in Figure 1. The elements Li are known time invariant linear dynamic systems while Ni are unknown static nonlinearities to be identified. The problem of identifying the unknown static nonlinearities is called the structured identification problem and is described in Hsu et al. (2008). The data structure used in the structured identification problem is a Linear Fractional Transformation (LFT) e.g. ⋆ Supported in part by NSF under Grants ECS 03-02554 and ECS 01-34132.
978-3-902661-47-0/09/$20.00 © 2009 IFAC
168
u1
u2
L1
N1
N3
L3
L2
N2
y
e2 Fig. 1. Example of an interconnected system Packard and Doyle (1993). As shown in Figure 2, all the linear dynamics are collected into the block L. The static nonlinear functions are gathered into the block N .
Fig. 2. LFT model structure. The linear block L is partitioned comformably as Lyu Lye Lyw . L= Lzu Lze Lzw
The nonlinear block N has a block-diagonal structure that we represent as
10.3182/20090706-3-FR-2004.0320
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009
N =
N1
..
.
.
Nn This notation signifies a specific association of some components of z as inputs to nonlinearity Ni . Without loss of generality, we assume that each component Nk is singleoutput. The case when the signal z, which is the input to the nonlinear block, is either directly measurable or can be reconstructed from knowledge of u and y alone has been shown to yield a particularly tractable identification problem in Hsu et al. (2008). However, this assumes that z satisfies appropriate conditions of persistence of excitation. In order to obtain suitable data, some care must be taken in the generation of an appropriate input for the experiment. This is the subject of this paper. Throughout the development of system identification theory, the experiment design problem has remained prominent, e.g. Mehra (1974b); Goodwin and Payne (1977); Ljung (1999); Rivera et al. (2003); Gevers (2005); Jansson and Hjalmarsson (2005); Rolain et al. (2006); Wahlberg et al. (2006). A typical approach for input design for linear dynamic systems when utilizing a parameterized model structure is to consider either maximizing a suitable measure of the Fisher information matrix or minimizing a measure of its inverse, which defines the CramerRao lower bound on the parameter error covariance for unbiased estimators. This optimization is done subject to constraints such as maximum amplitude or power of the input sequence. Since the trace of the information matrix is a quadratic function of the input trajectory for linear parameterizations, the input design can be posed as an linear-quadratic optimal control problem as in Mehra (1974a). However, it is usually preferred to minimize the Cramer-Rao bound directly, which is non-convex in terms of the input trajectory. An alternative is then to utilize a two step method, by which the Cramer-Rao bound is optimized with respect to the input distribution or spectrum, and then generate an input that achieves the desired distribution or spectrum in a separate step, e.g. Mehra (1974b); Jansson and Hjalmarsson (2005). Input design for nonlinear system identification presents special challenges, as both spatial and frequency content of input signals must be considered (Leontaritis and Billings (1987)). Kalaba and Spingarn (1977); Espie and Macchietto (1989) examined extensions of the optimal control approach to nonlinear systems, but this is practical only for systems with a few parameters. In Hjalmarsson and M˚ artensson (2007), optimization of the Cramer-Rao bound and the use of two step process, similar to the linear case, is suggested. For nonlinear systems higher order moments of the input are needed to characterize the bound, however it is shown in Hjalmarsson and M˚ artensson (2007) that by assuming a Gaussian input sequence with the autocorrelation as the design parameters, or by using a periodic deterministic sequence with the sequence values as the design parameters, the design problem can be expressed via a polynomial matrix inequality, and sum of squares relaxations used to find the solution. Another line of investigation relevant to the Hammerstein/Wiener model structure is the use of so called separable inputs
169
as in Billings and Fakhouri (1982); Hunter and Korenberg (1986); Enqvist and Ljung (2005), which have properties that allow the identification of linear and nonlinear blocks to be considered separately. However, since we consider the linear block known, our problem is distinct from Hammerstein/Wiener. There are two distinct issues that need to be addressed when considering input design for the structured nonlinear system identification problem. The first is characterization of an optimal sequence for z, that is, one that is “highly exciting”. The second is the specification of u in order to achieve this z. The problem of developing an optimal sequence u for the special case of Lzw and Lze equal to zero is addressed in Vincent et al. (2009). In the present paper, we propose a solution for the general case based on linear Internal Model Control (IMC) (see e.g. Goodwin et al. (2001)). The strategy is to design the optimal input u∗ using an approximate system, and then implement the input design by means of feedback, which can help to compensate for deviations of the actual plant from the design model. We first derive a sufficient condition for the stability of the feedback system and give a bound on the tracking error, defined as the difference between the optimal signal z ∗ and the actual one. Then, we describe how to design an IMC controller that minimizes this bound. The proposed LFT-IMC implementation strategy is applied to the problem of input design for the identification of a Nonlinear AutoRegressive with eXternal input (NARX) system, representing the semi-active suspension of a vehicle. 2. NOTATION ∆ ξ −1 I tr(X) k·k k · k∞
first difference unit-delay operator the identity trace of matrix X Euclidean vector norm or ℓ2 signal norm Induced ℓ2 operator norm or H∞ norm 3. STRUCTURED ID ALGORITHM
For convenience, the structured identification algorithm of Hsu et al. (2008) and the optimal input design method of Vincent et al. (2009) are summarized. The key to the identification algorithm is a regularization term that is a type of complexity metric for N , called the dispersion function. This term can be described after establishing some notation. Define Nu to be the structured nonlinearity that contains all elements of N without repeats, and let ζ [i] and ω [i] denote input and output sequences for the ith element of Nu . Definition 1. Given length L sequence z, define Πz to be the operator that maps z to ζ = ζ [1] × · · · × ζ [q] , with the elements of ζ sorted by magnitude and channel, and repeated elements removed. The signal ζ contains all of the locations in the input space of Nu that will be visited by z. Since Πz removes
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 only redundant information, there exists another operator, ¯ z such that z = Π ¯ z ζ. which will be denoted Π With Πz available, we can find ω = ω [1] × · · · × ω [q] , the outputs of Nu sorted compatibly with ζ, as ω = Πz w. If w was truly generated by a static nonlinearity, then any elements of w that are eliminated by Πz will also be repeated, as a static function cannot be multi-valued. On the other hand, if these elements are not repeated, we immediately know that w was not generated by a static nonlinearity. The set of signals w that contain repetitions compatible with z can be succinctly expressed as those vectors satisfying the constraint ¯ z Πz w = 0. (1) I −Π
We can now define a measure of complexity for (z, w), the dispersion function, as the sum of squared lengths of the line segments making up the linear interpolant of the scatter plot of (z, w), shifted so that the dispersion of (z, 0) is zero. For all signals z, and w, with w satisfying (1), the dispersion function is
¯ z w 2 D(z, w) = QΠ (2) ¯ = ∆ × · · · × ∆, ∆ is the first difference operator, where Q and k·k is the ℓ2 norm. Since Πz and ∆ are linear, the dispersion function is a quadratic function respect to w, e.g. D(z, w) = wT QT Qw ¯ z. with Q = QΠ
The principal assumptions for application of the identification algorithm are summarized below. A.1 The input-output data (uk , yk )L k=1 is known, and the signal y − Lyu u is bounded. A.2 L is known, Lye = I, and Lyw , Lzu are stable. A.3 z is measurable, i.e. it can be determined uniquely from input/output data alone. A.4 There is no undermodeling. That is, there exists a true nonlinearity N true and true signal etrue that generated the input-output data such that βD(z, N true (z)) < 1 and L1 ketrue k2 < 1. A.5 etrue and zk are independent. k Remark 1. We now offer some commentary on the assumptions. A.1 This assumption is standard. A.2 This assumption is crucial. The assumption Lye = I can be made without loss of generality by prefiltering. A.3 This assumption is essential, and represents a significant restriction on the types of interconnected systems that can be treated. However, we believe it is exactly these systems that can be readily identified. Note that the NARX model structure satisfies this condition. A.4 This assumption is standard. Note that this assumption is expected to hold for the particular L chosen for the experiment. A.5 This assumption is standard. Since z is measurable and known, a choice of N is equivalent to a choice for the sequence w = N (z). The signal e that is consistent with the input/output data follows immediately: e = y − Lyu u − Lyw w. (3)
170
Our preference is for estimates of w that achieve a small error while at the same time achieving a small dispersion. The identification algorithm relies on the optimization of an objective function that balances these two goals: 1 min kek2 + βD(z, w) (4) w L ¯ z Πz )w = 0. s.t. (I − Π A bound on the expected estimate error due to measurement noise is derived in Vincent et al. (2009). This bound can be used to evaluate the quality of the estimate obtained by solving optimization problem (4). For a slightly more compact notation, let H = Lyw . Let N be ¯ z. an orthonormal matrix that spans the nullspace of Az Π The bound is given by noise 2 1 E k˜ ω k ≤ trP −1 , (5) L where 1 ¯′ ′ ¯ ′¯ ′ ¯ Π H H Πz + β Q Q N. P =N L z In Vincent et al. (2009), an input design method is presented, which relies on the minimization of this bound with respect to the sequence z, with the constraint that z must be able to be created by u through Lzu . A technique is also proposed for deriving a sequence u such that Lzu u is equal to this optimal sequence z. The technique can be applied in the special case Lzw = 0 and Lze = 0. In the next section, we propose a solution for the general case Lzw 6= 0 and Lze 6= 0. 4. INPUT DESIGN IMPLEMENTATION As discussed in the previous sections, input design consists in finding an asymptotically realizable sequence z ∗ which allows an optimal (in some sense) estimation of N , and a corresponding sequence u∗ such that z ∗ = Lzu u∗ . However, the problem of implementation arises, since, for a general LFT stucture, the presence of non-zero Lzw and Lze leads always to a deviation of the actual z from the desired sequence z ∗ . In this section, we propose an approach to input design implementation for a general LFT stucture, aimed at reducing these interfering effects. The strategy can be used to implement any input design method and, in particular, the input design method of Vincent et al. (2009) when Lzw 6= 0 and Lze 6= 0.
The approach is to solve the input design problem with an approximate system, and then implement the input design using feedback, which can help to compensate for deviations of the actual plant from the design model. Let KN be a linear approximation to an a-priori estimate of N , which is to be improved by further experimentation. Define ∆N (z) = N (z)−KN z. In this way, ∆N represents a nonlinear deviation from the known linear gain KN . This is illustrated in Figure 3a, where the linear block KN and the nonlinear block N share the same signal path structure. Let L¯ be the system obtained by the feedback connection between L and KN . The subsystems of L¯ are given by L¯ya = Lya + Lyw KN L¯za and L¯za = (I − Lzw KN )−1 Lza , where a stands for u, e or w. By setting L¯zw = L¯ze = 0, the approximate system in Figure 3b can be used for the
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 design of an input sequence u∗ and associated z ∗ = Lzu u∗ leading to an optimal (in some sense) estimation of ∆N (i.e. of N ).
Fig. 4. LFT-IMC input generation. The IMC input generator is the structure inside the dashed lines.
Fig. 3. Approximation for design: (a) Original system with nonlinear block split into known linear approximation KN and unknown deviation ∆N . (b) Approximate system for design after closing the loop containing KN , and setting Lzw = Lze = 0. However, in general, u∗ cannot be used as the input of the system in Figure 3a, because of the interfering effect of Lzw and Lze . The appropriate input sequence could be found by inverting the LFT nonlinear system with respect to the output z. However, the inversion of such a nonlinear system appears to be a very hard problem. A simpler solution, which uses the knowledge of the system structure, is based on linear Internal Model Control (IMC) (see e.g. Goodwin et al. (2001)). The configuration for the IMC-based input generator is shown in Figure 4. A controller Q is chosen to generate u based on u∗ and the error between the desired and measured z. Letting d = L¯zw w ¯ + L¯ze e and w ¯ = ∆N (z), the system ∗ from u to z can be represented by the block diagram of Figure 5, which is a modified version of the typical linear IMC structure. Recall that z ∗ = L¯zu u∗ and let S be the system from d to z (sensitivity transfer function). As evident from Figure 5, if the closed loop system is stable (below we explain which kind of stability) and Q is a right-inverse of L¯zu , then S = 0,
z = z∗,
i.e. the interfering effect of Lzw and Lze is null and the actual sequence z is equal to the optimal sequence z ∗ . However, the exact inversion of L¯zu cannot be performed in many cases, since it may be tall and/or strictly proper and/or with some unstable zero. We thus look for an approximate right-inverse of L¯zu , i.e. for a controller Q such that the tracking error z ∗ − z is “small”.
171
Fig. 5. LFT-IMC block diagram for output z. We now introduce the notion of ℓ2 stable system, present a sufficient condition for the ℓ2 stability of the LFT-IMC system and, under this condition, give a bound on the tracking error. We then describe how to minimize this bound with respect to Q. Definition 2. A nonlinear system Snl : (z, y) = Snl (u, e) is ℓ2 stable if sup kSnl (u, e)k < ∞. k(u,e)k≤1
Let us define the following quantities: γ∆ = sup k∆N (z)k ,
kzk≤1
I − Lzu Q Lzw , γw =
∞ γe = I − Lzu Q Lze ∞ . Theorem 1. Consider the LFT-IMC system of Figure 4. If γw γ∆ < 1, then (i) the LFT-IMC system is ℓ2 stable, and (ii) the tracking error is bounded as
kz ∗ − zk ≤ γw γ∆ Lzu ku∗ k + γe kek . ∞
(6)
Proof.
i) Considering that d = Lzw w + Lze e, the control system of Figure 5 is equivalent to the control system of Figure 6, where a linear dynamic block Lb is connected in feedback with the nonlinear static block ∆N .
The Small Gain Theorem (see e.g. Zhou et al. (1996)) asserts that the overall feedback system is ℓ2 stable if γγ∆ < 1,
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 where γ = kM k∞ and M is the system from w to z. Since this system is I − Lzu Q Lzw , we have that γ = γw and claim i) follows.
these systems allows the design of performing controllers an thus the improvement of the suspensions performances in terms of vehicle comfort and stability.
ii) The output z of the control system of Figure 6 is given by z = Lzu u∗ + I − Lzu Q Lzw w + Lze e . The tracking error is then given by z ∗ − z = − I − Lzu Q Lzw w + Lze e . It follows that
The system considered for input design is a semi-active suspension model identified from experimental data. This model is a simplified version of those obtained in Milanese et al. (2004) from the same experimental data, and is of the NARX form yk = N yk−1 , yk−2 , u1k−1 , u2k−1 , u2k−2 + ek (8)
kz ∗ − zk ≤ γw k∆N (z)k + γe kek
≤ γw γ∆ kzk + γe kek
≤ γw γ∆ Lzu ku∗ k + γe kek . ∞
where yk is the suspension force, u1k is the semi-active damper control current, u2k is the suspension differential speed, and ek is noise. The nonlinearity N is a Nonlinear Set Membership function, see Milanese and Novara (2004), estimated from experimental data. The sampling time is 1/512 s. A NARX system can always be represented in LFT form. In particular, an LFT representation of the suspension system (8) is given by: uk = [u1k u2k ]T yk = wk + ek , T zk = wk−1 + ek−1 , wk−2 + ek−2 , u1k−1 , u2k−1 , u2k−2 , L ] , Lyw = Lye =1, yu = [ 0 0 0 0 1/ξ 0 0 1/ξ 2 1/ξ 0 Lzu = , Lzw = Lze = 0 . 0 1/ξ 0 0 1/ξ 2 0
Fig. 6. LFT-IMC block diagram for output z. Inequality (6) holds for both infinite and finite length signals (a consequence of the fact that the Small Gain Theorem holds for both these classes of signals). In the experiment design setting considered in this paper, we deal with signals of length L. Inequality (6) can thus be expressed as
kz ∗ − zkpow ≤ γw γ∆ Lzu ∞ ku∗ kpow + γe kekpow √ where kλkpow = kλk / L is the power norm (semi-norm for L → ∞). The power (semi-) norm is more suitable than the ℓ2 norm when L is large. Also, when the noise is white and L tends to infinity, kekpow is the standard deviation of the noise e. As a consequence of Theorem 1, it is clear how the controller Q should be designed to minimize the effect of L¯zw and L¯ze . Specifically, the controller to implement, QF , can be found from the the following convex H∞ optimization problem: (7) QF = arg min k(I − L¯zu Q) L¯zw Fw L¯ze Fe k∞ Q
where Fw and Fe are filters that capture the a-priori knowledge of the frequency spectra of w ¯ and e. 5. EXAMPLE This example is about input design for the identification of semi-active suspensions. These suspensions are characterized by a damping coefficient which is controllable by means of a current. The availability of accurate models for
172
According to the implementation strategy of Section 4, a linear approximation KN of N has been obtained from the static characteristic of the suspension provided by the constructor. Considering the feedback connection between L and KN and setting Lzw = Lze = 0, the NARX system (8) can be represented in the approximate LFT form of Figure 3b where ∆N (z) = N (z) − KN z and
Lyu = [ 0 0 ] , Lyw = Lye = 1, Lzu = (I − Lzw KN )−1 Lzu . The following sequence u∗ has been considered for the input generator: T u∗k = [u∗1 u∗2 k = 0, 1, ..., L − 1, k k ] , ∗1 where L = 5000, u is a random current signal uniformly distributed in the interval [0,1.8] A, and u∗2 is a random differential speed signal with zero mean and standard deviation 1 m/s. We can suppose that this sequence had been obtained by means of any input design technique. A white noise sequence e with zero mean and standard deviation 0.1 kN (kilonewton) has also been used to simulate a real experiment context. The IMC-based input generator shown in Figure 4 has been implemented, where Q = QF has been designed solving the optimization problem (7). The filters Fw = (0.08897ξ + 0.08897)2/(ξ − 0.8221)2 and Fe (ξ) = 1 have been used to capture the knowledge on the frequency spectra of the signals w ¯ and e. In particular, the filter Fw has been chosen by analyzing the spectra of the signal KN z obtained from a preliminary simulation of L and assuming a similar spectra for w ¯ = ∆N (z). Note that such a simulation does not require any experiment on the true system.
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 In order to evaluate the performances of the IMC-based input generator, the following indices have been considered: T EN C = kzN C − z ∗ kpow T EIMC = kzIMC − z ∗ kpow where z ∗ = Lzu u∗ is the desired signal to use for the identification of N , zN C = Lzu u∗ + Lze e + Lzw w is the output of the LFT system without feedback control, and zIMC = Lzu (u∗ − QF d) + Lze e + Lzw w is the output of the LFT system with IMC control, see Figures 4 and 5. A simulation has been performed using the signals u∗ and e described above. The values T EN C = 0.034 kN and T EIMC = 0.012 kN have been obtained for the performance indices. These values indicate that an improvement of about 65% in the performance index has been obtained by the IMC-based input generator with respect to the uncotrolled LFT system. Such an improvement can also be observed in Figure 7, where the first component (suspension force) of three signals z ∗ , zN C and zIMC is reported for a part of the simulation. 4 3
z1 (Force) [kN]
2 1 0 −1 −2 −3
0
0.5
1
1.5
Time [s]
Fig. 7. LFT-IMC input generation. First component of z (suspension force). Bold (black): z ∗1 . Dashed (blue): 1 1 zN C . Solid (red): zIMC . REFERENCES S. A. Billings and S. Y. Fakhouri. Identification of systems containing linear dynamics and static nonlinear elements. Automatica, 18(1):15–26, 1982. M. Enqvist and L. Ljung. Linear approximations of nonlinear FIR systems for separable input processes. Automatica, 41(3):459–473, 2005. D. Espie and S. Macchietto. The optimal design of dynamic experiments. AIChE Journal, 35(2):223–229, 1989. Michel Gevers. Identification for control: From early achievements to the revival of experiment design. European Journal of Control, 11:335–352, 2005. G. C. Goodwin and R. L. Payne. Dynamic System Identification: Exeriment Design and Data Analyisis. Academic Press, New York, 1977. Graham C. Goodwin, Stefan E. Graebe, and Mario E. Salgado. Control System Design. Prentice Hall, New Jersey, 2001.
173
H˚ akan Hjalmarsson and Jonas M˚ artensson. Optimal input design for identification of non-linear systems. In Proc. American Control Conference, pages 1572–1576. 2007. Kenneth Hsu, Tyrone L. Vincent, and Kameshwar Poolla. A kernel based approach to structured nonlinear system identification part I: Algorithms. In 14th IFAC Symposium on System Identification, Newcastle, Australia, 2006. Kenneth Hsu, Kameshwar Poolla, and Tyrone L. Vincent. Identification of structured nonlinear systems. IEEE Transactions on Automatic Control, 58(11):2497–2513, 2008. IW Hunter and MJ Korenberg. The identification of nonlinear biological systems: Wiener and Hammerstein cascade models. Biological Cybernetics, 55(2):135–144, 1986. H. Jansson and H. Hjalmarsson. Input design via LMIs admitting frequency-wise model specifications in confidence regions. Automatic Control, IEEE Transactions on, 50(10):1534–1549, 2005. RE Kalaba and K. Spingarn. Optimal input system identification for nonlinear dynamic systems. Journal of Optimization Theory and Applications, 21(1):91–102, 1977. IJ Leontaritis and SA Billings. Experimental design and identifiability for non-linear systems. International Journal of Systems Science, 18(1):189–202, 1987. Lennart Ljung. System Identification: Theory for the User. Prentice Hall, New Jersey, 2nd edition, 1999. R. Mehra. Optimal inputs for linear system identification. Automatic Control, IEEE Transactions on, 19(3):192– 200, 1974a. R. Mehra. Optimal input signals for parameter estimation in dynamic systems–Survey and new results. Automatic Control, IEEE Transactions on, 19(6):753–768, 1974b. M. Milanese and C. Novara. Set membership identification of nonlinear systems. Automatica, 40/6:957–975, 2004. M. Milanese, C. Novara, F. Mastronardi, and D. Amoroso. Experimental modeling of vertical dynamics of vehicles with controlled suspensions. In SAE World Congress, Detroit, Michigan, 2004. A. Packard and J.C. Doyle. The complex structured singular value. Automatica, 29:71–110, 1993. D.E. Rivera, H. Lee, M.W. Braun, and H.D. Mittelmann. Plant-friendly system identification: a challenge for the process industries. 13th IFAC Symposium on System Identification (SYSID 2003), pages 917–922, 2003. Y. Rolain, R. Pintelon, and J. Schoukens. Generation and processing of robust periodic excitations. In 14th IFAC Symposium on System Identification, Newcastle, Australia, pages 1347–1351, 2006. T. Vincent, C. Novara, K. Hsu, and K. Poolla. Experiment design for structured nonlinear system identification. In 15th IFAC Symposium on System Identification, SaintMalo, France, 2009. Bo Wahlberg, H˚ akan Hjalmarsson, and M¨arta Barenthin. On optimal input design in system identification. In 14th IFAC Symposium on System Identification, Newcastle, Australia, pages 499–504, 2006. Kemin Zhou, John C. Doyle, and Keith Glover. Robust and Optimal Control. Prentice-Hall, Upper Saddle River, NJ, 1996.