Modeling Complex Nonminimum Phase Zeros in Flexure Mechanisms

This paper presents a model to explain complex nonminimum phase (CNMP) zeros seen in the noncollocated frequency response of a large-displacement XY flexure mechanism, which employs multiple double parallelogram flexure modules (DPFMs) as buildingblocks. Geometric nonlinearities associated with large displacement along with the kinematic under-constraint in the DPFM lead to a coupling between the X and Y direction displacements. Via a lumped-parameter model that captures the most relevant geometric nonlinearity, it is shown that specific combinations of the operating point (i.e., flexure displacement) and mass asymmetry (due to manufacturing tolerances) give rise to CNMP zeros. This model demonstrates the merit of an intentionally asymmetric design over an intuitively symmetric design in avoiding CNMP zeros. Furthermore, a study of how the eigenvalues and eigenvectors of the flexure mechanism vary with the operating point and mass asymmetry indicates the presence of curve veering when the system transitions from minimum phase to CNMP. Based on this, the hypothesis of an inherent correlation between CNMP zeros and curve veering is proposed. [DOI: 10.1115/1.4036032]


Introduction and Motivation
This research investigation is motivated by the need to achieve large range, high precision, and high-speed-all simultaneouslyin multi-axis flexure-based motion systems [1,2].Such capability is of practical importance in various applications such as compact and affordable motion stages for semiconductor wafer inspection [3] and microelectromechanical system scanners for high-speed imaging [4].Flexure mechanisms are well-suited for these applications because of their joint-less construction and inherently high precision due to lack of friction and backlash, but present significant tradeoffs between large displacement and dynamic performance [5].
Large displacement generally implies transverse deformation of the constituent beams in the flexure mechanism greater than 5% of the beam length.This corresponds to several millimeters of displacement or motion range for a desktop size flexure-based motion system.The relevant system dynamics include natural frequencies, mode shapes, and transfer functions between the points of actuation and sensing.The closed-loop dynamic performance objectives include high bandwidth, good noise and disturbance rejection, good command tracking, small steady-state error, fast point-to-point positioning and settling, stability robustness, low sensitivity to plant variations, etc.
While recent results have demonstrated large range as well as high precision in multi-axis flexure mechanisms, achieving dynamic performance remains a challenge [5,6].Figure 1(a) shows an XY nanopositioning system based on a parallel-kinematic flexure mechanism, designed to achieve a range of 10 mm and precision of 625 nm per axis.This mechanism employs a systematic and symmetric layout of eight double parallelogram flexure module (DPFM).This design provides: a high degree of geometric decoupling between the X and Y motions of the motion stage resulting in large unconstrained motion range; actuator isolation that allows the use of large-stroke single-axis actuators (X actuator and Y actuator); and a complementary endpoint sensing using commonly available sensors (e.g., sensors 1 and 2 for the X direction).For reference, relevant dimensions of this flexure mechanism are summarized in Table 3 in the Appendix.
There are many factors that make the dynamics of such a flexure mechanism challenging.Large displacements lead to geometric nonlinearities in flexure mechanics.Given their dependence on displacement, these nonlinearities (and their impact on system dynamics) vary with the operating point of the flexure mechanism.Furthermore, large displacements require relatively low stiffness in the motion directions and therefore low natural frequencies of the first few modes.Any attempt to achieve a closed-loop bandwidth that is greater than these low natural frequencies requires a proper understanding of higher-order dynamics, which is complicated by the above-mentioned geometric nonlinearities.Furthermore, while the symmetric layouts (e.g., Fig. 1(a)) help provide large range, cancel undesired motions, improve space utilization, and enhance quasi-static performance, they also result in multiple closely spaced modes that are highly sensitive to manufacturing tolerances.This results in parametric uncertainty in the system dynamics.
Figure 1(b) shows an experimental measurement of the noncollocated X direction frequency response from force input F x to the displacement output X ms , for different values of Y displacement (Y ms ¼ 0, 1.5, 3 mm) [5,7].One may notice that there are multiple closely spaced modes around 150 Hz.These correspond to the natural frequency of the secondary mass in each DPFM (discussed further in Sec. 3).It is also noteworthy that the X direction frequency response changes with the Y direction operating point.At the operating points y ms ¼ 1.5 and 3 mm, the frequency response shows an additional phase drop of 360 deg and 720 deg, respectively, near 150 Hz compared to the nominal operating point y ms ¼ 0. The magnitude and phase below and above 150 Hz remains the same for all the operating points.Such observation cannot be explained by minimum phase zero pairs.Thus, such phase drop is due to complex nonminimum phase (CNMP) zero pairs on the right half plane.This dynamic response was unexpected, and the existence and number of CNMP zeros seemed arbitrary.From closed-loop performance stand-point, it is well known that nonminimum phase (NMP) zeros severely limit bandwidth, stability robustness, and positioning speed [8,9].When and why do these CNMP zeros appear?Can they be analytically predicted?Do they have physical meaning?Can they be avoided via physical system design?Addressing these questions is the motivation behind this investigation.
The paper is organized as follows: Section 2 provides an overview of the relevant literature on modeling geometric nonlinearities in flexure mechanics and NMP zeros in the dynamics of flexible systems.Section 3 presents closed-form kinematic relations to capture the relevant nonlinearity and coupling in the DPFM.In Sec. 4, these relations are employed in initially investigating the X direction dynamics of a simple representative XY flexure mechanism for different Y operating points.This simple mechanism has all the essential attributes of the more complex XY flexure mechanism of Fig. 1(a), but is more conducive for an initial investigation.A closed-form, parametric dynamic model helps predict the range of operating points and parametric asymmetry where CNMP zeros appear in a noncollocated transfer function of this simple flexure mechanism.Similar modeling and CNMP prediction is then extended to the XY flexure mechanism of Fig. 1(a).Section 5 explores a potential correlation between these CNMP zeros and the phenomenon of curve veering, which lends some physical insight into the former.The paper concludes in Sec.6 with a list of contributions and future tasks.One of the key findings is that the intentional use of specific parametric asymmetry, which is counter-intuitive, helps avoid the CNMP zeros altogether.

Literature Review
In recent years, there has been a growing body of research literature on the dynamics of flexure mechanisms.Lan et al. [10] presented a distributed-parameter dynamic modeling approach of elastic flexure mechanisms.The resulting equations of motion in the time-domain were solved using numerical methods, which are not very suitable for frequency domain analysis of complex NMP zeros.Akano and Fakinlede [11] used finite element-based nonlinear analysis to predict the effect of design parameters on the dynamic performance of flexure mechanisms.While accurate, these methods are computationally intensive and provide limited physical insights in the frequency domain.Alternatively, lumpedparameter closed-form modeling approaches have also been investigated.Shilpiekandula and Youcef-Toumi [12] derived a lumped-parameter dynamic model of a diaphragm flexure using Timoshenko beam theory, but did not include geometric nonlinearities.Awtar and Parmar [5] captured the nonlinear variation in the stiffness of flexure building-blocks to create a lumpedparameter dynamic model of a XY flexure mechanism (Fig. 1(a)), but did not capture the nonlinear coupling between X and Y directions in a DPFM and therefore were unable to predict the NMP behavior seen experimentally (Fig. 1(b)).The pseudorigid-body approach has also been used for modeling the nonlinear dynamics of flexure mechanisms [13][14][15].While this approach leads to simple lumped-parameter closed-form models, the model parameters are computed via numerical optimization and depend on the boundary conditions of each beam, thereby increasing the modeling complexity in flexure mechanisms that have a large number of beams.
Dynamic modeling of rigid link mechanisms with inherent flexibilities, e.g., robotic manipulators, has also been an active area of research.An overview and classification of various modeling approaches is found in the review paper by Dwivedy and Eberhard [16].Research in this area includes the study of manipulators with one or more flexible links as well as one or more flexible joints.Various methods including finite elements, assumed modes, lumped parameter, and inverse dynamics have been adopted to study the relevant dynamics.This body of work assumes small deformation of the links, compared to rigid body motion, which is justified since the links are designed to be stiff.However, this assumption fails for flexure mechanisms that provide large deformation in their motion directions.
The large deformation of constituent elements or beams in a flexure mechanism results in geometric nonlinearities arising from arc-length conservation, cross-sectional warping, trapeze effect, and Wagner's effect in beam mechanics [17][18][19][20][21][22][23].The impact of these nonlinearities on the dynamics of flexible beams and structures has been studied extensively, as reported in the review papers by Modi [24] and Pandalai [25].Furthermore, the dynamics and control of flexible beams with an end-mass [26] as well as rotating beams [27] have also been investigated.In further generalization, DaSilva formulated the nonlinear differential equations of motion for Euler-Bernoulli beams experiencing flexure along two principal directions, along with torsion and extension [28].Jonker has formulated a highly generalized model for spatial beams taking into account relevant nonlinearities, using finite element-based multibody dynamics computations [19,29].Nayfeh modeled the nonlinear transverse vibration of beams with properties that vary along the length [30].Zavodney and Nayfeh studied the nonlinear response of a slender beam with a tip mass to a principal parametric excitation [31].Moeenfard and Awtar studied the in-plane flexural and axial vibration of a flexure beam with a tip mass while accounting for the nonlinearity associated with arclength conservation [32].While the resulting nonlinear equations Transactions of the ASME of dynamics are solved in time-domain via perturbation, homotopy, or computational methods, this prior work [32] does not pursue the frequency domain investigation relevant to the present work.
Separately, there exists a significant body of work in the frequency domain dynamics of lumped or distributed-parameter flexible systems [33][34][35].It has been shown that lightly damped flexible systems with collocated sensor and actuator have alternating poles and zeros along the imaginary axis and are easy to stabilize in closed-loop [36,37].Noncollocated systems do not share these attributes and, under certain conditions, exhibit real nonminimum phase (NMP) zeros in the right half plane [38][39][40].Spector and Flashner [38] studied the sensitivities of beam cross section, material properties, and sensor placement on the locations of poles and zeros in flexible systems.They showed that as the sensor placement is moved away from the actuator, the conjugate zeros, originally located along the imaginary axis, migrate toward infinity and then reappear along the real axis.Miu [39] provided a physical explanation for these real NMP zeros stating that they are related to the nonpropagation of energy within the structural subsystem confined by the actuator and sensor.Unlike real NMP zeros, CNMP zeros are relatively rare and have been reported in the context of a noncollocated acoustical transfer function of a room [41], as well as in a noncollocated transfer function of a lumped-parameter spring-mass system [42,43].Awtar and Craig identified CNMP zeros arising due to an electromagnetic coupling between a direct current motor and tachometer used in a servosystem [44].These studies on CNMP zeros simply report a mathematical or experimental observation, without providing further insight into when or why the zeros appear.

Modeling Geometric Nonlinearity in Double Parallelogram Flexure Modules
A DPFM, shown in Fig. 2, comprises a reference stage, a secondary stage, and a primary stage connected via two parallelogram flexures in series.This arrangement uses geometric reversal to cancel the X direction kinematic error motion (i.e., parabolic trajectory) of one parallelogram with that of the second parallelogram.In quasi-static operation, these two error motions exactly cancel out, resulting in a straight-line motion along the Y direction between the primary and reference stages while providing high stiffness in the X and H axes. Furthermore, the DPFM provides large-displacement range and low stiffness in the Y direction for a given footprint.For a building block, these are desirable attributes because they help minimize various error motions and enable large motion range in the resultant flexure mechanisms [5].
One of the limitations of the DPFM is that it presents a kinematic under-constraint associated with the Y displacement of its secondary stage.In quasi-static operation, this under-constraint adversely impacts the X direction stiffness, as reported previously [23].In dynamic operation, this leads to an additional degrees-offreedom (DOF) in the DPFM even though it is intended to be a single-DOF building block.Referring to the X and Y displacements of the primary and secondary stages with respect to the reference stage in Fig. 2, the following quasi-static relations hold [7,23]: (1) Here, L stands for the beam length (Fig. 2).These kinematic relations, which define the X direction error motion of the parallelogram flexures, are the result of the geometric nonlinearity associated with beam arc-length conservation.Relative to the reference stage, there are four displacement coordinates in the DPFM (i.e., X i , Y i , X j , and Y j ); of these, only two are independent because of the above kinematic relations.This illustrates the 2DOFs of the DPFM and its under-constraint.
The above quasi-static relations are based on certain core assumptions that may be extended to low-frequency dynamics spanning the resonance of the secondary stage of a DPFM with respect to the reference and primary stages.In the flexure mechanism of Fig. 1, this is the frequency range in which complex nonminimum phase zeros appear.Since the secondary stage mass is small compared to the other stages and the X direction stiffness of the beams is relatively high, modes associated with beam stretching appear at much higher frequencies.This observation justifies the first assumption that the individual flexure beams may be treated as inextensible.Beam inextensibility also eliminates H displacement coordinates in the DFPM.Furthermore, in the frequency range of interest, the beams are also assumed to be massless, given the relatively larger masses of the various stages.This implies that the beams deform in the quasi-static S-shape with amplitude dictated by the relative Y displacement of the relevant parallelogram flexure [22,23].At higher frequencies, there will be other beam shapes dictated by the resonance of individual beams, in which case the above relations will no longer hold.
Because of large out-of-plane stiffness, the out-of-plane displacements (and associated modes) are neglected in the dynamic modeling.Also, damping is neglected since the flexure mechanism is monolithic with no sliding joints.Experiments confirm that the dissipation from material hysteresis and viscous effects is very low (damping ratios < 0.001).Beyond arc-length conservation, there are several other nonlinearities in beam mechanics [20,21], but these are systematically estimated and neglected using an order of magnitude analysis [7,45].

Dynamic Modeling of XY Mechanisms
The kinematic nonlinearity of Eqs. ( 1) and ( 2) may be incorporated in deriving the equations of motion for a flexure mechanism involving the parallelogram or double parallelogram flexure modules.To investigate how the frequency domain dynamic response varies with the operating point, these nonlinear relations may be linearized about an arbitrary operating point.But this potentially holds the risk of premature linearization.To test this possibility, we retain the nonlinear kinematic relations throughout the derivation of the dynamic equations of motion and linearize the latter in the end.The results prove to be the same as when the kinematic relations themselves are linearized at the onset [45].
Thus, defining Y ij ¼ Y i À Y j , using the subscript "o" to denote nominal values at an operating point, and lower case letters to denote deviation from these nominal values, Eq. ( 1) becomes Since the nominal values are still related by Eq. ( 1), Eq. ( 3) may be linearized for small deviations about the nominal values Here, a is a coupling coefficient that depends on the operating point Y ijo and captures the coupling between the X and Y axis displacement coordinates.
Next, for the purpose of investigating CNMP zeros in XY flexure mechanisms that have multiple DPFM building-blocks, we select a simple representative flexure mechanism shown in Fig. 3 to initially limit modeling complexity and enable physical insights into the observed dynamic phenomena.Yet, this design captures all the essential attributes of more complex flexure mechanisms (e.g., Fig. 1).The layout comprises the smallest number of symmetrically placed DPFMs (i.e., two) needed to produce multiple (i.e., two) closely spaced modes associated with the kinematic under-constraint of the secondary stages (two and three).This mechanism allows stage ‹ to displace in the X and Y directions.The former is due to the X direction bearings (indicated by the rollers) at stages fl and °, and the latter is due to large bending deformations of the beams in the two DPFM.Large deformation leads to the geometric nonlinearity and associated coupling between X and Y displacements, mentioned earlier.For this mechanism, the operating point is given by a static displacement of stage 1 in the Y direction with respect to ground (Y 1o ) caused by a constant force F Y1o .Thus, the noncollocated transfer from an X direction force P on stage fl to the X displacement of stage ‹ (X 1 ) can be investigated for different values of Y 1o .In fact, this simple mechanism is representative of a portion (indicated by the larger dashed rectangle) of the more complex XY mechanism of Fig. 1(a).
The five stages have eight displacement coordinates, as shown in Fig. 3; lower case versions of these coordinates represent the respective deviations with respect to an operating point.Furthermore, these coordinates are related by four kinematic relations (Eq.( 4)), one for each parallelogram.Thus, this mechanism has 4DOF and is therefore referred to as the simple representative 4DOF (SR4DOF) mechanism in this paper.The displacement coordinates x 1 , y 1 , y 2 , and y 3 are chosen for this analysis, to be able to study the X displacement of stage ‹ and the Y displacements of stages › and fi with respect to stage ‹.
Assuming a lumped-parameter Y direction bending stiffness for each of the parallelograms (k 12 , k 24 , k 13 , and k 35 ), the equations of motion for the SR4DOF may be derived to be the following: It should be noted that the lower case displacement coordinates in the above equations of motion represent the small deviations about the respective nominal operating point values.Based on these equations, one can derive the transfer function G(s) from the input force P to output displacement x 1 for different values of a, which depends on Y 1o In the above stiffness matrix, the Y direction bending stiffness for each parallelogram (k 12 , k 24 , k 13 , and k 35 ) is nominally 24EI/L 3 , where E is the bending modulus, L is the beam length, and I is the second moment of area about the Z axis [23].To obtain numerical results, we use the same dimensions as those for the XY flexure mechanism of Fig. 1(a) (see Table 3 in the Appendix).Furthermore, although not included in the above derivation, small nominal values for damping are assumed to avoid singularities in the numerical simulation.Although the SR4DOF is intended to be symmetric, there is the possibility for parametric asymmetry between (m 4 and m 5 ), (m 2 and m 3 ), (k 4 and k 5 ), or (k 12 , k 24 , k 13 , and k 35 ) resulting from finite manufacturing tolerances.Initially, we assume the parameters be perfectly symmetric; but the above lumped-parameter model allows us to study the impact of asymmetries in Sec.4.2.
The last three of the four predicted modes are shown in Fig. 4 (relative to the displaced configuration/operating point of Fig. 3), while all modes are quantified in Table 1.There are three key observations: (1) The first mode (not shown in Fig. 4) is the "rigid body" mode in which all the stages vibrate together in the X direction due to springs k 4 and k 5 ; (2) the second mode is associated The third and fourth modes arise due to the under-constrained secondary stages in the DPFM.When Y 1o ¼ 0 (i.e., a ¼ 0), the vibration of the secondary stages does not cause any X direction motion of stage ‹.Therefore, these two modes are unobservable in the G(s) transfer function.However, when Y 1o 6 ¼ 0, the X and Y displacements of the DPFMs get coupled, which affects the third and fourth modes differently.For the third mode, the Y vibration of the two secondary stages is coupled to the X vibration of stage ‹.However, in the fourth mode, the two secondary stages have the same vibration magnitude and phase in the Y direction, which results in a cancelation of the coupling at stage ‹ in the X direction.Instead, the coupling results in X direction vibrations at stages fl and °(Fig.4(c)).Thus, when parameters are symmetric but Y 1o 6 ¼ 0, the third mode shows up in the G(s) transfer function, while the fourth mode remains unobservable.

Parametric Asymmetry and Complex Nonminimum
Phase Zeros.While the SR4DOF was assumed to be perfectly symmetric so far, some parametric asymmetry is inevitable in a practical situation due to manufacturing tolerances.Using the above lumped-parameter model, we varied parametric asymmetry over a 65% range for (m 4 and m 5 ) and (k 4 and k 5 ).This did not have much of an effect on the flexure mechanism dynamics in terms of mode shapes and transfer functions.But an asymmetry in (k 12 , k 24 , k 13 , and k 35 ) or (m 2 and m 3 ) impacts the vibrations of the secondary stages (i.e., third and fourth modes) and therefore the overall flexure mechanism.Of the two sets of parameters, the mass parameter is more sensitive.When the DPFM is used as a building block, its secondary stage size and mass are minimized to reduce footprint and raise the resonance frequency at which it vibrates.For example, in the XY flexure mechanisms of Figs.1(a) and 3, the nominal mass of the secondary stage is 18 g.Therefore, even a small additional mass such as 0.9 g results in a relatively large variation (5%).Therefore, in this section, we investigate how an asymmetry in masses, Dm 23 (¼m 2 /m 3 À 1), affects the dynamics of SR4DOF.
As seen via the respective eigenvectors in Table 1, the impact of nonzero Dm 23 on the first and second modes is minimal.This mass asymmetry primarily impacts the vibration of the two secondary stages, which directly influence the third and fourth modes.As noted earlier, the fourth mode of the SR4DOF flexure mechanism is unobservable in the transfer function G(s) for any Y 1o (zero or nonzero) when Dm 23 ¼ 0. However, for a small parametric asymmetry, e.g., Dm 23 ¼ 5%, the two secondary stages have different vibration magnitudes as seen in the eigenvector in Table 1.Thus, the X direction coupling no longer cancels out at stage ‹, and the fourth mode appears in G(s).Similarly, the impact of Dm 23 ¼ 5% on the third mode is significant.
G(s) is plotted in Fig. 5 as Y 1o varies from 0% to 5%, for Dm 23 ¼ 0% and 5%.Key observations are: (1) As expected, the fourth mode is unobservable when Dm 23 ¼ 0 but appears when Dm 23 ¼ 5%.(2) The natural frequencies of the third and fourth modes drop as Y 1o increases.As the X vibration of stage ‹ gets increasingly coupled with the Y vibration of the secondary stages, the modal mass increases more than the modal stiffness, resulting in reduced natural frequencies.(3) A 360 deg phase drop is observed at around 150 Hz in the asymmetric case (Dm 23 ¼ 5%) when Y 1o ¼ 5% but not in a symmetric case (Dm 23 ¼ 0%).In the latter case, the phase drop due to the complex pole pair (third mode) is offset by a phase rise due to the complex zero pair ("valley"), resulting in no net phase drop.For the asymmetric case, there are two stable complex pole pairs (third and fourth modes), each contributing 180 deg phase drop.But since the overall phase drop is 360 deg, this implies that there is no phase rise or drop at the valley even though there are two pairs of complex zeros at in this frequency region ($153 Hz).This indicates the presence of a quartet of complex zeros, with one pair in the left half plane and the second pair in the right half plane, thereby contributing no net  phase change.This is an important observation because it suggests that CNMP zeros can arise at certain combinations of operating points and parametric asymmetry.
4.2 Existence of Complex Nonminimum Phase Zeros.Next, we proceed to analytically determine the conditions under which CNMP zeros arise.Based on modal decomposition, the transfer function G(s) may be written as follows, where b i is the modal residue and x i is the corresponding natural frequency: The decomposed form can also be expressed via a numerator N(u) and a denominator D(u), each a polynomial.Note that there are no odd power s terms because damping is ignored.If N(u), which is a cubic polynomial, has two complex conjugate roots, then G(s) will have a quartet of complex zeros.Two of these zeros will be in the right half plane (i.e., CNMP zeros).For this to happen, the following inequality in the coefficients of N(u) has to hold [46]: Therefore, this is the mathematical condition for the existence of CNMP zeros in the G(s) transfer function.Shown in Fig. 6, D is plotted in a contour map against a range of operating points and parametric asymmetry values for the SR4DOF flexure mechanism.The color in the contour map represents the magnitude of D: red represents higher positive values, blue represents lower positive values, and the black region represents the conditions for which D becomes negative, indicating the presence of CNMP zeros.This particular mechanism is seen to be very sensitive to positive asymmetry, i.e., if m 2 is greater than m 3 even by a small amount, then CNMP zeros arise in specific ranges of Y 1o .However, if m 2 < m 3 , then the entire operating range is free of CNMP zeros.The reason for such asymmetric behavior is due to the physical asymmetry introduced by the actuator placement in Fig. 3.
With this finding, we are able to replicate via modeling some aspects of the NMP phenomenon previously observed experimentally (Fig. 1(b)).Although, for this study, we intentionally chose the SR4DOF mechanism to keep modeling complexity and assumptions minimal, it is representative of the more complex designs in that it incorporates the key attributes of DPFM building-blocks (with their under-constrained secondary stages), geometrically symmetric design, large displacements leading to nonlinear coupling between axes, noncollocated transfer functions, and parametric asymmetry.
In the design of multi-DOF flexure mechanisms, it is a common guideline to employ symmetric and/or periodic geometries to cancel undesired motion, improve space utilization, and enhance quasi-static performance [1].However, the above dynamic model for the SR4DOF mechanism indicates that a perfectly symmetric layout is sensitive to parametric asymmetry, which is likely to occur due to manufacturing tolerances and can give rise to CNMP zeros.But if the design is intentionally made asymmetric (i.e., if Dm 23 is sufficiently negative), then CNMP zeros can be avoided even if there are finite manufacturing tolerances.
This also shows that one can choose the springs in the flexure mechanism (i.e., beam flexures dimensions and layout) to be symmetric to achieve the desired quasi-static performance, while choose certain masses to be asymmetric which provides the desired dynamic performance without impacting quasi-static performance.This combination of symmetry in certain attributes and asymmetry in others is rather counter-intuitive but helps meet both quasi-static and dynamic performance goals.
With CNMP zeros thus eliminated via physical design, we also create the possibility of achieving closed-loop bandwidth higher than 150 Hz, while maintaining robustness, in the SR4DOF flexure mechanism.This would have been impossible in the presence of CNMP zeros at around 150 Hz.
4.3 Modeling a Complex XY Mechanism.So far, we modeled the SR4DOF flexure mechanism to predict CNMP zeros at around 150 Hz.Next, we extend this modeling approach to the more complex (and practically relevant) XY flexure mechanism of Fig. 1(a).The transfer function X actuator force F x to the X displacement of the motion stage X ms is investigated.The Y actuator is used to provide a constant force F yo to achieve various Y direction operating points Y mso .There are 13 rigid stages in this case, each with an X and Y displacement coordinate.Each of the 16 parallelogram flexures provides one kinematic relation between relative X and Y coordinates.Therefore, the model has ten independent DOFs, which results in the same number of equations of motion, natural frequencies, and mode shapes.There is a rigid body X vibration mode, a rigid body Y vibration mode, and eight modes associated with the vibration of the secondary stages (all around 150 Hz) in the eight DPFM.Next, we arbitrarily vary these secondary stage masses with respect to their nominal value and use the model to predict the 0 deg/360 deg/720 deg phase drop  We present three cases in Table 2 with different combinations of secondary stage mass variations that result in the three different phase drops.Thus, we are able to analytically predict the seemingly unexplained phenomena observed experimentally in Fig. 1(b).

Possible Relation Between Complex Nonminimum Phase Zeros and Curve Veering
As noted earlier, the SR4DOF mechanism considered in this paper has a symmetric and repetitive geometry, which leads to multiple closely spaced modes (i.e., the third and fourth modes shown in Table 1).Furthermore, these modes vary with the operating point and parametric asymmetry.The operating point determines the extent of the cross-axis coupling between X and Y displacements, thus building a connection between the third and fourth modes.All these features make the phenomenon of curve veering (or mode veering) potentially relevant in this study.
Curve veering occurs when the eigenvalue loci of two closely spaced modes in a system approach each other and then diverge, as a result of parameter variation [47].The point in the parameter space where the two modes are the closest is called the veering point.The special case when the two modes intersect at the veering point is called mode crossing [48].In the vicinity of the veering point, eigenvectors undergo dramatic changes.As a result, the system could become so-called "critically configured" meaning that small changes in a system parameter could cause large changes in system response [49].
In the SR4DOF flexure mechanism dynamics, we observe curve veering close to the operating point and parametric asymmetry value at which the CNMP zeros arise.Figure 7(a) shows the eigenvalue loci of the third and fourth modes as a function of the operating point (Y 1o ) and parametric asymmetry (Dm 23 ).When the structure is completely symmetric (i.e., Dm 23 ¼ 0), the two loci intersect with each other at a frequency of 153 Hz and Y 1o ¼ 6.6%.Graphically (Fig. 7(a)) as well as qualitatively, this is mode crossing.When the structure is asymmetric (e.g., Dm 23 > 0), the two loci do not intersect.Instead, they approach each other and then diverge as Y 1o increases.Thus, the point when Y 1o ¼ 6.6% is the veering point.
In the third and fourth modes, the Y direction vibrations of the two secondary stages (i.e., y 2 and y 3 ) are dominant in the respective eigenvectors (see Table 1).Therefore, only the evolution of y 2 and y 3 are plotted in Fig. 7. Figure 7(b) shows such evolution of the third mode.When Dm 23 ¼ 0 and Y 1o ¼ 0, y 2 is 0.7 and y 3 is À0.7 (i.e., opposite directions).When Y 1o reaches the crossing point, the value of y 3 suddenly changes from À0.7 to 0.7 and stays constant with y 2 unchanged (i.e., y 2 and y 3 become in the same direction); in other words, modes 3 and 4 "swap" at the crossing point.This transition is more gradual (i.e., without a discontinuity) for nonzero Dm 23 .Similarly, as shown in Fig. 7(c), in the fourth mode, y 2 and y 3 change from being in the same direction to opposite directions.Note that mode shapes shown in Fig. 4 correspond to those before curve veering has happened.
To summarize, the eigenvalues and eigenvectors of the third and fourth modes exhibit veering.Furthermore, for the asymmetric case  (Dm 23 > 0) at the veering point, y 3 becomes zero while y 2 becomes one in the third mode, while opposite happens in the fourth mode.This is recognized to be mode localization [50].
The above analysis shows that curve veering exists in the SR4DOF flexure mechanism dynamics.However, the more important observation here is that the curve veering happens at about the same operating point (Y 1o ¼ 6.6%), positive parametric asymmetry (Dm 23 > 0), and frequency (153 Hz) at which CNMP zeros appear (see Figs. 5 and 6).Moreover, the CNMP phenomenon and curve veering share the same key factors such as closely spaced modes, parametric asymmetry, and mode coupling (caused by operating point variation in this case).Therefore, we hypothesize that this is not merely a coincidence, and that there exists a fundamental relationship between these two phenomena.If established, this would allow a new physical perspective and interpretation of the CNMP phenomenon and may help guide physical system design.

Contributions and Conclusion
The key contributions of this paper are: (1) A lumpedparameter modeling approach is proposed to analytically model the dynamics of flexure mechanisms comprising the parallelogram or double parallelogram modules.This model captures the key relevant geometric nonlinearity in large-displacement flexure mechanics.Linearization about any arbitrary operating point enables frequency domain analysis.( 2) Based on this model, we are able to predict previously unexplained CNMP zeros seen experimentally.The model establishes the existence of CNMP zeros under certain combinations of operating point and parametric asymmetry in the noncollocated transfer function of a simple representative XY flexure mechanism.(3) This finding helps generate the design insight that, rather than an intuitively symmetric design, an intentional asymmetry in mass can avoid CNMP zeros and make the system conducive to better dynamic performance.
In addition, there are several new questions posed by this work that are currently being addressed: (1) Experimental validation of the analytical predictions such as the CNMP map (Fig. 6) that can validate the modeling simplifications and assumptions.(2) While CNMP zeros were predicted for the SR4DOF as well as the full XY flexure mechanisms, these results are mathematical; greater physical insight into what causes the CNMP zeros is desirable.(3) The potential correlation between CNMP zeros and curve veering was based on observations in this paper but needs to be investigated scientifically.(4) Based on the findings of this paper, physical and control system design can be explored to achieve the originally stated dynamic performance goals of large range, high precision, and high-speed.
This overall investigation potentially has relevance not just to the XY flexure mechanisms considered in this paper but also to a broader range of flexible systems.

Fig. 4
Fig. 4 Mode shapes of the SR4DOF.The black arrows indicate the relative motion direction of each stage.(a) Second mode, (b) third mode, and (c) fourth mode.

Fig. 6
Fig. 6 Contour map of D function

Fig. 5 G
Fig. 5 G(s) transfer function for different operating points (Y 1o ) and mass asymmetry (Dm 23 )

Fig. 7 (
Fig. 7 (a) Eigenvalue loci of the third and fourth modes and (b) and (c) eigenvector elements of the third and fourth modes.The sign of magnitude indicates the phase.

Table 2
Deviations of masses that cause different phase drops