Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Identifying Position-Dependent Mechanical Systems: A Modal Approach Applied to a Flexible Wafer Stage

Identifying Position-Dependent Mechanical Systems: A Modal Approach Applied to a Flexible Wafer... Identifying Position-Dependent Mechanical Systems: A Modal Approach Applied to a Flexible Wafer Stage Robbert Voorhoeve, Robin de Rozario, Wouter Aangenent, and Tom Oomen, Senior Member, IEEE Abstract—Increasingly stringent performance requirements for motion control necessitate the use of increasingly detailed models of the system behavior. Motion systems inherently move, there- Point of interest fore, spatio-temporal models of the flexible dynamics are essen- tial. In this paper, a two-step approach for the identification of Measured output the spatio-temporal behavior of mechanical systems is developed and applied to a lightweight prototype industrial wafer stage. The proposed approach exploits a modal modeling framework and combines recently developed powerful linear time invari- ant (LTI) identification tools with a spline-based mode-shape Fixed world Flexible stage interpolation approach to estimate the spatial system behavior. The experimental results for the wafer stage application confirm the suitability of the proposed approach for the identification of complex position-dependent mechanical systems, and its potential Fig. 1. Schematic flexible wafer-stage system. for motion control performance improvements. Index Terms—system-identification, precision mechatronics include over-actuation and over-sensing [6], spatial vibration control [7], multivariable robust control [1], and inferential I. I NTRODUCTION control of unmeasured performance variables [8]. Invariably, Increasingly stringent performance requirements for pre- these approaches are characterized by an increased reliance cision motion systems lead to a situation where the flex- on model-based control design procedures, necessitating the ible dynamics of moving machine components need to be development of control-relevant, efficient, and numerically actively modeled and controlled. Typical examples include reliable identification algorithms capable of dealing with the the wafer stages in lithographic wafer scanners [1], [2]. complex spatio-temporal system behavior [2], [9], [10]. Traditionally, these stages can be accurately approximated The flexible dynamics of these systems in conjunction with as a rigid body in the frequency range relevant for control the fact that these motion systems move lead to position- [3], [4], thereby enabling static decoupling of the rigid-body dependent system behavior [4], [11]. In this paper, mechanical dynamics and subsequent decentralized control design [5]. systems consisting of a single flexible moving body are Furthermore, when this rigid-body approximation is used, considered and deformations are assumed to be small. As an the spatial system behavior follows directly from the stage example, consider the schematic flexible wafer-stage system geometry. Due to increasing accuracy and speed requirements, in Figure 1. Here, the flexible wafer-stage moves in relation the flexible dynamics in future systems can no longer be to the sensors, which are connected to the fixed world. As a neglected and need to be explicitly addressed. Approaches to result of this relative motion, the sensors measure the position address the resulting complex spatio-temporal system behavior at different points on the flexible structure, and therefore the spatial behavior of this flexible system is observed differently This research is supported by ASML Research as part of the TU/e Impulse I depending on the relative location of the sensors. However, program and is part of the research program VIDI with project number 15698, due to the fact that there is only a single moving body and financed by the Netherlands Organization for Scientific Research (NWO). Robbert Voorhoeve, Robin de Rozario, and Tom Oomen are with the because deformations are small, the structural dynamics of the Eindhoven University of Technology, Department of Mechanical Engineering, flexible body do not change as the system moves. Still, as a Control Systems Technology group, Eindhoven, The Netherlands (e-mail result of the position dependency in the measurement system, addresses: r.j.voorhoeve@gmail.com, r.d.rozario@tue.nl, t.a.e.oomen@tue.nl). Wouter Aangenent is with ASML Research Mechatronics, Veldhoven, The the system dynamics as observed by the sensors are no longer Netherlands (e-mail address: wouter.aangenent@asml.com). Linear-Time-Invariant (LTI), necessitating a deviation from the This article has been accepted for publication in IEEE Transactions on standard LTI control design techniques used in the context of Control Systems Technology (DOI: 10.1109/TCST.2020.2974140). 2020 IEEE. Personal use of this material is permitted. Permission from high-precision motion systems. IEEE must be obtained for all other uses, in any current or future media, Even though standard linear control design approaches including reprinting/republishing this material for advertising or promotional are no longer sufficient to control the position dependent purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. dynamics, the particular properties of the considered class of arXiv:1807.06942v2 [cs.SY] 28 Feb 2020 2 systems can be exploited to utilize control approaches that remain relatively close to LTI control theory. Approaches in literature that utilize additional system information to retain some aspects of LTI control theory include the gain-scheduling control approach [12]–[14] and the Linear Parameter-Varying (LPV) control framework [15]–[18], which formalizes the gain-scheduling method by ensuring stability and performance through a rigorous model-based mathematical approach. While the class of systems considered in this paper can be readily treated within the LPV control framework, the present paper proposes an approach that explicitly utilizes the additional assumptions of a single moving body to simplify the problem. A key challenge for systematic LPV control is the avail- Fig. 2. The experimental wafer-stage setup. ability of accurate LPV models. The need for accurate LPV models spurred the development of LPV system identification with a strong focus on black-box parametric models [19]–[22]. complexity of such a state-of-the-art industrial system. The This resulted in a well-developed theoretical framework that novelty of this paper therefore lies in the formulation and val- categorizes the identification techniques in local [23], [24], idation of the full position-dependent identification framework, and global approaches, e.g., [25], [26]. Depending on the which, while making use of previous results, including but not application, both approaches have been reported to effectively limited to previous results by the authors, e.g., in [10], [31], support the identification of practically relevant systems [27], has not been previously published. albeit that the identification of systems with high dynamic The identification methods used in this paper have parallels order remains challenging [21], [28]. Due to the high model with the field of experimental modal analysis. Research in complexity associated with general LPV modeling, the success this field has seen significant developments in the past decade of black-box LPV approaches is limited for the identification [32]–[34], in part due to recent consolidation efforts of the of mechanical systems with a large number of resonant modes approaches used in modal analysis and system identification [29], [30], showing a need to reduce the modeling complexity [35]–[37]. Contrary to modal analysis, the modeling goal in by using additional prior system knowledge. this paper is to obtain position-dependent models for control. Although many important developments have progressed This is reflected by the emphasis on accurate mode-shape LPV identification for control, the continuously increasing interpolation and the possibility to incorporate control-relevant complexity of motion systems necessitates a practical iden- identification criteria as in, e.g., [2]. tification approach of reduced complexity that is systematic, The outline of this paper is as follows. In Section II, the accurate, and user-friendly. A key step in this paper is to experimental setup is introduced and the control challenges as utilize the knowledge that the system consists of a single well as the position dependent modeling problem for this setup moving body with small deformations to derive a parsimonious are formulated. In Section III, the proposed two-step position- model-set, significantly reducing the modeling complexity dependent identification approach is explained. In Section IV, compared to a full LPV approach. Furthermore, recently the LTI identification approach and the obtained results for developed efficient and reliable LTI identification tools are the prototype wafer stage system are presented. In Section employed to obtain accurate and coherent local models which V, the approach and results for the mode shape interpolation are particularly suited to a subsequent interpolation step to are presented. In Section VI, a discussion is presented on the obtain the desired position-dependent model. The aim of this applicability of the proposed approach for control. In Section paper is to develop an effective and practical approach for the VII the conclusions of this paper are formulated as well as an identification of position-dependent mechanical systems. outlook on ongoing research. The main contributions of this paper are the following. 1) A two-step modal identification framework for position- II. E XPERIM ENTAL S ETUP AND PROBLEM FORMULATION dependent mechanical systems, including In this section the experimental wafer-stage setup and a) a flexible framework of parameterizations and al- related control challenges are outlined. In Section II-C, the gorithms aimed at obtaining accurate modal LTI position-dependent modeling problem is formulated, as is models of complex mechanical systems, considered in this paper. b) an approach for the interpolation of identified mode-shapes to obtain position-dependent models for control of flexible mechanical systems. A. Experimental Setup 2) Application and validation of the developed approach on The experimental setup considered in this paper is the Over- a state-of-the-art industrial wafer stage setup. Actuated-Testrig (OAT), as is shown in Figure 2. The system These contributions are inextricably linked as the proposed is considered here as an analogue to the flexible wafer-stage as two-step position-dependent identification framework is ex- depicted in Figure 1. The system is controlled in six motion plicitly developed with the goal of being able to handle the degrees of freedom and is magnetically levitated having no 3 w z P () c c K() Fig. 4. The LPV standard plant framework. See, for example, [8], [31], [38], for control design approaches for such problems. Fig. 3. Actuator and sensor configuration of the considered experimental setup. The locations of actuators used for control and identification are marked The LPV standard plant framework, as depicted in Figure with red crosses ( ) actuators used only for identification as blue crosses ( ), 4, can be used to describe both these control problems. sensors used for identification and control as red circles ( ), and sensors used Here, u and y are the output and input signals available only for identification as green circles ( ). c c for control, w is the generalized disturbance signal and z p p is the generalized performance signal, which in this case involves the positioning error of the point of interest. This mechanical connection to the fixed world. Furthermore, the control problem, with the LPV standard plant P () in a experimental system is equipped with additional actuators generalized feedback interconnection with an LPV controller and sensors to facilitate the spatio-temporal identification of K (), has been studied extensively in, e.g., [15]–[17] and, for the flexible dynamics. While these additional actuators and appropriately bounded sets of scheduling variable trajectories, sensors are available on the experimental setup, they are not i.e., (t) 2 D , efficient algorithms exist for various robust and available in the considered wafer-stage system as depicted in optimal control problems defined in this framework. Figure 1 for which the experimental system is an analogue. The main difficulty for the practical application of these Therefore, these additional actuators and sensors are only used methods concerns the availability of an accurate LPV system- for identification and are considered to be unavailable for the model G(), which is the part of the standard plant P () control of the system. In this paper, only the out-of-plane pertaining to the physical system that is to be controlled. That motions are considered, i.e., the motions perpendicular to the this is a difficult problem is evidenced by the fact that accurate surface of the wafer-stage, both for visualization purposes and modeling of LTI precision systems is already considered to be due to the availability of multiple spatially distributed sensors a challenging problem, see, e.g., [9], [10]. Accurate modeling in this direction. for LPV systems is generally significantly more challenging, The out-of-plane sensors and actuators used for the exper- since the incorporation of a scheduling parameter dependency iments in this paper are shown in Figure 3. Seven actuators, typically leads to a highly increased model complexity [21]. depicted by crosses in Figure 3, four of which, shown in Due to the additional constraints on the considered class of red, are used for closed-loop control and the remaining three, systems, consisting of one flexible moving body and assum- shown in blue, are used to apply additional spatially distributed ing small deformations, the complexity of the problem can excitation for identification. Sixteen sensors are available for be significantly reduces. This problem of obtaining accurate identification, as shown in Figure 3 by circles. Similarly, a position-dependent models for such mechanical systems, such distinction is made between the sensors used for closed-loop as the wafer-stage example considered here, is addressed in control, shown in red, and those only used for the spatially the remainder of this paper. distributed identification, which are shown in green. C. Position-Dependent Modeling Problem B. Control Challenges Consider the schematic wafer-stage setup shown in Figure In the considered wafer-stage example, the scheduling vari- 1. For this setup, two key control challenges are recognized ables  relate to the position of the wafer stage. Due to the which are related to the position-dependent system behavior. fact that there is only a single moving body, the structural First, as previously outlined, the relative motion of the wafer dynamics of the flexible body are invariant to the position stage with respect to the sensors leads to an input-output be- of the wafer stage with respect to the sensors. Furthermore, havior that is no longer LTI, necessitating a position/dependent deformations are considered to be small in the sense that the modeling and control perspective. Second, the point of interest, structure does not deform in such a way that it influences i.e., the point on the wafer that is being exposed in the photo- the dynamics of the systems. As a result, the A matrix of lithographic process, also changes as the wafer stage moves. the state-space description for this class of systems does not This involves the control of an unmeasured performance vari- depend on the scheduling variable , see, e.g., [39], [40]. able since the point of interest is not directly measured. This Furthermore, for the considered system the positions of the results in a position-dependent inferential control problem. actuators are fixed relative to the wafer-table and therefore the 4 input matrix B is also not position dependent. Therefore, the that is applicable for any geometrically complex domain D is scheduling variables have no influence on the system states the Finite Element Method (FEM). This approach uses many and as a result there is no memory in the system pertaining to localized basis functions to accurately approximate the spatial the past-trajectory of the scheduling variables. The system is system behavior [41]. therefore only dependent on the current, instantaneous, value The temporal system behavior for this basis function ap- of the scheduling parameters. In state-space the considered proach is determined by the dynamics of the generalized system-model is described as, coordinates, q(t) = [q (t) : : : q (t)] . Under the assumption 1 n of small rotations and strains, and assuming that the material is linear elastic obeying Hooke’s Law, the equations of motion x _ (t) = Ax(t) + Bu(t) (1) that govern the temporal input-output behavior of a mechanical G() : y(t) = C ((t))x(t) + D ((t))u(t) (2) y y system are given by the set of coupled second order ordinary z(t) = C ((t))x(t) + D ((t))u(t): (3) z z differential equations [39, Section 2.2] In the wafer-stage example, both the system outputs y(t) Mq (t) +Dq_(t) +Kq(t) = Qu(t) ; (7) and the performance variables z(t), i.e., the point of interest n n n n position, can be considered as specific local instances of the q q q q where M 2 R is the mass matrix, D 2 R the n n out-of-plane deflection of the surface of the wafer stage. In this q q damping matrix, K 2 R the stiffness matrix, and Q 2 interpretation, the relevant output for this system is this out-of- n n q u R the input distribution matrix. plane deflection, which is denoted here as z(%; t), where % is The set of coupled equations of motion (7) can be decoupled the in-plane coordinate of the point for which the deflection is for the undamped case by transforming to a modal descrip- considered. The modeling problem is then to identify from tion, which is obtained by solving the generalized eigenvalue experimental data of the system, a model for the system problem behavior of the form K ! M  = 0 ; i = 1; : : : ; n : (8) i q x _ (t) = Ax(t) + Bu(t) (4) i G(%) : z(%; t) = C(%)x(t) + D(%)u(t) : (5) 2 The eigenvalues, ! , are the squared undamped resonance- frequencies of the modes, and the eigenvectors,  , are the The model in the form (1)–(3) is recovered from this descrip- i associated mode shapes as parameterized in the basis W (%) = tion by including the static geometric relations between the position of the wafer stage  and the coordinates %, at which [w (%) : : : w (%)]. By applying the substitution q = , 1 n 1 1 the sensors view the wafer-stage as well as the coordinate where  = [ : : :  ], and multiplying (7) with  M 1 n of the point of interest. This is a simple affine relation yields dependent on the definitions of the origins for two coordinate I (t) +D  _ (t) + (t) = Ru(t) ; (9) systems and the sensor locations, i.e., for a certain sensor y , G (%) : z(%; t) = L(%)(t) ; (10) C ((t)) = C(% ((t)) with % ((t)) = % + (t). The y y y y ;0 1 1 1 1 modeling problem considered in the remainder of this paper 1 1 2 2 2 where D =  M D, = diag([! : : : ! ]), 1 n is therefore to model the system of the form (4)–(5). 1 1 L(%) = W (%), and R =  M Q. In the context of identification for control, low-order models III. POSITION-D EPENDENT M ODELING APPROACH are desired that are accurate in a limited frequency band of In this section the proposed position-dependent modeling interest [2]. This means that only a limited number of modes, approach is described. First, the modal modeling framework n < n , are required to model the relevant temporal system m q for mechanical systems is outlined. Next, the proposed two- behavior [39]. Modeling the spatial system behavior, using step identification approach is developed. a generic set of basis functions L(%) = W (%), typically requires a large number of basis functions leading to a high modeling complexity. In this paper, a two-step identification A. Modal Models of Mechanical Systems approach is proposed to directly identify the mode-shapes, i.e., The quantity of interest is the out-of-plane deflection, the columns of L(%) from measured data. z(%; t) : D T 7! R, of the surface of a flexible body, which is modeled here as a continuum. The domain D 2 R of the coordinate % is the surface of the considered structure, e.g., B. Two-Step Identification Approach the (x; y) surface of the wafer-stage. To model the spatio- To identify the spatio-temporal system behavior, measure- temporal evolution of z(%; t), a basis-function expansion is ment data is first obtained in experiments with fixed sensor used for time–space separation, see, e.g., [41], i.e., locations % . As a result of the fixed sensor locations the input-output system behavior is linear time invariant, similar z(%; t) = w (%)q (t) : (6) i i to the local approach in LPV identification. Experiments are i=1 performed with various fixed sensor locations covering the For n ! 1 this expansion converges as long as fw (%)g domain D. An LTI system model is then identified from the q i i=1 is a convergent set of functions for the class of continuous obtained experimental data where the model is parameterized functions on the spatial domain D [41]. A widely used method in modal form, i.e., (9)–(10), with n modes. Instead of m 5 nc v the position-dependent output equation (10) the measured, spatially sampled, outputs z (t) are modeled as s u G 2 3 2 3 z(% ; t) L(% ) 1 1 6 7 6 7 . . . . z (t) = = L (t); L  (11) s 4 5 s s 4 5 . . Fig. 5. Closed-loop identification scheme. z(% ; t) L(% ) n n % % n n % m where the parameters in L 2 R are considered as spatially sampled estimates of the mode shapes L(%). identification scheme is shown in Figure 5. A distinction is This first step requires the LTI identification of a complex made between inputs which are used in the control loop u (t) mechanical system with a high model order and many inputs and those that are not used in the control loop u (t). The nc and outputs. The identification of such complex mechanical control inputs u (t) also include the in-plane actuator signals systems requires the use of efficient and numerically reliable that are used to stabilize the in-plane rigid-body modes. The identification approaches, as have been developed and inves- excitation signals used for system identification are the non- tigated in, e.g., [9], [10]. control inputs u (t) and the additive perturbation signals nc In the second step, the spatial mode shapes L(%) are r (t) as shown in Figure 5. The applied excitation signals estimated from the identified parameters in L . In this step, are all random-phase multisines with a flat amplitude spectrum interpolation techniques are used to reconstruct continuous which are successively applied to the single inputs in separate mode shapes based on these spatially sampled estimates. Since experiments. this step involves the interpolation of spatial functions in % and With a total of 16 out-of-plane sensors, 8 control inputs, not of systems that dynamically depended on a scheduling including 4 in the in-plane direction, and 3 non-control inputs, variable , the interpolation pitfalls as shown in [42] are the identification problem involves first identifying a 24 11 avoided. In Section V, a promising robust and physically closed-loop FRF given by motivated interpolation approach is proposed, but several other " # interpolation techniques, which might be more suitable for ~ ~ P ( ) P ( z r k z u k s u s nc ~ c other applications, can be used in this second step of the P ( ) = ; (12) CL k ~ ~ P ( ) P ( u r k u u k c u c nc proposed two-step approach. c In summary, the proposed two-step approach aims to: where the notation P ( ) is used to denote the identified y x k 1) Identify the modal mechanical LTI model given by (9) empirical transfer function estimate (ETFE) from the input and (11), i.e., estimate the parameters in L ; ; D ; s m signal x to output signal y at frequency point . To obtain R. the FRF of the open-loop system G from this closed-loop FRF 2) Estimate the mode-shapes L(%) based on the spatially the following relation is used, sampled mode-shapes L " # h i ~ ~ P P u r u u c u c nc ~ c IV. LTI I DENTIFICATION OF SPATIALLY SAMPLED ~ ~ G( ) = P P ; z r z u s u s nc ~ ~ P P u r u u S YSTEMS nc u nc nc (13) In this section, the first step of the proposed identification where the arguments have been omitted on the right hand approach is outlined, which is the LTI identification of the side of this expression for brevity. In (13) P = 0 and u r nc u spatially sampled system G . P = I , see, e.g., [2, Appendix A] for additional detail u u nc nc on this closed-loop identification approach. By removing the A. Methods in-plane input directions, the 16  7 FRF of the system G , as given by (9) and (11), is obtained. The LTI identification approach considered here involves a number of key aspects. First, a non-parametric identification In this paper, the delays from the hold circuit in the digital approach is considered, aimed at obtaining accurate Frequency measurement environment are first carefully determined using Response Function (FRF) estimates of the spatially sampled a combination of readily available automatic methods and system G . Second, the modal parametrization as used in this visual inspection of the data while manually tweaking the paper is defined. Third, a black-box Matrix Fraction Descrip- delay value. Subsequently the the FRF measurements are tion (MFD) parametrization is employed, which is parameter- compensated for these delays such that the obtained delay- ized such that it closely matches the modal parametrization. compensated FRF can be modeled in continuous time, i.e., Fourth, the identification algorithms used to estimate the by G (s ), with s = j! , and the remaining identification s k k k system models from the measured data are explained. procedure can be performed in the s-domain. For additional 1) Non-Parametric Identification: The non-parametric fre- details on such a pseudo continuous time modeling approach, quency response function for the wafer-stage system is esti- see, e.g., [44], [43, Section 8.5]. Note that when the model is to mated using the robust multisine approach as explained in, e.g., be used in a control context, the importance of determining the [43, Section 3.7]. The rigid-body motions of the system need correct delay increases significantly and one might consider to be controlled for stable operation, therefore all experiments changing to a fully discrete-time z-domain approach which is are performed in a closed-loop configuration. The closed-loop generally better suited to deal with such unknown delays. 6 2) Modal Parametrization: As outlined in Section III, the identification of increasingly complex system where the the modal model for the spatially sampled system G is numerical conditioning becomes an important limiting factor given by (9) and (11), with parameters contained in the for the performance of the identification algorithms. matrices L ; ; D ; and R. The matrices L and R are s m s Furthermore, this general parametrization enables the use 2 2 fully parameterized while = diag(!  ) with parameters of more structured LMFD parameterizations, which enable 2 2 2 T !  = [! : : : ! ] . The damping matrix D is either fully m the parametrization of system with arbitrary McMillan degree 1 n parameterized, which constitutes a general viscous damping n instead of only being able to parameterize systems where model, or, when using the modal damping model, this is equal the degree n is a multiple of the number of outputs p, as to D = diag( ) with  = [ : : :  ] . m;mod 1 n is the case when using the fully parameterized unstructured Using the modal damping model, the set of differential LMFD as in, e.g., [32], [33]. Here, a generic (pseudo)- equations (9) describing the systems temporal behavior be- observable-canonical LMFD parametrization with a quasi- comes fully decoupled, meaning the system can be considered constant degree structure is used, see, e.g., [34], [45], [46]. as a superposition of independently evolving modes [39]. This parametrization is both identifiable, in the sense that it This representation is extensively used in modal analysis and is not over-parameterized, and generic, meaning that it can design as it simplifies the physical interpretation of the modal approximate all proper LTI systems of the given order up to parameters and the incurred modeling errors by assuming arbitrary precision, see [45]. modal damping is generally small for lightly damped sys- This LMFD parametrization is often used for black-box tems, see [39, Section 2.4]. For the wafer-stage example the identification of LTI systems, whereas in this paper the goal modal damping assumption is used to facilitate a parsimonious is to identify spatio-temporal mechanical systems by utilizing parametrization. Note that the modeling framework used in the modal form, i.e., (9) and (11). Therefore, in this paper this paper enables general linear damping models. a number of additional constraints are incorporated in the The complexity of this parametrization is determined by LMFD parametrization such that it more closely resembles the number of modes that are modeled, n . While a limited the mechanical system model. Here, the following properties number of modes usually dominate the behavior in a given are enforced, frequency range of interest, for some applications the com- 1) an even McMillan degree, by taking n = 2 n , x m bined low-frequency compliance effect of unmodeled higher 2) a relative degree r  2, by constraining the column order modes also needs to be taken into account. In such a degrees of the numerator polynomial matrix N (s; ) to case it is relevant to model an additional compliance term, be 2 lower than the corresponding column degrees of e.g., as direct feed-through such as D(%) in (5), to describe the denominator polynomial matrix D(s; ), the quasi-static deformation resulting from each input signal 3) a prescribed number of rigid-body modes n such that rb u(t), see, e.g., [38]. In this paper such a feed-through term n = 2 n poles are located at s = 0, by factoring out 0 rb is not modeled; this can be straightforwardly incorporated in the rigid body dynamics, see [10] for details. the proposed modeling framework when required for a given application. 4) Identification Algorithms: The identification problem The modal parametrization used in this paper is defined by considered here is to find the parameter vector  that min- (9) and (11) with the modal parameters given by, imizes the identification cost, which is in this paper is a h i weighted least-squares cost function, i.e., T 2 = vec : (14) m L R ! 3) MFD Parametrization: For certain identification algo- = arg min V () = "(s ; ) "(s ; ) ; (17) k k rithms, it is necessary that the model parametrization can be k=1 written as a polynomial matrix fraction description. In this pa- where per, a Left Matrix Fraction Description (LMFD) parametriza- tion is used, given by ~ ^ "(s ; ) = W (k) vec G(s ) G(s ; ) ; (18) k k k ^ ^ ^ G(s; ) = D(s; ) N (s; ) ; (15) pqpq with weighting matrix W (k) 2 C . This general cost pq pp ^ ^ where N (s; ) 2 R [s] and D(s; ) 2 R [s] are real function V () includes other commonly used identification polynomial matrices in the Laplace variable s. Furthermore, criteria [10], such as the sample maximum likelihood criterion these polynomial matrices are linearly parameterized with [43], control relevant identification criteria [2], and the input- respect to the parameter vector  2 R using a set of basis output and element-wise weighted criterion used in [47]. functions such that, Minimization of the cost function (17) is a nonlinear least- h i squares optimization problem. Suitable algorithms to solve this ^ ^ vec = (s) = (s) : (16) D(s; ) N (s; ) j j problem are the Sanathanan-Koerner algorithm and the Gauss- j=1 Newton algorithm, or closely related methods such as the This general linear parametrization allows the use of data- Levenberg-Marquart algorithm. These algorithms are defined dependent orthogonal vector polynomials as basis functions, as follows, where the Sanathanan-Koerner algorithm is only p(q+p)1 2 R [s], see, e.g., [9], which is a key aspect for defined for MFD parameterizations, i.e., using (15)–(16). j 7 h0i Algorithm 1 (Sanathanan-Koerner [48]): Let  be given. LM FD: SK LM FD: LM In iteration i = 0; 1; : : : , solve the linear least squares problem Modal: initial estimate Modal: LM hi+1i hii = arg min W (s ;  ) (s ) ; (19) SK k k k=1 with h i hii hii 1 T ^ W (k;  ) = W (k) D(s ;  ) : G (s ) I sk k k q (20) Algorithm 2 (Gauss-Newton [49]): Given an initial estimate 0 5 10 15 20 25 30 35 40 45 50 h0i , compute for i = 0; 1; : : : iteration number hi+1i hii hii hii Fig. 6. Evolution of the cost function during various steps of the LTI =  + arg min J (s ;  ) + "(s ;  ) ; k  k identification approach. k=1 (21) with are beyond the scope of this paper and will be reported @"(s ; ) @ vec(G(s ; )) k k hii elsewhere. Due to the modal damping assumption used in this J (s ;  ) = = W (k) : T T @ @ hii paper for the modal model, and since this modal damping hii (22) assumption is not enforced in the MFD parametrization, the The Gauss-Newton algorithm, and related algorithms such transformation in step 4 of this identification approach is as the Levenberg-Marquardt algorithm, generally provide fast approximate. This approximate transformation is performed monotonic convergence to a minimum of the cost function by calculating the n = 2n pole locations by solving x m V (). However, these algorithms often converge to local- det[D(s; )] = 0, separating the poles into pole pairs such 2 2 minima that are far from optimal, and therefore their per- that (s + p )(s + p ) = s +  s + ! , with  ; ! 2 R for 1;i 2;i i i i formance is strongly dependent on the quality of the ini- i = 1; : : : ; n , and estimating a model of the form, h0i tial estimate  . The Sanathanan-Koerner algorithm on the other hand does not generally converge monotonically, and, ^ G = ; (23) m;trans s +  s + ! if convergent, its stationary points are generally not optima i i=1 of the cost function [50]. However, the Sanathanan-Koerner pq with R 2 R . In this estimation the denominator param- algorithm often yields adequate, albeit suboptimal, results eters are fixed to the values obtained from the pole-pairs of irrespective of the quality of the initial estimate. Therefore, the the LMFD model. This model is fitted based on the FRF data Sanathanan-Koerner algorithm is often used to provide initial using the same cost function as the other identification steps, estimates that are subsequently refined using a gradient based i.e, using (17). The parameters in L and R are then obtained optimization algorithm such as the Gauss-Newton algorithm, from the singular value decomposition of the residue matrices see, e.g., [10]. R = U S V such that i i i 5) Identification Approach: To identify the modal system i 1 1 H model defined by (9), (11) and parameterized by (14) from [L ] = [U ] [S ] ; [R] = [V ] ; i = 1; : : : ; n ; (24) s i i i 1 m 1 i the identified FRF G (s ), the following steps are followed. s k where [X ] and [X ] denote, respectively, the i-th and column 1) Define weighting filters W (k) as in (18). and the j-th row of matrix X . This transformation performs 2) Perform i iterations of the SK algorithm using the SK well for the considered system, as shown by the results in mechanical LMFD model with constraints as defined in Section IV-B. For more general approaches to transform black- Section IV-A3. box models to a gray-box models, see, e.g., [51]. 3) Perform i iterations of the GN or LM algorithm GN for the LMFD model using the parameters SK;min B. Results corresponding to the lowest cost function value during For the identification of the wafer-stage system, the weight- SK iterations as initial estimate. ing function in (18) is chosen as a weighting with the element- 4) Transform the identified LMFD model to an initial wise inverse of the identified FRF G(s ), i.e., estimate for the modal model as defined by parameters k (14). W (k) = diag(vec(j(s )j)) ; (25) inv k 5) Perform a maximum of i iterations of the GN GN;mod or LM algorithm to converge to a optimum of the cost [(s )] = : (26) [G(s )] function as in (17) for the identified modal model. j When considering generally damped modal models, the This choice reflects the goal of minimizing the relative error ^ ~ fourth step of this identification approach can be performed between the model G(s ; ) and the FRF G(s ). For more k k using an exact transformation, i.e., relating two realizations advanced weighting choices that take into account the control of the same system. Details of this exact transformation objective, see [2]. To emphasize the accurate estimation of the V 8 -100 -150 -200 -100 -150 -200 -100 -150 -200 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 Frequency [Hz] Fig. 7. Bode diagrams of the identified FRFs (dotted blue) and fitted modal model (solid red) for a 3 7 part of the full 16 7 identified system. first few lower-frequency resonances the weighting function is (14), which leads to a slight increase in the cost function truncated, i.e., value. In the final step the cost function is again minimized using the monotonically convergent LM algorithm with the W (k) = min (W (k); w ) ; (27) inv max modal parametrization, which in this case only yields a small decrease of the cost function value showing that the initial where w is chosen such that clipping of the weight max modal estimate is already of a high quality. generally occurs only in the high frequency range, after the The small increase in cost function value when transforming first few resonances. the LMFD modal to the modal model is expected, as this step Steps 2-5 of the identification approach as proposed in reduces the model complexity by enforcing modal damping Section IV-A5 are at first only performed for a 3  7 part and through the elimination of computational modes identified of the full 16 7 identified FRF. This is done both to improve in the LMFD model. This is done by visually comparing the computational efficiency and to simplify the implementation identified pole-locations with the resonances of the identified of rigid-body mode constraints, see, e.g., [10]. After the FRFs, similar to the use of stabilization diagrams in experi- successful identification of the modal parameters for the 3 7 mental modal analysis, see, e.g., [32], [33]. The model order system, the model for the remaining 13 outputs is determined of the identified LMFD model is n = 2n = 46 while the x m by estimating the sampled mode-shape parameters in L for model order of the modal model is n = 24. This shows these outputs while all other parameters remain fixed. This that the transformation yields a significant decrease in model is again done by minimizing the cost function (17) for these complexity with a modest increase in cost function value. additional output, which in this case is simply a linear least squares problem. Figure 7 shows the FRF and identified modal model for the In Figure 6, the evolution of the cost function, V () in 37 part of the system on which the identification procedure is (17), is shown during steps 2-5 of the identification approach performed and Figure 8 shows a more detailed view including of Section IV-A5. As can be seen from this figure, the SK the phase of a single element of the FRF and identified modal algorithm in step 2 is not monotonically convergent, but model. These figures show a good agreement between the yields an appropriate initial estimate for subsequent refinement model and the FRF in the low frequency region as well as using the LM algorithm. The LM algorithms used in step for the dominant resonances in the frequency region up to 3 does show monotonic convergence and yields an LMFD about 1 kHz. In the frequency region beyond 1 kHz, there estimate with a cost function value approximately one order are an additional few accurately identified modes but also a of magnitude below that of the initial estimate provided by the number of unmodeled resonances. The identification of these SK approach. In the next step the LMFD model is transformed high frequency modes is not the main focus in this paper since to the modally damped model defined by the parameters the spatial behavior for such high frequency modes is typically Magnitude [dB] 9 such that U , which is generally interpreted as a measure of the bending energy of the functions, exists for all functions in -100 the space. The functions W that minimize (28) are given by -150 W (x; y; #) = # + x# + y# + # G (x; y); (30) s 0 x y j j -200 j=1 180 2 2 2 G (x; y) = r ln(r ); r = (x  x) + (y  y) ; (31) j j j j j see, e.g., [52]. The number of parameters in (30) is n + 3, where the three additional parameters are related to the -90 monomials up to the first degree which represent the set of functions in W for which U = 0, i.e., the kernel of U . To -180 constrain this underdetermined set of equations, the following 1 2 3 10 10 10 three additional constraints are added which make sure that Frequency [Hz] the function space parameterized using the Green’s functions G (x; y) is orthogonal to the space of first order polynomial, Fig. 8. Detailed Bode diagram including phase of the identified FRF (dotted blue) and fitted modal model (solid red) for the (1; 2) element of the identified n n n % % % X X X system as depicted in Figure 7. # = 0 ; # x  = 0 ; # y  = 0 : (32) j j j j j j=1 j=1 j=1 The solution to (28) using (30) and (32) is given by also of a higher spatial frequency, meaning a higher spatial resolution, than what is available, is required to identify the z  = W (x  ; y  ; #) + # ; (33) k s k k k associated mode-shapes. see, e.g., [52]. In explicit matrix-form this yields h i V. M ODE-SHAPE INTERPOLATION T # = X ; (34) z  : : : z  0 1 n 13 In this section, the interpolation of the spatially sampled h i mode shapes, as given by the columns of L , is considered. T with # = , and where # # # # : : : # 0 x y 1 n This interpolation step is performed to obtain a position- % dependent model that is continuous in the spatial variable %. 2 3 " # 1 x  y 1 1 X X + I 6 7 0 G . X = ; X = ; (35) 0 4 5 A. Methods 0 X 33 A popular method for the interpolation of various types of 1 x  y n n % % 2 3 data at arbitrary spatially distributed points is the smoothed G (x  ; y  ) : : : G (x  ; y  ) 1 1 1 n 1 1 thin-plate-spline interpolation approach. The use of thin-plate- 6 . . 7 . . X = : (36) 4 5 . . splines is physically motivated by the fact that the spline functions are derived as the functions that minimize the G (x  ; y  ) : : : G (x  ; y  ) 1 n n n n n % % % % % bending energy of a thin sheet of elastic material. Therefore, For the interpolation of the spatially-sampled mode shapes this approach is particularly well-suited for the considered as identified in L , the points (x  ; y  ) are equal to % , i.e., the s j j j application of interpolating the structural mode-shapes of (x; y) positions of the sensors, and the values for z  are given motion systems which are thin in one dimension relative to by the identified parameters in the columns of L . This inter- the other dimensions, such as the wafer-stage example. polation is carried out independently for each mode shape, i.e., The smoothed thin-plate-spline interpolating function W for each column of L , where for each mode shape a different for a single mode-shape is derived as follows. Given a set of smoothing parameters  is used. These smoothing parameters n points f(x  ; y  ; z  ) 2 R g and a user-defined smoothing % j j j provide a trade-off between robustness to estimation errors in parameter  2 [0;1), find an interpolating functionW 2 W L and interpolation accuracy at the data-points % . In this s j such that, paper, the values of the smoothing parameters are determined using a Leave-One-Out-Cross-Validation (LOOCV) approach, min jW (x  ; y  ) z  j + U; (28) s j j j i.e., the value of  is used which minimizes the LOOCV error. W 2W j=1 For details on this cross-validation approach, see, e.g., [53]. with 1 1 Z Z B. Results U =  W (x; y) dx dy ; (29) In Figure 9, four of the identified flexible mode-shapes are 11 shown. In total 12 mode-shapes are identified including the Here, the function space W is the space of continuously dif- three out-of-plane rigid-body modes. Up to the ninth mode, as ferentiable functions with square-integrable second derivatives shown in Figure 9d, the identified mode shapes agree well with Phase[deg.] Magnitude [dB] 10 Eigenfrequency: 139.23 Hz Eigenfrequency: 437.46 Hz 1.5 1.5 1 1 0.5 0.5 0 0 -0.5 -0.5 -1 -1 -1.5 -1.5 0.5 0.5 0.5 0.5 x x y y -0.5 -0.5 -0.5 -0.5 x x (b) Fifth mode (saddle). (a) Fourth mode (torsion). Eigenfrequency: 506.08 Hz Eigenfrequency: 879.46 Hz 0.6 1.2 0.4 1 0.2 0.8 0 0.6 -0.2 0.4 -0.4 0.2 -0.6 0 -0.8 -0.2 -1 -0.4 -1.2 -0.6 0.5 0.5 0.5 0.5 x x 0 0 0 0 y y -0.5 -0.5 -0.5 -0.5 x x (c) Sixth mode (umbrella). (d) Ninth mode. Fig. 9. Top views and 3-D surface plots of the identified mode shapes with red dots indicating the points of the identified spatially-sampled mode shapes. theoretical mode-shapes for a thin-plate or the mode-shapes By utilizing the identified position-dependent model of the as obtained by means of a Finite-Element-Method analysis of wafer-stage system, G (%), such global and local control the system, a detailed comparison is omitted for brevity. For approaches can be described in the LPV standard plant frame- the higher-order modes the spatial resolution of the sensors is work, as shown in Figure 4, and can be solved using a range of insufficient to accurately reconstruct the smooth mode-shapes. approaches. In this section, several global and local approaches The results in Figure 9 show the viability of the proposed are considered for both feedback and feedforward control. approach to obtain accurate position-dependent models of flexible mechanical systems. In the following section, the A. Global Spatio-Temporal Feedforward Control potential of the proposed approach is discussed in enabling For the wafer-stage example, the problem of global feedfor- various position-dependent control approaches. ward control aims to minimize the error between the reference signal r (t) and the out-of-plane deflection of the surface VI. OUTLOOK FOR POSITION-D EPENDENT M OTION of the wafer stage z(%; t) over the entire spatial domain D. CONTROL APPLICATIONS More precisely, the global approach is aimed at minimizing In this section, several control approaches are explored the following weighted spatial-norm of the error e(%; t) = that are enabled by the availability of accurate position- r (t) z(%; t) dependent models. For flexible motion systems both feedback and feedforward control problems become more complex as Z Z the point that should track the reference is often not directly kek = e (%; t)W (%)e(%; t) d% dt : (37) 2(D ) S measured, such as the point-of-interest of the wafer stage example in Figure 1. Furthermore, the location of this point- of-interest can change over time, increasing the complexity of This problem can be effectively formulated and solved as an the control problem. H optimal control problem, for details see [31]. Approaches to enhance the control performance for such In Figures 10 and 11, simulation results are shown for flexible motion systems can generally be classified as either, the wafer stage example system. Here, Figure 10 shows the 1) global approaches, aimed at preventing or mitigating the reference profile for r as well as the (x; y) coordinates of the flexible deformations in the entire system, such that the point of interest % (t) and Figure 11 shows the local errors deformation related errors are small; or at the point of interest, i.e., e (t) = r (t) z(% (t); t), for a z z z 2) local approaches, aimed at controlling the position of classical feedforward controller, that minimizes the error at the the point-of-interest of the deformed system sensor locations, and the global feedforward approach. These z z y y z z y y 11 % G() z r e u y x z z^ 0.8 % + zy 0.6 0.4 z^ O() 0.2 obs 0 0.005 0.01 0.015 0.02 0.025 0.03 Fig. 12. Position-dependent observer-based feedback control scheme. t [s] Fig. 10. Reference profile r , and the (x; y)-coordinates of the point of interest % over time. In [38], an observer-based inferential control approach is proposed which is especially suited to minimize the influence of disturbance-induced compliant deformations, i.e., the quasi- 0.05 static deformations induced by a locally applied force, often modeled using an additional feedthrough term, see, e.g., [35]. Combined with a moving disturbance source and point-of- -0.05 interest location, the position dependency of this compliant Classic control effect necessitates a position-dependent control approach to -0.1 Global control obtain the desired performance. This position-dependent con- 0 0.005 0.01 0.015 0.02 0.025 0.03 troller can be effectively and intuitively realized by combining t [s] an observer containing a position-dependent system model with an LTI controller that is robust to the remaining position- Fig. 11. Inferential error e at point of interest over time, showing signifi- dependency, as shown in Figure 12. cantly improved performance for the global feedforward controller relative to the classical feedforward controller. VII. C ONCLUSIONS AND OUTLOOK results show that the global approach yields a significantly A. Conclusions improved inferential performance as opposed to the classical This paper provides a general procedure for the identifi- approach. These results and the practical potential of this cation of position-dependent precision mechatronic systems global feedforward approach for future motion systems are which consist of a single flexible moving body and with small enabled by the modeling procedure developed in this paper. deformations. This is an essential step for the control of future high-precision motion systems. A key step in the proposed B. Local Inferential Feedback and Feedforward Control approach is to utilize prior mechanical systems knowledge The general problem of local inferential control aims to as embedded in the modal modeling framework to obtain a directly optimize the performance at the location where per- parsimonious model set. formance is required, such as the point-of-interest in the wafer- In Section IV, a flexible framework of parametrizations stage example, see Figure 1. Both for feedback and feedfor- and identification algorithms is proposed that is especially ward, these approaches can be described in the LPV stan- suited for the identification of modal models of mechanical dard plant framework when an accurately identified position- systems. For the considered state-of-the-art industrial wafer dependent model is available. Local inferential feedforward stage system, with a total of 7 considered inputs and 16 control typically involves the inversion of the system dynamics outputs, the proposed identification approach yields a very from the inputs to the time-varying or parameter-varying accurate modal system model with 12 identified modes. The performance variables which, apart from explicit inversion, spline-based interpolation approach proposed in Section V can be realized by solving an optimal control problem or provides a robust and effective method to reconstruct the using iterative learning control, see, e.g., [18], [54]. Inferential spatial mode shapes, and is successfully applied to reconstruct feedback control, generally requires the use of two degree of 9 of the identified mode shapes. freedom controller structures as opposed to the single degree Potential applications of the proposed position-dependent of freedom controller structure used in traditional feedback modeling approach for control are numerous, including, e.g., control, see, e.g., [8], [38], this can be directly incorporated the use in global spatio-temporal feedforward control and in the standard-plant approach. observer based inferential feedback control. Inferential feedback control is especially relevant when significant disturbance forces are present in the system. In B. Outlook [38], it is shown that a disturbance-observer can be effectively utilized to estimate the inferential performance variable in In this paper, systems are considered that can be written the presence of significant disturbance forces that are non- as (4)–(5). Although this description is less general than (1)– collocated with the actuator forces. In [38], it is also shown (3), it is envisaged that the proposed framework of identify- that in such a case it is essential to include disturbance models ing modal models of mechanical systems and subsequently in the standard-plant description to obtain accurate results. interpolating the spatial system behavior can be extended normalized signals Inferential error e z 12 to be more broadly applicable. In particular, extending the [18] C. Hoffmann and H. Werner, “A survey of linear parameter-varying con- trol applications validated by experiments or high-fidelity simulations,” proposed framework to consider the modeling of interact- IEEE Transactions on Control Systems Technology, vol. 23, no. 2, pp. ing mechanical subsystems provides the ability to model a 416–433, March 2015. variety of relevant mechanical systems. Examples of such [19] B. Bamieh and L. Giarr, “Identification of linear parameter varying models,” International Journal of Robust and Nonlinear Control, vol. 12, interacting mechanical subsystems include the much used H- no. 9, pp. 841–853, 2002. bridge systems as considered in, e.g., [29], [30]. By separately [20] M. Lovera and G. Mercere, ` “Identification for gain-scheduling: a bal- considering the spatio-temporal behavior of the subsystems, anced subspace approach,” in Proceedings of the 2007 American Control Conference, July 2007, pp. 858–863. such as the beam and carriage in an H-bridge system, and [21] J.-W. van Wingerden and M. Verhaegen, “Subspace identification of modeling the full system behavior as a general interconnection bilinear and LPV systems for open- and closed-loop data,” Automatica, of these component models, a more general class of systems vol. 45, no. 2, pp. 372–381, 2009. can indeed be described, including systems with position [22] R. Toth, ´ Modeling and identification of linear parameter-varying sys- tems, ser. Lecture Notes in Control and Information Sciences. Berlin, dependent state-matrices. Validating the practical applicability Germany: Springer-Verlag Heidelberg, 2010, vol. 403. and performance of this approach, as well as the control [23] J. De Caigny, J. F. Camino, and J. Swevers, “Interpolation-based mod- approaches discussed in Section VI, is a subject of ongoing eling of MIMO LPV systems,” IEEE Transactions on Control Systems Technology, vol. 19, no. 1, pp. 46–63, 2011. research. [24] D. Vizer, G. Mercere, O. Prot, E. Laroche, and M. Lovera, “Linear fractional LPV model identification from local experiments: an H - REFERENCES based optimization technique,” in Proceedings of the 52nd Conference on Decision and Control, Florence, Italy, 2013, pp. 4559–4564. [1] M. van de Wal, G. van Baars, F. Sperling, and O. Bosgra, “Multivariable [25] F. Felici, J. van Wingerden, and M. Verhaegen, “Subspace identification H = feedback control design for high-precision wafer stage motion,” of MIMO LPV systems using a periodic scheduling sequence,” Auto- Control Engineering Practice, vol. 10, no. 7, pp. 739–755, 2002. matica, vol. 43, no. 10, pp. 1684–1697, 2007. [2] T. Oomen, R. Van Herpen, S. Quist, M. Van De Wal, O. Bosgra, and [26] J. Goos and R. Pintelon, “Continuous-time identification of periodically M. Steinbuch, “Connecting system identification and robust control for parameter-varying state space models,” Automatica, vol. 71, pp. 254– next-generation motion control of a wafer stage,” IEEE Transactions on 263, 2016. Control Systems Technology, vol. 22, no. 1, pp. 102–118, 2014. [3] R. Munnig Schmidt, G. Schitter, and J. van Eijk, The Design of High [27] D. Turk, J. Gillis, G. Pipeleers, and J. Swevers, “Identification of lin- Performance Mechatronics, 1st ed. Amsterdam, The Netherlands: IOS ear parameter-varying systems: A reweighted ` -norm regularization 2;1 Press, 2011. approach,” Mechanical Systems and Signal Processing, vol. 100, pp. [4] T. Oomen, “Advanced motion control for precision mechatronics: Con- 729–742, 2018. trol, identification, and learning of complex systems,” IEEJ Journal of [28] P. B. Cox, “Towards efficient identification of linear parameter-varying Industry Applications, vol. 7, no. 2, pp. 127–140, 2018. state-space models,” Ph.D. dissertation, Eindhoven University of Tech- [5] A. J. Fleming and K. K. Leang, Design, Modeling and Control of nology, 2018. Nanopositioning Systems, 1st ed., ser. Advances in Industrial Control. [29] M. Groot Wassink, M. van de Wal, C. Scherer, and O. Bosgra, “LPV Switzerland: Springer International Publishing, 2014. control for a wafer stage: beyond the theoretical solution,” Control [6] R. van Herpen, T. Oomen, E. Kikken, M. van de Wal, W. Aangenent, Engineering Practice, vol. 13, no. 2, pp. 231–245, 2005. and M. Steinbuch, “Exploiting additional actuators and sensors for nano- [30] M. Steinbuch, R. van de Molengraft, and A. van der Voort, “Experi- positioning robust motion control,” IFAC Mechatronics, vol. 24, no. 6, mental modelling and LPV control of a motion system,” in Proceedings pp. 619–631, 2014. of the 2003 American Control Conference, vol. 2, Denver, Colorado, [7] S. O. R. Moheimani, D. Halim, and A. J. Fleming, Spatial Control of United States, June 2003, pp. 1374–1379. Vibration: Theory and Experiments, ser. Series on Stability, Vibration [31] R. de Rozario, R. Voorhoeve, W. Aangenent, and T. Oomen, “Global and Control of Systems, Series A. Singapore: World scientific, 2003, feedforward control of spatio-temporal mechanical systems: With ap- vol. 10. plication to a prototype wafer stage,” IFAC-PapersOnLine, vol. 50, [8] T. Oomen, E. Grassens, and F. Hendriks, “Inferential motion con- no. 1, pp. 14 575–14 580, 2017, IFAC 20th Triennial World Congress. trol: identification and robust control framework for positioning an Toulouse, France. unmeasurable point of interest,” IEEE Transactions on Control Systems [32] P. Guillaume, P. Verboven, S. Vanlanduit, H. Van Der Auweraer, Technology, vol. 23, no. 4, pp. 1602–1610, 2015. and B. Peeters, “A poly-reference implementation of the least-squares [9] R. van Herpen, T. Oomen, and M. Steinbuch, “Optimally conditioned complex frequency-domain estimator,” in Proceedings of the 21st Inter- instrumental variable approach for frequency-domain system identifica- national Modal Analysis Conference, Kissimmee, Florida, United States, tion,” Automatica, vol. 50, no. 9, pp. 2281–2293, 2014. 2003, pp. 183–192. [10] R. Voorhoeve, R. de Rozario, and T. Oomen, “Identification for mo- [33] B. Cauberghe, “Applied frequency-domain system identification in the tion control: Incorporating constraints and numerical considerations,” field of experimental and operational modal analysis,” Ph.D. dissertation, in Proceedings of the 2016 American Control Conference, Boston, Vrije Universiteit Brussel, 2004. Massachusetts, United States, July 2016, pp. 6209–6214. [34] J. Vayssettes, G. Mercere, ` P. Vacher, and R. De Callafon, “Frequency- [11] M. M. da Silva, W. Desmet, and H. V. Brussel, “Design of mechatronic domain identification of aircraft structural modes from short-duration systems with configuration-dependent dynamics: Simulation and opti- flight tests,” International Journal of Control, vol. 87, no. 7, pp. 1352– mization,” IEEE/ASME Transactions on Mechatronics, vol. 13, no. 6, 1372, 2014. pp. 638–646, Dec 2008. [35] A. J. Fleming and S. O. R. Moheimani, “Spatial system identification [12] J. S. Shamma and M. Athans, “Gain scheduling: potential hazards and of a simply supported beam and a trapezoidal cantilever plate,” IEEE possible remedies,” IEEE Control Systems Magazine, vol. 12, no. 3, pp. Transactions on Control Systems Technology, vol. 11, no. 5, pp. 726– 101–107, June 1992. 736, 2003. [13] W. J. Rugh and J. S. Shamma, “Research on gain scheduling,” Auto- [36] R. Pintelon, P. Guillaume, and J. Schoukens, “Uncertainty calculation matica, vol. 36, no. 10, pp. 1401–1425, 2000. in (operational) modal analysis,” Mechanical Systems and Signal Pro- [14] D. J. Leith and W. E. Leithead, “Survey of gain-scheduling analysis and cessing, vol. 21, no. 6, pp. 2359–2373, 2007. design,” International Journal of Control, vol. 73, no. 11, pp. 1001– 1025, 2000. [37] E. Reynders, “System identification methods for (operational) modal [15] A. Packard, “Gain scheduling via linear fractional transformations,” analysis: Review and comparison,” Archives of Computational Methods Systems & Control Letters, vol. 22, no. 2, pp. 79–92, 1994. in Engineering, vol. 19, no. 1, pp. 51–124, Mar 2012. [16] P. Apkarian and P. Gahinet, “A convex characterization of gain- [38] R. Voorhoeve, N. Dirkx, T. Melief, W. Aangenent, and T. Oomen, scheduled H controllers,” IEEE Transactions on Automatic Control, “Estimating structural deformations for inferential control: a disturbance vol. 40, no. 5, pp. 853–864, May 1995. observer approach,” IFAC-PapersOnLine, vol. 49, no. 21, pp. 642–648, [17] C. W. Scherer, “LPV control and full block multipliers,” Automatica, 2016, 7th IFAC Symposium on Mechatronic Systems. Loughborough vol. 37, no. 3, pp. 361–375, 2001. University, United Kingdom. 13 [39] W. K. Gawronski, Dynamics and Control of Structures: A Modal Robin de Rozario received the M.Sc. degree in Me- Approach, 1st ed., ser. Mechanical Engineering Series. Springer-Verlag chanical Engineering (cum laude) in 2015 from the New York, 1998. Eindhoven University of Technology, Eindhoven, the [40] A. Preumont, Vibration Control of Active Structures: An Introduction, Netherlands, where he is currently pursuing a Ph.D. 4th ed., ser. Solid Mechanics and Its Applications. Springer, 2018, vol. degree in the Control Systems Technology group. 246. His research interests include system identification [41] H.-X. Li and C. Qi, “Modeling of distributed parameter systems for and learning control for high-performance motion applications—a synthesized review from time–space separation,” Jour- systems. nal of Process Control, vol. 20, no. 8, pp. 891–901, 2010. [42] R. Toth, ´ F. Felici, P. S. C. Heuberger, and P. M. J. Van den Hof, “Discrete time LPV I/O and state space representations, differences of behavior and pitfalls of interpolation,” in Proceedings of the 2007 European Control Conference, Kos, Greece, July 2007, pp. 5418–5425. [43] R. Pintelon and J. Schoukens, System Identification: A Frequency Domain Approach, 2nd ed. Hoboken, New Jersey, United States: Wiley- IEEE press, 2012. [44] C. H. Houpis and G. B. Lamont, Digital Control Systems: Theory, Hardware, Software, 2nd ed., ser. McGraw-Hill series in electrical engineering: Control theory. McGraw-Hill, 1992. [45] K. Glover and J. C. Willems, “Parametrizations of linear dynamical systems: Canonical forms and identifiability,” IEEE Transactions on Automatic Control, vol. 19, no. 6, pp. 640–646, 1974. [46] M. de Mathelin and M. Bodson, “Canonical vs pseudo-canonical forms for the structural and parametric identification of multivariable systems,” in Proceedings of the 1st European Control Conference, vol. 2, Greno- Wouter Aangenent received the M.Sc. degree (cum ble, France, 1991, pp. 1282–1287. laude) and the Ph.D. degree in Mechanical Engineer- [47] R. A. de Callafon, D. de Roover, and P. M. J. Van den Hof, “Multi- ing from the Eindhoven University of Technology, variable least squares frequency domain identification using polynomial Eindhoven, The Netherlands, in 2004 and 2008, matrix descriptions,” in Proceedings of the 35th Conference on Decision respectively. In 2008, he joined the Research Depart- and Control, Kobe, Japan, 1996, pp. 2030–2035. ment of ASML, Veldhoven, The Netherlands, where [48] C. K. Sanathanan and J. Koerner, “Transfer function synthesis as a ratio he currently leads the Mechatronics and Control of two complex polynomials,” IEEE Transactions on Automatic Control, Research Group. vol. 8, no. 1, pp. 56–58, 1963. [49] D. S. Bayard, “High-order multivariable transfer function curve fitting: Algorithms, sparse matrix methods and experimental results,” Automat- ica, vol. 30, no. 9, pp. 1439–1444, 1994. [50] A. H. Whitfield, “Asymptotic behaviour of transfer function sythesis methods,” International Journal of Control, vol. 45, no. 3, pp. 1083– 1092, 1987. [51] G. Mercere, ` O. Prot, and J. A. Ramos, “Identification of parameterized gray-box state-space systems: From a black-box linear time-invariant representation to a structured one,” IEEE Transactions on Automatic Control, vol. 59, no. 11, pp. 2873–2885, Nov 2014. [52] G. Wahba, Spline Models for Observational Data, ser. CBMS-NSF Regional Conference Series in Applied Mathematics. Philadelphia, Pennsylvania, United States: Society for Industrial and Applied Mathe- matics, 1990, vol. 59. [53] G. Wahba and J. Wendelberger, “Some new mathematical methods for variational objective analysis using splines and cross validation,” Tom Oomen (SM’06) received the M.Sc. degree Monthly weather review, vol. 108, no. 8, pp. 1122–1143, 1980. (cum laude) and Ph.D. degree from the Eindhoven [54] R. de Rozario, T. Oomen, and M. Steinbuch, “Iterative learning control University of Technology, Eindhoven, The Nether- and feedforward for LPV systems: Applied to a position-dependent mo- lands. He held visiting positions at KTH, Stockholm, tion system,” in Proceedings of the 2017 American Control Conference, Sweden, and at The University of Newcastle, Aus- Seattle, Washington, United States, May 2017, pp. 3518–3523. tralia. Presently, he is associate professor with the Department of Mechanical Engineering at the Eind- hoven University of Technology. He is a recipient of the Corus Young Talent Graduation Award, the IFAC 2019 TC 4.2 Mechatronics Young Research Award, the 2015 IEEE Transactions on Control Systems Technology Outstanding Paper Award, the 2017 IFAC Mechatronics Best Paper Award, the 2019 IEEJ Journal of Industry Applications Best Paper Award, and recipient of a Veni and Vidi personal grant. He is Associate Editor of the IEEE Control Systems Letters (L-CSS), IFAC Mechatronics, Robbert Voorhoeve received the M.Sc. degree in and IEEE Transactions on Control Systems Technology. He is a member of Mechanical Engineering (cum laude) and Applied the Eindhoven Young Academy of Engineering. His research interests are in Physics in 2013 and his Ph.D. degree in Mechanical the field of data-driven modeling, learning, and control, with applications in Engineering in 2018 from the Eindhoven University precision mechatronics. of Technology, Eindhoven, The Netherlands. His research interests include system identification, iden- tification for advanced motion control, and control of complex mechatronic systems. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Electrical Engineering and Systems Science arXiv (Cornell University)

Identifying Position-Dependent Mechanical Systems: A Modal Approach Applied to a Flexible Wafer Stage

Loading next page...
 
/lp/arxiv-cornell-university/identifying-position-dependent-mechanical-systems-a-modal-approach-xYsJ2f8nl6

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

ISSN
1063-6536
eISSN
ARCH-3348
DOI
10.1109/TCST.2020.2974140
Publisher site
See Article on Publisher Site

Abstract

Identifying Position-Dependent Mechanical Systems: A Modal Approach Applied to a Flexible Wafer Stage Robbert Voorhoeve, Robin de Rozario, Wouter Aangenent, and Tom Oomen, Senior Member, IEEE Abstract—Increasingly stringent performance requirements for motion control necessitate the use of increasingly detailed models of the system behavior. Motion systems inherently move, there- Point of interest fore, spatio-temporal models of the flexible dynamics are essen- tial. In this paper, a two-step approach for the identification of Measured output the spatio-temporal behavior of mechanical systems is developed and applied to a lightweight prototype industrial wafer stage. The proposed approach exploits a modal modeling framework and combines recently developed powerful linear time invari- ant (LTI) identification tools with a spline-based mode-shape Fixed world Flexible stage interpolation approach to estimate the spatial system behavior. The experimental results for the wafer stage application confirm the suitability of the proposed approach for the identification of complex position-dependent mechanical systems, and its potential Fig. 1. Schematic flexible wafer-stage system. for motion control performance improvements. Index Terms—system-identification, precision mechatronics include over-actuation and over-sensing [6], spatial vibration control [7], multivariable robust control [1], and inferential I. I NTRODUCTION control of unmeasured performance variables [8]. Invariably, Increasingly stringent performance requirements for pre- these approaches are characterized by an increased reliance cision motion systems lead to a situation where the flex- on model-based control design procedures, necessitating the ible dynamics of moving machine components need to be development of control-relevant, efficient, and numerically actively modeled and controlled. Typical examples include reliable identification algorithms capable of dealing with the the wafer stages in lithographic wafer scanners [1], [2]. complex spatio-temporal system behavior [2], [9], [10]. Traditionally, these stages can be accurately approximated The flexible dynamics of these systems in conjunction with as a rigid body in the frequency range relevant for control the fact that these motion systems move lead to position- [3], [4], thereby enabling static decoupling of the rigid-body dependent system behavior [4], [11]. In this paper, mechanical dynamics and subsequent decentralized control design [5]. systems consisting of a single flexible moving body are Furthermore, when this rigid-body approximation is used, considered and deformations are assumed to be small. As an the spatial system behavior follows directly from the stage example, consider the schematic flexible wafer-stage system geometry. Due to increasing accuracy and speed requirements, in Figure 1. Here, the flexible wafer-stage moves in relation the flexible dynamics in future systems can no longer be to the sensors, which are connected to the fixed world. As a neglected and need to be explicitly addressed. Approaches to result of this relative motion, the sensors measure the position address the resulting complex spatio-temporal system behavior at different points on the flexible structure, and therefore the spatial behavior of this flexible system is observed differently This research is supported by ASML Research as part of the TU/e Impulse I depending on the relative location of the sensors. However, program and is part of the research program VIDI with project number 15698, due to the fact that there is only a single moving body and financed by the Netherlands Organization for Scientific Research (NWO). Robbert Voorhoeve, Robin de Rozario, and Tom Oomen are with the because deformations are small, the structural dynamics of the Eindhoven University of Technology, Department of Mechanical Engineering, flexible body do not change as the system moves. Still, as a Control Systems Technology group, Eindhoven, The Netherlands (e-mail result of the position dependency in the measurement system, addresses: r.j.voorhoeve@gmail.com, r.d.rozario@tue.nl, t.a.e.oomen@tue.nl). Wouter Aangenent is with ASML Research Mechatronics, Veldhoven, The the system dynamics as observed by the sensors are no longer Netherlands (e-mail address: wouter.aangenent@asml.com). Linear-Time-Invariant (LTI), necessitating a deviation from the This article has been accepted for publication in IEEE Transactions on standard LTI control design techniques used in the context of Control Systems Technology (DOI: 10.1109/TCST.2020.2974140). 2020 IEEE. Personal use of this material is permitted. Permission from high-precision motion systems. IEEE must be obtained for all other uses, in any current or future media, Even though standard linear control design approaches including reprinting/republishing this material for advertising or promotional are no longer sufficient to control the position dependent purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. dynamics, the particular properties of the considered class of arXiv:1807.06942v2 [cs.SY] 28 Feb 2020 2 systems can be exploited to utilize control approaches that remain relatively close to LTI control theory. Approaches in literature that utilize additional system information to retain some aspects of LTI control theory include the gain-scheduling control approach [12]–[14] and the Linear Parameter-Varying (LPV) control framework [15]–[18], which formalizes the gain-scheduling method by ensuring stability and performance through a rigorous model-based mathematical approach. While the class of systems considered in this paper can be readily treated within the LPV control framework, the present paper proposes an approach that explicitly utilizes the additional assumptions of a single moving body to simplify the problem. A key challenge for systematic LPV control is the avail- Fig. 2. The experimental wafer-stage setup. ability of accurate LPV models. The need for accurate LPV models spurred the development of LPV system identification with a strong focus on black-box parametric models [19]–[22]. complexity of such a state-of-the-art industrial system. The This resulted in a well-developed theoretical framework that novelty of this paper therefore lies in the formulation and val- categorizes the identification techniques in local [23], [24], idation of the full position-dependent identification framework, and global approaches, e.g., [25], [26]. Depending on the which, while making use of previous results, including but not application, both approaches have been reported to effectively limited to previous results by the authors, e.g., in [10], [31], support the identification of practically relevant systems [27], has not been previously published. albeit that the identification of systems with high dynamic The identification methods used in this paper have parallels order remains challenging [21], [28]. Due to the high model with the field of experimental modal analysis. Research in complexity associated with general LPV modeling, the success this field has seen significant developments in the past decade of black-box LPV approaches is limited for the identification [32]–[34], in part due to recent consolidation efforts of the of mechanical systems with a large number of resonant modes approaches used in modal analysis and system identification [29], [30], showing a need to reduce the modeling complexity [35]–[37]. Contrary to modal analysis, the modeling goal in by using additional prior system knowledge. this paper is to obtain position-dependent models for control. Although many important developments have progressed This is reflected by the emphasis on accurate mode-shape LPV identification for control, the continuously increasing interpolation and the possibility to incorporate control-relevant complexity of motion systems necessitates a practical iden- identification criteria as in, e.g., [2]. tification approach of reduced complexity that is systematic, The outline of this paper is as follows. In Section II, the accurate, and user-friendly. A key step in this paper is to experimental setup is introduced and the control challenges as utilize the knowledge that the system consists of a single well as the position dependent modeling problem for this setup moving body with small deformations to derive a parsimonious are formulated. In Section III, the proposed two-step position- model-set, significantly reducing the modeling complexity dependent identification approach is explained. In Section IV, compared to a full LPV approach. Furthermore, recently the LTI identification approach and the obtained results for developed efficient and reliable LTI identification tools are the prototype wafer stage system are presented. In Section employed to obtain accurate and coherent local models which V, the approach and results for the mode shape interpolation are particularly suited to a subsequent interpolation step to are presented. In Section VI, a discussion is presented on the obtain the desired position-dependent model. The aim of this applicability of the proposed approach for control. In Section paper is to develop an effective and practical approach for the VII the conclusions of this paper are formulated as well as an identification of position-dependent mechanical systems. outlook on ongoing research. The main contributions of this paper are the following. 1) A two-step modal identification framework for position- II. E XPERIM ENTAL S ETUP AND PROBLEM FORMULATION dependent mechanical systems, including In this section the experimental wafer-stage setup and a) a flexible framework of parameterizations and al- related control challenges are outlined. In Section II-C, the gorithms aimed at obtaining accurate modal LTI position-dependent modeling problem is formulated, as is models of complex mechanical systems, considered in this paper. b) an approach for the interpolation of identified mode-shapes to obtain position-dependent models for control of flexible mechanical systems. A. Experimental Setup 2) Application and validation of the developed approach on The experimental setup considered in this paper is the Over- a state-of-the-art industrial wafer stage setup. Actuated-Testrig (OAT), as is shown in Figure 2. The system These contributions are inextricably linked as the proposed is considered here as an analogue to the flexible wafer-stage as two-step position-dependent identification framework is ex- depicted in Figure 1. The system is controlled in six motion plicitly developed with the goal of being able to handle the degrees of freedom and is magnetically levitated having no 3 w z P () c c K() Fig. 4. The LPV standard plant framework. See, for example, [8], [31], [38], for control design approaches for such problems. Fig. 3. Actuator and sensor configuration of the considered experimental setup. The locations of actuators used for control and identification are marked The LPV standard plant framework, as depicted in Figure with red crosses ( ) actuators used only for identification as blue crosses ( ), 4, can be used to describe both these control problems. sensors used for identification and control as red circles ( ), and sensors used Here, u and y are the output and input signals available only for identification as green circles ( ). c c for control, w is the generalized disturbance signal and z p p is the generalized performance signal, which in this case involves the positioning error of the point of interest. This mechanical connection to the fixed world. Furthermore, the control problem, with the LPV standard plant P () in a experimental system is equipped with additional actuators generalized feedback interconnection with an LPV controller and sensors to facilitate the spatio-temporal identification of K (), has been studied extensively in, e.g., [15]–[17] and, for the flexible dynamics. While these additional actuators and appropriately bounded sets of scheduling variable trajectories, sensors are available on the experimental setup, they are not i.e., (t) 2 D , efficient algorithms exist for various robust and available in the considered wafer-stage system as depicted in optimal control problems defined in this framework. Figure 1 for which the experimental system is an analogue. The main difficulty for the practical application of these Therefore, these additional actuators and sensors are only used methods concerns the availability of an accurate LPV system- for identification and are considered to be unavailable for the model G(), which is the part of the standard plant P () control of the system. In this paper, only the out-of-plane pertaining to the physical system that is to be controlled. That motions are considered, i.e., the motions perpendicular to the this is a difficult problem is evidenced by the fact that accurate surface of the wafer-stage, both for visualization purposes and modeling of LTI precision systems is already considered to be due to the availability of multiple spatially distributed sensors a challenging problem, see, e.g., [9], [10]. Accurate modeling in this direction. for LPV systems is generally significantly more challenging, The out-of-plane sensors and actuators used for the exper- since the incorporation of a scheduling parameter dependency iments in this paper are shown in Figure 3. Seven actuators, typically leads to a highly increased model complexity [21]. depicted by crosses in Figure 3, four of which, shown in Due to the additional constraints on the considered class of red, are used for closed-loop control and the remaining three, systems, consisting of one flexible moving body and assum- shown in blue, are used to apply additional spatially distributed ing small deformations, the complexity of the problem can excitation for identification. Sixteen sensors are available for be significantly reduces. This problem of obtaining accurate identification, as shown in Figure 3 by circles. Similarly, a position-dependent models for such mechanical systems, such distinction is made between the sensors used for closed-loop as the wafer-stage example considered here, is addressed in control, shown in red, and those only used for the spatially the remainder of this paper. distributed identification, which are shown in green. C. Position-Dependent Modeling Problem B. Control Challenges Consider the schematic wafer-stage setup shown in Figure In the considered wafer-stage example, the scheduling vari- 1. For this setup, two key control challenges are recognized ables  relate to the position of the wafer stage. Due to the which are related to the position-dependent system behavior. fact that there is only a single moving body, the structural First, as previously outlined, the relative motion of the wafer dynamics of the flexible body are invariant to the position stage with respect to the sensors leads to an input-output be- of the wafer stage with respect to the sensors. Furthermore, havior that is no longer LTI, necessitating a position/dependent deformations are considered to be small in the sense that the modeling and control perspective. Second, the point of interest, structure does not deform in such a way that it influences i.e., the point on the wafer that is being exposed in the photo- the dynamics of the systems. As a result, the A matrix of lithographic process, also changes as the wafer stage moves. the state-space description for this class of systems does not This involves the control of an unmeasured performance vari- depend on the scheduling variable , see, e.g., [39], [40]. able since the point of interest is not directly measured. This Furthermore, for the considered system the positions of the results in a position-dependent inferential control problem. actuators are fixed relative to the wafer-table and therefore the 4 input matrix B is also not position dependent. Therefore, the that is applicable for any geometrically complex domain D is scheduling variables have no influence on the system states the Finite Element Method (FEM). This approach uses many and as a result there is no memory in the system pertaining to localized basis functions to accurately approximate the spatial the past-trajectory of the scheduling variables. The system is system behavior [41]. therefore only dependent on the current, instantaneous, value The temporal system behavior for this basis function ap- of the scheduling parameters. In state-space the considered proach is determined by the dynamics of the generalized system-model is described as, coordinates, q(t) = [q (t) : : : q (t)] . Under the assumption 1 n of small rotations and strains, and assuming that the material is linear elastic obeying Hooke’s Law, the equations of motion x _ (t) = Ax(t) + Bu(t) (1) that govern the temporal input-output behavior of a mechanical G() : y(t) = C ((t))x(t) + D ((t))u(t) (2) y y system are given by the set of coupled second order ordinary z(t) = C ((t))x(t) + D ((t))u(t): (3) z z differential equations [39, Section 2.2] In the wafer-stage example, both the system outputs y(t) Mq (t) +Dq_(t) +Kq(t) = Qu(t) ; (7) and the performance variables z(t), i.e., the point of interest n n n n position, can be considered as specific local instances of the q q q q where M 2 R is the mass matrix, D 2 R the n n out-of-plane deflection of the surface of the wafer stage. In this q q damping matrix, K 2 R the stiffness matrix, and Q 2 interpretation, the relevant output for this system is this out-of- n n q u R the input distribution matrix. plane deflection, which is denoted here as z(%; t), where % is The set of coupled equations of motion (7) can be decoupled the in-plane coordinate of the point for which the deflection is for the undamped case by transforming to a modal descrip- considered. The modeling problem is then to identify from tion, which is obtained by solving the generalized eigenvalue experimental data of the system, a model for the system problem behavior of the form K ! M  = 0 ; i = 1; : : : ; n : (8) i q x _ (t) = Ax(t) + Bu(t) (4) i G(%) : z(%; t) = C(%)x(t) + D(%)u(t) : (5) 2 The eigenvalues, ! , are the squared undamped resonance- frequencies of the modes, and the eigenvectors,  , are the The model in the form (1)–(3) is recovered from this descrip- i associated mode shapes as parameterized in the basis W (%) = tion by including the static geometric relations between the position of the wafer stage  and the coordinates %, at which [w (%) : : : w (%)]. By applying the substitution q = , 1 n 1 1 the sensors view the wafer-stage as well as the coordinate where  = [ : : :  ], and multiplying (7) with  M 1 n of the point of interest. This is a simple affine relation yields dependent on the definitions of the origins for two coordinate I (t) +D  _ (t) + (t) = Ru(t) ; (9) systems and the sensor locations, i.e., for a certain sensor y , G (%) : z(%; t) = L(%)(t) ; (10) C ((t)) = C(% ((t)) with % ((t)) = % + (t). The y y y y ;0 1 1 1 1 modeling problem considered in the remainder of this paper 1 1 2 2 2 where D =  M D, = diag([! : : : ! ]), 1 n is therefore to model the system of the form (4)–(5). 1 1 L(%) = W (%), and R =  M Q. In the context of identification for control, low-order models III. POSITION-D EPENDENT M ODELING APPROACH are desired that are accurate in a limited frequency band of In this section the proposed position-dependent modeling interest [2]. This means that only a limited number of modes, approach is described. First, the modal modeling framework n < n , are required to model the relevant temporal system m q for mechanical systems is outlined. Next, the proposed two- behavior [39]. Modeling the spatial system behavior, using step identification approach is developed. a generic set of basis functions L(%) = W (%), typically requires a large number of basis functions leading to a high modeling complexity. In this paper, a two-step identification A. Modal Models of Mechanical Systems approach is proposed to directly identify the mode-shapes, i.e., The quantity of interest is the out-of-plane deflection, the columns of L(%) from measured data. z(%; t) : D T 7! R, of the surface of a flexible body, which is modeled here as a continuum. The domain D 2 R of the coordinate % is the surface of the considered structure, e.g., B. Two-Step Identification Approach the (x; y) surface of the wafer-stage. To model the spatio- To identify the spatio-temporal system behavior, measure- temporal evolution of z(%; t), a basis-function expansion is ment data is first obtained in experiments with fixed sensor used for time–space separation, see, e.g., [41], i.e., locations % . As a result of the fixed sensor locations the input-output system behavior is linear time invariant, similar z(%; t) = w (%)q (t) : (6) i i to the local approach in LPV identification. Experiments are i=1 performed with various fixed sensor locations covering the For n ! 1 this expansion converges as long as fw (%)g domain D. An LTI system model is then identified from the q i i=1 is a convergent set of functions for the class of continuous obtained experimental data where the model is parameterized functions on the spatial domain D [41]. A widely used method in modal form, i.e., (9)–(10), with n modes. Instead of m 5 nc v the position-dependent output equation (10) the measured, spatially sampled, outputs z (t) are modeled as s u G 2 3 2 3 z(% ; t) L(% ) 1 1 6 7 6 7 . . . . z (t) = = L (t); L  (11) s 4 5 s s 4 5 . . Fig. 5. Closed-loop identification scheme. z(% ; t) L(% ) n n % % n n % m where the parameters in L 2 R are considered as spatially sampled estimates of the mode shapes L(%). identification scheme is shown in Figure 5. A distinction is This first step requires the LTI identification of a complex made between inputs which are used in the control loop u (t) mechanical system with a high model order and many inputs and those that are not used in the control loop u (t). The nc and outputs. The identification of such complex mechanical control inputs u (t) also include the in-plane actuator signals systems requires the use of efficient and numerically reliable that are used to stabilize the in-plane rigid-body modes. The identification approaches, as have been developed and inves- excitation signals used for system identification are the non- tigated in, e.g., [9], [10]. control inputs u (t) and the additive perturbation signals nc In the second step, the spatial mode shapes L(%) are r (t) as shown in Figure 5. The applied excitation signals estimated from the identified parameters in L . In this step, are all random-phase multisines with a flat amplitude spectrum interpolation techniques are used to reconstruct continuous which are successively applied to the single inputs in separate mode shapes based on these spatially sampled estimates. Since experiments. this step involves the interpolation of spatial functions in % and With a total of 16 out-of-plane sensors, 8 control inputs, not of systems that dynamically depended on a scheduling including 4 in the in-plane direction, and 3 non-control inputs, variable , the interpolation pitfalls as shown in [42] are the identification problem involves first identifying a 24 11 avoided. In Section V, a promising robust and physically closed-loop FRF given by motivated interpolation approach is proposed, but several other " # interpolation techniques, which might be more suitable for ~ ~ P ( ) P ( z r k z u k s u s nc ~ c other applications, can be used in this second step of the P ( ) = ; (12) CL k ~ ~ P ( ) P ( u r k u u k c u c nc proposed two-step approach. c In summary, the proposed two-step approach aims to: where the notation P ( ) is used to denote the identified y x k 1) Identify the modal mechanical LTI model given by (9) empirical transfer function estimate (ETFE) from the input and (11), i.e., estimate the parameters in L ; ; D ; s m signal x to output signal y at frequency point . To obtain R. the FRF of the open-loop system G from this closed-loop FRF 2) Estimate the mode-shapes L(%) based on the spatially the following relation is used, sampled mode-shapes L " # h i ~ ~ P P u r u u c u c nc ~ c IV. LTI I DENTIFICATION OF SPATIALLY SAMPLED ~ ~ G( ) = P P ; z r z u s u s nc ~ ~ P P u r u u S YSTEMS nc u nc nc (13) In this section, the first step of the proposed identification where the arguments have been omitted on the right hand approach is outlined, which is the LTI identification of the side of this expression for brevity. In (13) P = 0 and u r nc u spatially sampled system G . P = I , see, e.g., [2, Appendix A] for additional detail u u nc nc on this closed-loop identification approach. By removing the A. Methods in-plane input directions, the 16  7 FRF of the system G , as given by (9) and (11), is obtained. The LTI identification approach considered here involves a number of key aspects. First, a non-parametric identification In this paper, the delays from the hold circuit in the digital approach is considered, aimed at obtaining accurate Frequency measurement environment are first carefully determined using Response Function (FRF) estimates of the spatially sampled a combination of readily available automatic methods and system G . Second, the modal parametrization as used in this visual inspection of the data while manually tweaking the paper is defined. Third, a black-box Matrix Fraction Descrip- delay value. Subsequently the the FRF measurements are tion (MFD) parametrization is employed, which is parameter- compensated for these delays such that the obtained delay- ized such that it closely matches the modal parametrization. compensated FRF can be modeled in continuous time, i.e., Fourth, the identification algorithms used to estimate the by G (s ), with s = j! , and the remaining identification s k k k system models from the measured data are explained. procedure can be performed in the s-domain. For additional 1) Non-Parametric Identification: The non-parametric fre- details on such a pseudo continuous time modeling approach, quency response function for the wafer-stage system is esti- see, e.g., [44], [43, Section 8.5]. Note that when the model is to mated using the robust multisine approach as explained in, e.g., be used in a control context, the importance of determining the [43, Section 3.7]. The rigid-body motions of the system need correct delay increases significantly and one might consider to be controlled for stable operation, therefore all experiments changing to a fully discrete-time z-domain approach which is are performed in a closed-loop configuration. The closed-loop generally better suited to deal with such unknown delays. 6 2) Modal Parametrization: As outlined in Section III, the identification of increasingly complex system where the the modal model for the spatially sampled system G is numerical conditioning becomes an important limiting factor given by (9) and (11), with parameters contained in the for the performance of the identification algorithms. matrices L ; ; D ; and R. The matrices L and R are s m s Furthermore, this general parametrization enables the use 2 2 fully parameterized while = diag(!  ) with parameters of more structured LMFD parameterizations, which enable 2 2 2 T !  = [! : : : ! ] . The damping matrix D is either fully m the parametrization of system with arbitrary McMillan degree 1 n parameterized, which constitutes a general viscous damping n instead of only being able to parameterize systems where model, or, when using the modal damping model, this is equal the degree n is a multiple of the number of outputs p, as to D = diag( ) with  = [ : : :  ] . m;mod 1 n is the case when using the fully parameterized unstructured Using the modal damping model, the set of differential LMFD as in, e.g., [32], [33]. Here, a generic (pseudo)- equations (9) describing the systems temporal behavior be- observable-canonical LMFD parametrization with a quasi- comes fully decoupled, meaning the system can be considered constant degree structure is used, see, e.g., [34], [45], [46]. as a superposition of independently evolving modes [39]. This parametrization is both identifiable, in the sense that it This representation is extensively used in modal analysis and is not over-parameterized, and generic, meaning that it can design as it simplifies the physical interpretation of the modal approximate all proper LTI systems of the given order up to parameters and the incurred modeling errors by assuming arbitrary precision, see [45]. modal damping is generally small for lightly damped sys- This LMFD parametrization is often used for black-box tems, see [39, Section 2.4]. For the wafer-stage example the identification of LTI systems, whereas in this paper the goal modal damping assumption is used to facilitate a parsimonious is to identify spatio-temporal mechanical systems by utilizing parametrization. Note that the modeling framework used in the modal form, i.e., (9) and (11). Therefore, in this paper this paper enables general linear damping models. a number of additional constraints are incorporated in the The complexity of this parametrization is determined by LMFD parametrization such that it more closely resembles the number of modes that are modeled, n . While a limited the mechanical system model. Here, the following properties number of modes usually dominate the behavior in a given are enforced, frequency range of interest, for some applications the com- 1) an even McMillan degree, by taking n = 2 n , x m bined low-frequency compliance effect of unmodeled higher 2) a relative degree r  2, by constraining the column order modes also needs to be taken into account. In such a degrees of the numerator polynomial matrix N (s; ) to case it is relevant to model an additional compliance term, be 2 lower than the corresponding column degrees of e.g., as direct feed-through such as D(%) in (5), to describe the denominator polynomial matrix D(s; ), the quasi-static deformation resulting from each input signal 3) a prescribed number of rigid-body modes n such that rb u(t), see, e.g., [38]. In this paper such a feed-through term n = 2 n poles are located at s = 0, by factoring out 0 rb is not modeled; this can be straightforwardly incorporated in the rigid body dynamics, see [10] for details. the proposed modeling framework when required for a given application. 4) Identification Algorithms: The identification problem The modal parametrization used in this paper is defined by considered here is to find the parameter vector  that min- (9) and (11) with the modal parameters given by, imizes the identification cost, which is in this paper is a h i weighted least-squares cost function, i.e., T 2 = vec : (14) m L R ! 3) MFD Parametrization: For certain identification algo- = arg min V () = "(s ; ) "(s ; ) ; (17) k k rithms, it is necessary that the model parametrization can be k=1 written as a polynomial matrix fraction description. In this pa- where per, a Left Matrix Fraction Description (LMFD) parametriza- tion is used, given by ~ ^ "(s ; ) = W (k) vec G(s ) G(s ; ) ; (18) k k k ^ ^ ^ G(s; ) = D(s; ) N (s; ) ; (15) pqpq with weighting matrix W (k) 2 C . This general cost pq pp ^ ^ where N (s; ) 2 R [s] and D(s; ) 2 R [s] are real function V () includes other commonly used identification polynomial matrices in the Laplace variable s. Furthermore, criteria [10], such as the sample maximum likelihood criterion these polynomial matrices are linearly parameterized with [43], control relevant identification criteria [2], and the input- respect to the parameter vector  2 R using a set of basis output and element-wise weighted criterion used in [47]. functions such that, Minimization of the cost function (17) is a nonlinear least- h i squares optimization problem. Suitable algorithms to solve this ^ ^ vec = (s) = (s) : (16) D(s; ) N (s; ) j j problem are the Sanathanan-Koerner algorithm and the Gauss- j=1 Newton algorithm, or closely related methods such as the This general linear parametrization allows the use of data- Levenberg-Marquart algorithm. These algorithms are defined dependent orthogonal vector polynomials as basis functions, as follows, where the Sanathanan-Koerner algorithm is only p(q+p)1 2 R [s], see, e.g., [9], which is a key aspect for defined for MFD parameterizations, i.e., using (15)–(16). j 7 h0i Algorithm 1 (Sanathanan-Koerner [48]): Let  be given. LM FD: SK LM FD: LM In iteration i = 0; 1; : : : , solve the linear least squares problem Modal: initial estimate Modal: LM hi+1i hii = arg min W (s ;  ) (s ) ; (19) SK k k k=1 with h i hii hii 1 T ^ W (k;  ) = W (k) D(s ;  ) : G (s ) I sk k k q (20) Algorithm 2 (Gauss-Newton [49]): Given an initial estimate 0 5 10 15 20 25 30 35 40 45 50 h0i , compute for i = 0; 1; : : : iteration number hi+1i hii hii hii Fig. 6. Evolution of the cost function during various steps of the LTI =  + arg min J (s ;  ) + "(s ;  ) ; k  k identification approach. k=1 (21) with are beyond the scope of this paper and will be reported @"(s ; ) @ vec(G(s ; )) k k hii elsewhere. Due to the modal damping assumption used in this J (s ;  ) = = W (k) : T T @ @ hii paper for the modal model, and since this modal damping hii (22) assumption is not enforced in the MFD parametrization, the The Gauss-Newton algorithm, and related algorithms such transformation in step 4 of this identification approach is as the Levenberg-Marquardt algorithm, generally provide fast approximate. This approximate transformation is performed monotonic convergence to a minimum of the cost function by calculating the n = 2n pole locations by solving x m V (). However, these algorithms often converge to local- det[D(s; )] = 0, separating the poles into pole pairs such 2 2 minima that are far from optimal, and therefore their per- that (s + p )(s + p ) = s +  s + ! , with  ; ! 2 R for 1;i 2;i i i i formance is strongly dependent on the quality of the ini- i = 1; : : : ; n , and estimating a model of the form, h0i tial estimate  . The Sanathanan-Koerner algorithm on the other hand does not generally converge monotonically, and, ^ G = ; (23) m;trans s +  s + ! if convergent, its stationary points are generally not optima i i=1 of the cost function [50]. However, the Sanathanan-Koerner pq with R 2 R . In this estimation the denominator param- algorithm often yields adequate, albeit suboptimal, results eters are fixed to the values obtained from the pole-pairs of irrespective of the quality of the initial estimate. Therefore, the the LMFD model. This model is fitted based on the FRF data Sanathanan-Koerner algorithm is often used to provide initial using the same cost function as the other identification steps, estimates that are subsequently refined using a gradient based i.e, using (17). The parameters in L and R are then obtained optimization algorithm such as the Gauss-Newton algorithm, from the singular value decomposition of the residue matrices see, e.g., [10]. R = U S V such that i i i 5) Identification Approach: To identify the modal system i 1 1 H model defined by (9), (11) and parameterized by (14) from [L ] = [U ] [S ] ; [R] = [V ] ; i = 1; : : : ; n ; (24) s i i i 1 m 1 i the identified FRF G (s ), the following steps are followed. s k where [X ] and [X ] denote, respectively, the i-th and column 1) Define weighting filters W (k) as in (18). and the j-th row of matrix X . This transformation performs 2) Perform i iterations of the SK algorithm using the SK well for the considered system, as shown by the results in mechanical LMFD model with constraints as defined in Section IV-B. For more general approaches to transform black- Section IV-A3. box models to a gray-box models, see, e.g., [51]. 3) Perform i iterations of the GN or LM algorithm GN for the LMFD model using the parameters SK;min B. Results corresponding to the lowest cost function value during For the identification of the wafer-stage system, the weight- SK iterations as initial estimate. ing function in (18) is chosen as a weighting with the element- 4) Transform the identified LMFD model to an initial wise inverse of the identified FRF G(s ), i.e., estimate for the modal model as defined by parameters k (14). W (k) = diag(vec(j(s )j)) ; (25) inv k 5) Perform a maximum of i iterations of the GN GN;mod or LM algorithm to converge to a optimum of the cost [(s )] = : (26) [G(s )] function as in (17) for the identified modal model. j When considering generally damped modal models, the This choice reflects the goal of minimizing the relative error ^ ~ fourth step of this identification approach can be performed between the model G(s ; ) and the FRF G(s ). For more k k using an exact transformation, i.e., relating two realizations advanced weighting choices that take into account the control of the same system. Details of this exact transformation objective, see [2]. To emphasize the accurate estimation of the V 8 -100 -150 -200 -100 -150 -200 -100 -150 -200 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 Frequency [Hz] Fig. 7. Bode diagrams of the identified FRFs (dotted blue) and fitted modal model (solid red) for a 3 7 part of the full 16 7 identified system. first few lower-frequency resonances the weighting function is (14), which leads to a slight increase in the cost function truncated, i.e., value. In the final step the cost function is again minimized using the monotonically convergent LM algorithm with the W (k) = min (W (k); w ) ; (27) inv max modal parametrization, which in this case only yields a small decrease of the cost function value showing that the initial where w is chosen such that clipping of the weight max modal estimate is already of a high quality. generally occurs only in the high frequency range, after the The small increase in cost function value when transforming first few resonances. the LMFD modal to the modal model is expected, as this step Steps 2-5 of the identification approach as proposed in reduces the model complexity by enforcing modal damping Section IV-A5 are at first only performed for a 3  7 part and through the elimination of computational modes identified of the full 16 7 identified FRF. This is done both to improve in the LMFD model. This is done by visually comparing the computational efficiency and to simplify the implementation identified pole-locations with the resonances of the identified of rigid-body mode constraints, see, e.g., [10]. After the FRFs, similar to the use of stabilization diagrams in experi- successful identification of the modal parameters for the 3 7 mental modal analysis, see, e.g., [32], [33]. The model order system, the model for the remaining 13 outputs is determined of the identified LMFD model is n = 2n = 46 while the x m by estimating the sampled mode-shape parameters in L for model order of the modal model is n = 24. This shows these outputs while all other parameters remain fixed. This that the transformation yields a significant decrease in model is again done by minimizing the cost function (17) for these complexity with a modest increase in cost function value. additional output, which in this case is simply a linear least squares problem. Figure 7 shows the FRF and identified modal model for the In Figure 6, the evolution of the cost function, V () in 37 part of the system on which the identification procedure is (17), is shown during steps 2-5 of the identification approach performed and Figure 8 shows a more detailed view including of Section IV-A5. As can be seen from this figure, the SK the phase of a single element of the FRF and identified modal algorithm in step 2 is not monotonically convergent, but model. These figures show a good agreement between the yields an appropriate initial estimate for subsequent refinement model and the FRF in the low frequency region as well as using the LM algorithm. The LM algorithms used in step for the dominant resonances in the frequency region up to 3 does show monotonic convergence and yields an LMFD about 1 kHz. In the frequency region beyond 1 kHz, there estimate with a cost function value approximately one order are an additional few accurately identified modes but also a of magnitude below that of the initial estimate provided by the number of unmodeled resonances. The identification of these SK approach. In the next step the LMFD model is transformed high frequency modes is not the main focus in this paper since to the modally damped model defined by the parameters the spatial behavior for such high frequency modes is typically Magnitude [dB] 9 such that U , which is generally interpreted as a measure of the bending energy of the functions, exists for all functions in -100 the space. The functions W that minimize (28) are given by -150 W (x; y; #) = # + x# + y# + # G (x; y); (30) s 0 x y j j -200 j=1 180 2 2 2 G (x; y) = r ln(r ); r = (x  x) + (y  y) ; (31) j j j j j see, e.g., [52]. The number of parameters in (30) is n + 3, where the three additional parameters are related to the -90 monomials up to the first degree which represent the set of functions in W for which U = 0, i.e., the kernel of U . To -180 constrain this underdetermined set of equations, the following 1 2 3 10 10 10 three additional constraints are added which make sure that Frequency [Hz] the function space parameterized using the Green’s functions G (x; y) is orthogonal to the space of first order polynomial, Fig. 8. Detailed Bode diagram including phase of the identified FRF (dotted blue) and fitted modal model (solid red) for the (1; 2) element of the identified n n n % % % X X X system as depicted in Figure 7. # = 0 ; # x  = 0 ; # y  = 0 : (32) j j j j j j=1 j=1 j=1 The solution to (28) using (30) and (32) is given by also of a higher spatial frequency, meaning a higher spatial resolution, than what is available, is required to identify the z  = W (x  ; y  ; #) + # ; (33) k s k k k associated mode-shapes. see, e.g., [52]. In explicit matrix-form this yields h i V. M ODE-SHAPE INTERPOLATION T # = X ; (34) z  : : : z  0 1 n 13 In this section, the interpolation of the spatially sampled h i mode shapes, as given by the columns of L , is considered. T with # = , and where # # # # : : : # 0 x y 1 n This interpolation step is performed to obtain a position- % dependent model that is continuous in the spatial variable %. 2 3 " # 1 x  y 1 1 X X + I 6 7 0 G . X = ; X = ; (35) 0 4 5 A. Methods 0 X 33 A popular method for the interpolation of various types of 1 x  y n n % % 2 3 data at arbitrary spatially distributed points is the smoothed G (x  ; y  ) : : : G (x  ; y  ) 1 1 1 n 1 1 thin-plate-spline interpolation approach. The use of thin-plate- 6 . . 7 . . X = : (36) 4 5 . . splines is physically motivated by the fact that the spline functions are derived as the functions that minimize the G (x  ; y  ) : : : G (x  ; y  ) 1 n n n n n % % % % % bending energy of a thin sheet of elastic material. Therefore, For the interpolation of the spatially-sampled mode shapes this approach is particularly well-suited for the considered as identified in L , the points (x  ; y  ) are equal to % , i.e., the s j j j application of interpolating the structural mode-shapes of (x; y) positions of the sensors, and the values for z  are given motion systems which are thin in one dimension relative to by the identified parameters in the columns of L . This inter- the other dimensions, such as the wafer-stage example. polation is carried out independently for each mode shape, i.e., The smoothed thin-plate-spline interpolating function W for each column of L , where for each mode shape a different for a single mode-shape is derived as follows. Given a set of smoothing parameters  is used. These smoothing parameters n points f(x  ; y  ; z  ) 2 R g and a user-defined smoothing % j j j provide a trade-off between robustness to estimation errors in parameter  2 [0;1), find an interpolating functionW 2 W L and interpolation accuracy at the data-points % . In this s j such that, paper, the values of the smoothing parameters are determined using a Leave-One-Out-Cross-Validation (LOOCV) approach, min jW (x  ; y  ) z  j + U; (28) s j j j i.e., the value of  is used which minimizes the LOOCV error. W 2W j=1 For details on this cross-validation approach, see, e.g., [53]. with 1 1 Z Z B. Results U =  W (x; y) dx dy ; (29) In Figure 9, four of the identified flexible mode-shapes are 11 shown. In total 12 mode-shapes are identified including the Here, the function space W is the space of continuously dif- three out-of-plane rigid-body modes. Up to the ninth mode, as ferentiable functions with square-integrable second derivatives shown in Figure 9d, the identified mode shapes agree well with Phase[deg.] Magnitude [dB] 10 Eigenfrequency: 139.23 Hz Eigenfrequency: 437.46 Hz 1.5 1.5 1 1 0.5 0.5 0 0 -0.5 -0.5 -1 -1 -1.5 -1.5 0.5 0.5 0.5 0.5 x x y y -0.5 -0.5 -0.5 -0.5 x x (b) Fifth mode (saddle). (a) Fourth mode (torsion). Eigenfrequency: 506.08 Hz Eigenfrequency: 879.46 Hz 0.6 1.2 0.4 1 0.2 0.8 0 0.6 -0.2 0.4 -0.4 0.2 -0.6 0 -0.8 -0.2 -1 -0.4 -1.2 -0.6 0.5 0.5 0.5 0.5 x x 0 0 0 0 y y -0.5 -0.5 -0.5 -0.5 x x (c) Sixth mode (umbrella). (d) Ninth mode. Fig. 9. Top views and 3-D surface plots of the identified mode shapes with red dots indicating the points of the identified spatially-sampled mode shapes. theoretical mode-shapes for a thin-plate or the mode-shapes By utilizing the identified position-dependent model of the as obtained by means of a Finite-Element-Method analysis of wafer-stage system, G (%), such global and local control the system, a detailed comparison is omitted for brevity. For approaches can be described in the LPV standard plant frame- the higher-order modes the spatial resolution of the sensors is work, as shown in Figure 4, and can be solved using a range of insufficient to accurately reconstruct the smooth mode-shapes. approaches. In this section, several global and local approaches The results in Figure 9 show the viability of the proposed are considered for both feedback and feedforward control. approach to obtain accurate position-dependent models of flexible mechanical systems. In the following section, the A. Global Spatio-Temporal Feedforward Control potential of the proposed approach is discussed in enabling For the wafer-stage example, the problem of global feedfor- various position-dependent control approaches. ward control aims to minimize the error between the reference signal r (t) and the out-of-plane deflection of the surface VI. OUTLOOK FOR POSITION-D EPENDENT M OTION of the wafer stage z(%; t) over the entire spatial domain D. CONTROL APPLICATIONS More precisely, the global approach is aimed at minimizing In this section, several control approaches are explored the following weighted spatial-norm of the error e(%; t) = that are enabled by the availability of accurate position- r (t) z(%; t) dependent models. For flexible motion systems both feedback and feedforward control problems become more complex as Z Z the point that should track the reference is often not directly kek = e (%; t)W (%)e(%; t) d% dt : (37) 2(D ) S measured, such as the point-of-interest of the wafer stage example in Figure 1. Furthermore, the location of this point- of-interest can change over time, increasing the complexity of This problem can be effectively formulated and solved as an the control problem. H optimal control problem, for details see [31]. Approaches to enhance the control performance for such In Figures 10 and 11, simulation results are shown for flexible motion systems can generally be classified as either, the wafer stage example system. Here, Figure 10 shows the 1) global approaches, aimed at preventing or mitigating the reference profile for r as well as the (x; y) coordinates of the flexible deformations in the entire system, such that the point of interest % (t) and Figure 11 shows the local errors deformation related errors are small; or at the point of interest, i.e., e (t) = r (t) z(% (t); t), for a z z z 2) local approaches, aimed at controlling the position of classical feedforward controller, that minimizes the error at the the point-of-interest of the deformed system sensor locations, and the global feedforward approach. These z z y y z z y y 11 % G() z r e u y x z z^ 0.8 % + zy 0.6 0.4 z^ O() 0.2 obs 0 0.005 0.01 0.015 0.02 0.025 0.03 Fig. 12. Position-dependent observer-based feedback control scheme. t [s] Fig. 10. Reference profile r , and the (x; y)-coordinates of the point of interest % over time. In [38], an observer-based inferential control approach is proposed which is especially suited to minimize the influence of disturbance-induced compliant deformations, i.e., the quasi- 0.05 static deformations induced by a locally applied force, often modeled using an additional feedthrough term, see, e.g., [35]. Combined with a moving disturbance source and point-of- -0.05 interest location, the position dependency of this compliant Classic control effect necessitates a position-dependent control approach to -0.1 Global control obtain the desired performance. This position-dependent con- 0 0.005 0.01 0.015 0.02 0.025 0.03 troller can be effectively and intuitively realized by combining t [s] an observer containing a position-dependent system model with an LTI controller that is robust to the remaining position- Fig. 11. Inferential error e at point of interest over time, showing signifi- dependency, as shown in Figure 12. cantly improved performance for the global feedforward controller relative to the classical feedforward controller. VII. C ONCLUSIONS AND OUTLOOK results show that the global approach yields a significantly A. Conclusions improved inferential performance as opposed to the classical This paper provides a general procedure for the identifi- approach. These results and the practical potential of this cation of position-dependent precision mechatronic systems global feedforward approach for future motion systems are which consist of a single flexible moving body and with small enabled by the modeling procedure developed in this paper. deformations. This is an essential step for the control of future high-precision motion systems. A key step in the proposed B. Local Inferential Feedback and Feedforward Control approach is to utilize prior mechanical systems knowledge The general problem of local inferential control aims to as embedded in the modal modeling framework to obtain a directly optimize the performance at the location where per- parsimonious model set. formance is required, such as the point-of-interest in the wafer- In Section IV, a flexible framework of parametrizations stage example, see Figure 1. Both for feedback and feedfor- and identification algorithms is proposed that is especially ward, these approaches can be described in the LPV stan- suited for the identification of modal models of mechanical dard plant framework when an accurately identified position- systems. For the considered state-of-the-art industrial wafer dependent model is available. Local inferential feedforward stage system, with a total of 7 considered inputs and 16 control typically involves the inversion of the system dynamics outputs, the proposed identification approach yields a very from the inputs to the time-varying or parameter-varying accurate modal system model with 12 identified modes. The performance variables which, apart from explicit inversion, spline-based interpolation approach proposed in Section V can be realized by solving an optimal control problem or provides a robust and effective method to reconstruct the using iterative learning control, see, e.g., [18], [54]. Inferential spatial mode shapes, and is successfully applied to reconstruct feedback control, generally requires the use of two degree of 9 of the identified mode shapes. freedom controller structures as opposed to the single degree Potential applications of the proposed position-dependent of freedom controller structure used in traditional feedback modeling approach for control are numerous, including, e.g., control, see, e.g., [8], [38], this can be directly incorporated the use in global spatio-temporal feedforward control and in the standard-plant approach. observer based inferential feedback control. Inferential feedback control is especially relevant when significant disturbance forces are present in the system. In B. Outlook [38], it is shown that a disturbance-observer can be effectively utilized to estimate the inferential performance variable in In this paper, systems are considered that can be written the presence of significant disturbance forces that are non- as (4)–(5). Although this description is less general than (1)– collocated with the actuator forces. In [38], it is also shown (3), it is envisaged that the proposed framework of identify- that in such a case it is essential to include disturbance models ing modal models of mechanical systems and subsequently in the standard-plant description to obtain accurate results. interpolating the spatial system behavior can be extended normalized signals Inferential error e z 12 to be more broadly applicable. In particular, extending the [18] C. Hoffmann and H. Werner, “A survey of linear parameter-varying con- trol applications validated by experiments or high-fidelity simulations,” proposed framework to consider the modeling of interact- IEEE Transactions on Control Systems Technology, vol. 23, no. 2, pp. ing mechanical subsystems provides the ability to model a 416–433, March 2015. variety of relevant mechanical systems. Examples of such [19] B. Bamieh and L. Giarr, “Identification of linear parameter varying models,” International Journal of Robust and Nonlinear Control, vol. 12, interacting mechanical subsystems include the much used H- no. 9, pp. 841–853, 2002. bridge systems as considered in, e.g., [29], [30]. By separately [20] M. Lovera and G. Mercere, ` “Identification for gain-scheduling: a bal- considering the spatio-temporal behavior of the subsystems, anced subspace approach,” in Proceedings of the 2007 American Control Conference, July 2007, pp. 858–863. such as the beam and carriage in an H-bridge system, and [21] J.-W. van Wingerden and M. Verhaegen, “Subspace identification of modeling the full system behavior as a general interconnection bilinear and LPV systems for open- and closed-loop data,” Automatica, of these component models, a more general class of systems vol. 45, no. 2, pp. 372–381, 2009. can indeed be described, including systems with position [22] R. Toth, ´ Modeling and identification of linear parameter-varying sys- tems, ser. Lecture Notes in Control and Information Sciences. Berlin, dependent state-matrices. Validating the practical applicability Germany: Springer-Verlag Heidelberg, 2010, vol. 403. and performance of this approach, as well as the control [23] J. De Caigny, J. F. Camino, and J. Swevers, “Interpolation-based mod- approaches discussed in Section VI, is a subject of ongoing eling of MIMO LPV systems,” IEEE Transactions on Control Systems Technology, vol. 19, no. 1, pp. 46–63, 2011. research. [24] D. Vizer, G. Mercere, O. Prot, E. Laroche, and M. Lovera, “Linear fractional LPV model identification from local experiments: an H - REFERENCES based optimization technique,” in Proceedings of the 52nd Conference on Decision and Control, Florence, Italy, 2013, pp. 4559–4564. [1] M. van de Wal, G. van Baars, F. Sperling, and O. Bosgra, “Multivariable [25] F. Felici, J. van Wingerden, and M. Verhaegen, “Subspace identification H = feedback control design for high-precision wafer stage motion,” of MIMO LPV systems using a periodic scheduling sequence,” Auto- Control Engineering Practice, vol. 10, no. 7, pp. 739–755, 2002. matica, vol. 43, no. 10, pp. 1684–1697, 2007. [2] T. Oomen, R. Van Herpen, S. Quist, M. Van De Wal, O. Bosgra, and [26] J. Goos and R. Pintelon, “Continuous-time identification of periodically M. Steinbuch, “Connecting system identification and robust control for parameter-varying state space models,” Automatica, vol. 71, pp. 254– next-generation motion control of a wafer stage,” IEEE Transactions on 263, 2016. Control Systems Technology, vol. 22, no. 1, pp. 102–118, 2014. [3] R. Munnig Schmidt, G. Schitter, and J. van Eijk, The Design of High [27] D. Turk, J. Gillis, G. Pipeleers, and J. Swevers, “Identification of lin- Performance Mechatronics, 1st ed. Amsterdam, The Netherlands: IOS ear parameter-varying systems: A reweighted ` -norm regularization 2;1 Press, 2011. approach,” Mechanical Systems and Signal Processing, vol. 100, pp. [4] T. Oomen, “Advanced motion control for precision mechatronics: Con- 729–742, 2018. trol, identification, and learning of complex systems,” IEEJ Journal of [28] P. B. Cox, “Towards efficient identification of linear parameter-varying Industry Applications, vol. 7, no. 2, pp. 127–140, 2018. state-space models,” Ph.D. dissertation, Eindhoven University of Tech- [5] A. J. Fleming and K. K. Leang, Design, Modeling and Control of nology, 2018. Nanopositioning Systems, 1st ed., ser. Advances in Industrial Control. [29] M. Groot Wassink, M. van de Wal, C. Scherer, and O. Bosgra, “LPV Switzerland: Springer International Publishing, 2014. control for a wafer stage: beyond the theoretical solution,” Control [6] R. van Herpen, T. Oomen, E. Kikken, M. van de Wal, W. Aangenent, Engineering Practice, vol. 13, no. 2, pp. 231–245, 2005. and M. Steinbuch, “Exploiting additional actuators and sensors for nano- [30] M. Steinbuch, R. van de Molengraft, and A. van der Voort, “Experi- positioning robust motion control,” IFAC Mechatronics, vol. 24, no. 6, mental modelling and LPV control of a motion system,” in Proceedings pp. 619–631, 2014. of the 2003 American Control Conference, vol. 2, Denver, Colorado, [7] S. O. R. Moheimani, D. Halim, and A. J. Fleming, Spatial Control of United States, June 2003, pp. 1374–1379. Vibration: Theory and Experiments, ser. Series on Stability, Vibration [31] R. de Rozario, R. Voorhoeve, W. Aangenent, and T. Oomen, “Global and Control of Systems, Series A. Singapore: World scientific, 2003, feedforward control of spatio-temporal mechanical systems: With ap- vol. 10. plication to a prototype wafer stage,” IFAC-PapersOnLine, vol. 50, [8] T. Oomen, E. Grassens, and F. Hendriks, “Inferential motion con- no. 1, pp. 14 575–14 580, 2017, IFAC 20th Triennial World Congress. trol: identification and robust control framework for positioning an Toulouse, France. unmeasurable point of interest,” IEEE Transactions on Control Systems [32] P. Guillaume, P. Verboven, S. Vanlanduit, H. Van Der Auweraer, Technology, vol. 23, no. 4, pp. 1602–1610, 2015. and B. Peeters, “A poly-reference implementation of the least-squares [9] R. van Herpen, T. Oomen, and M. Steinbuch, “Optimally conditioned complex frequency-domain estimator,” in Proceedings of the 21st Inter- instrumental variable approach for frequency-domain system identifica- national Modal Analysis Conference, Kissimmee, Florida, United States, tion,” Automatica, vol. 50, no. 9, pp. 2281–2293, 2014. 2003, pp. 183–192. [10] R. Voorhoeve, R. de Rozario, and T. Oomen, “Identification for mo- [33] B. Cauberghe, “Applied frequency-domain system identification in the tion control: Incorporating constraints and numerical considerations,” field of experimental and operational modal analysis,” Ph.D. dissertation, in Proceedings of the 2016 American Control Conference, Boston, Vrije Universiteit Brussel, 2004. Massachusetts, United States, July 2016, pp. 6209–6214. [34] J. Vayssettes, G. Mercere, ` P. Vacher, and R. De Callafon, “Frequency- [11] M. M. da Silva, W. Desmet, and H. V. Brussel, “Design of mechatronic domain identification of aircraft structural modes from short-duration systems with configuration-dependent dynamics: Simulation and opti- flight tests,” International Journal of Control, vol. 87, no. 7, pp. 1352– mization,” IEEE/ASME Transactions on Mechatronics, vol. 13, no. 6, 1372, 2014. pp. 638–646, Dec 2008. [35] A. J. Fleming and S. O. R. Moheimani, “Spatial system identification [12] J. S. Shamma and M. Athans, “Gain scheduling: potential hazards and of a simply supported beam and a trapezoidal cantilever plate,” IEEE possible remedies,” IEEE Control Systems Magazine, vol. 12, no. 3, pp. Transactions on Control Systems Technology, vol. 11, no. 5, pp. 726– 101–107, June 1992. 736, 2003. [13] W. J. Rugh and J. S. Shamma, “Research on gain scheduling,” Auto- [36] R. Pintelon, P. Guillaume, and J. Schoukens, “Uncertainty calculation matica, vol. 36, no. 10, pp. 1401–1425, 2000. in (operational) modal analysis,” Mechanical Systems and Signal Pro- [14] D. J. Leith and W. E. Leithead, “Survey of gain-scheduling analysis and cessing, vol. 21, no. 6, pp. 2359–2373, 2007. design,” International Journal of Control, vol. 73, no. 11, pp. 1001– 1025, 2000. [37] E. Reynders, “System identification methods for (operational) modal [15] A. Packard, “Gain scheduling via linear fractional transformations,” analysis: Review and comparison,” Archives of Computational Methods Systems & Control Letters, vol. 22, no. 2, pp. 79–92, 1994. in Engineering, vol. 19, no. 1, pp. 51–124, Mar 2012. [16] P. Apkarian and P. Gahinet, “A convex characterization of gain- [38] R. Voorhoeve, N. Dirkx, T. Melief, W. Aangenent, and T. Oomen, scheduled H controllers,” IEEE Transactions on Automatic Control, “Estimating structural deformations for inferential control: a disturbance vol. 40, no. 5, pp. 853–864, May 1995. observer approach,” IFAC-PapersOnLine, vol. 49, no. 21, pp. 642–648, [17] C. W. Scherer, “LPV control and full block multipliers,” Automatica, 2016, 7th IFAC Symposium on Mechatronic Systems. Loughborough vol. 37, no. 3, pp. 361–375, 2001. University, United Kingdom. 13 [39] W. K. Gawronski, Dynamics and Control of Structures: A Modal Robin de Rozario received the M.Sc. degree in Me- Approach, 1st ed., ser. Mechanical Engineering Series. Springer-Verlag chanical Engineering (cum laude) in 2015 from the New York, 1998. Eindhoven University of Technology, Eindhoven, the [40] A. Preumont, Vibration Control of Active Structures: An Introduction, Netherlands, where he is currently pursuing a Ph.D. 4th ed., ser. Solid Mechanics and Its Applications. Springer, 2018, vol. degree in the Control Systems Technology group. 246. His research interests include system identification [41] H.-X. Li and C. Qi, “Modeling of distributed parameter systems for and learning control for high-performance motion applications—a synthesized review from time–space separation,” Jour- systems. nal of Process Control, vol. 20, no. 8, pp. 891–901, 2010. [42] R. Toth, ´ F. Felici, P. S. C. Heuberger, and P. M. J. Van den Hof, “Discrete time LPV I/O and state space representations, differences of behavior and pitfalls of interpolation,” in Proceedings of the 2007 European Control Conference, Kos, Greece, July 2007, pp. 5418–5425. [43] R. Pintelon and J. Schoukens, System Identification: A Frequency Domain Approach, 2nd ed. Hoboken, New Jersey, United States: Wiley- IEEE press, 2012. [44] C. H. Houpis and G. B. Lamont, Digital Control Systems: Theory, Hardware, Software, 2nd ed., ser. McGraw-Hill series in electrical engineering: Control theory. McGraw-Hill, 1992. [45] K. Glover and J. C. Willems, “Parametrizations of linear dynamical systems: Canonical forms and identifiability,” IEEE Transactions on Automatic Control, vol. 19, no. 6, pp. 640–646, 1974. [46] M. de Mathelin and M. Bodson, “Canonical vs pseudo-canonical forms for the structural and parametric identification of multivariable systems,” in Proceedings of the 1st European Control Conference, vol. 2, Greno- Wouter Aangenent received the M.Sc. degree (cum ble, France, 1991, pp. 1282–1287. laude) and the Ph.D. degree in Mechanical Engineer- [47] R. A. de Callafon, D. de Roover, and P. M. J. Van den Hof, “Multi- ing from the Eindhoven University of Technology, variable least squares frequency domain identification using polynomial Eindhoven, The Netherlands, in 2004 and 2008, matrix descriptions,” in Proceedings of the 35th Conference on Decision respectively. In 2008, he joined the Research Depart- and Control, Kobe, Japan, 1996, pp. 2030–2035. ment of ASML, Veldhoven, The Netherlands, where [48] C. K. Sanathanan and J. Koerner, “Transfer function synthesis as a ratio he currently leads the Mechatronics and Control of two complex polynomials,” IEEE Transactions on Automatic Control, Research Group. vol. 8, no. 1, pp. 56–58, 1963. [49] D. S. Bayard, “High-order multivariable transfer function curve fitting: Algorithms, sparse matrix methods and experimental results,” Automat- ica, vol. 30, no. 9, pp. 1439–1444, 1994. [50] A. H. Whitfield, “Asymptotic behaviour of transfer function sythesis methods,” International Journal of Control, vol. 45, no. 3, pp. 1083– 1092, 1987. [51] G. Mercere, ` O. Prot, and J. A. Ramos, “Identification of parameterized gray-box state-space systems: From a black-box linear time-invariant representation to a structured one,” IEEE Transactions on Automatic Control, vol. 59, no. 11, pp. 2873–2885, Nov 2014. [52] G. Wahba, Spline Models for Observational Data, ser. CBMS-NSF Regional Conference Series in Applied Mathematics. Philadelphia, Pennsylvania, United States: Society for Industrial and Applied Mathe- matics, 1990, vol. 59. [53] G. Wahba and J. Wendelberger, “Some new mathematical methods for variational objective analysis using splines and cross validation,” Tom Oomen (SM’06) received the M.Sc. degree Monthly weather review, vol. 108, no. 8, pp. 1122–1143, 1980. (cum laude) and Ph.D. degree from the Eindhoven [54] R. de Rozario, T. Oomen, and M. Steinbuch, “Iterative learning control University of Technology, Eindhoven, The Nether- and feedforward for LPV systems: Applied to a position-dependent mo- lands. He held visiting positions at KTH, Stockholm, tion system,” in Proceedings of the 2017 American Control Conference, Sweden, and at The University of Newcastle, Aus- Seattle, Washington, United States, May 2017, pp. 3518–3523. tralia. Presently, he is associate professor with the Department of Mechanical Engineering at the Eind- hoven University of Technology. He is a recipient of the Corus Young Talent Graduation Award, the IFAC 2019 TC 4.2 Mechatronics Young Research Award, the 2015 IEEE Transactions on Control Systems Technology Outstanding Paper Award, the 2017 IFAC Mechatronics Best Paper Award, the 2019 IEEJ Journal of Industry Applications Best Paper Award, and recipient of a Veni and Vidi personal grant. He is Associate Editor of the IEEE Control Systems Letters (L-CSS), IFAC Mechatronics, Robbert Voorhoeve received the M.Sc. degree in and IEEE Transactions on Control Systems Technology. He is a member of Mechanical Engineering (cum laude) and Applied the Eindhoven Young Academy of Engineering. His research interests are in Physics in 2013 and his Ph.D. degree in Mechanical the field of data-driven modeling, learning, and control, with applications in Engineering in 2018 from the Eindhoven University precision mechatronics. of Technology, Eindhoven, The Netherlands. His research interests include system identification, iden- tification for advanced motion control, and control of complex mechatronic systems.

Journal

Electrical Engineering and Systems SciencearXiv (Cornell University)

Published: Jul 18, 2018

There are no references for this article.