Sean Bryan
email: sean.a.bryan@asu.edu

James Aguirre

George Che

Simon Doyle

Daniel Flanigan

Christopher Groppi

Bradley Johnson

Glenn Jones

Philip Mauskopf

Heather Mccarrick

Alessandro Monfardini

Tony Mroczkowski

WSPEC: A Waveguide Filter-Bank Focal Plane Array Spectrometer for Millimeter Wave Astronomy and Cosmology

Keywords: millimeter waves, spectroscopy, kinetic inductance detectors

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Large area millimeter wave spectral surveys with moderate spectral resolution can be used to characterize the large scale structure and star formation history of the universe using intensity mapping of emission lines such as CO [START_REF] Lidz | Intensity mapping with carbon monoxide emission lines and the redshifted 21 cm line[END_REF] and CII [START_REF] Silva | Prospects for detecting cii emission during the epoch of reionization[END_REF] . Multi-object mm-wave wide-band spectroscopy with moderate spectral resolution will enable galaxy redshift surveys. Further studies of hot gas in galaxy clusters through the Sunyaev Zeldovich (SZ) effect will be enabled by high angular resolution and moderate spectral resolution instruments. Achieving these science goals requires a millimeter wave instrument with moderate resolving power, single-dish-class angular resolution, and high instantaneous sensitivity. Z-spec [START_REF] Inami | A broadband millimeter-wave spectrometer Z-spec: sensitivity and ULIRGs[END_REF] , an R ∼ 100 single-pixel waveguide grating spectrometer, has demonstrated photon noise limited performance in the atmospheric window from 200-300 GHz. Many programs, such as SuperSpec 4 , Micro-Spec 5 , DESHIMA [START_REF] Endo | Development of DESHIMA: a redshift machine based on a superconducting on-chip filterbank[END_REF] , and TIME [START_REF] Staniszewski | The Tomographic Ionized-Carbon Mapping Experiment (TIME) CII Imaging Spectrometer[END_REF] , are currently working towards on-chip spectrometer technology to enable focal plane arrays. In this paper, we describe a waveguide filter bank spectrometer, WSPEC, which can be fabricated with existing highprecision machining techniques and is an alternative technology for enabling filled focal plane array spectrometers. The detectors are single-layer aluminum LEKIDs 8 similar to those described by Heather McCarrick and Daniel Flanigan elsewhere in these proceedings. Aluminum LEKIDs have successfully been used on the IRAM telescope for science observations with the NIKA instrument 9 , and have also been sucessfully built and demonstrated for the NIKA2 instrument as described by Johannes Goupy elsewhere in these proceedings.

The spectrometer concept is shown in Figure 1, which is a drawing of a prototype WR10 device we built and tested [START_REF] Bryan | A compact filter-bank waveguide spectrometer for millimeter wavelengths[END_REF] . Broadband light from the sky is coupled into a main waveguide using a feed horn. At an individual spectrometer channel, there is a narrow coupling section of waveguide, which is narrow enough that its cutoff frequency is 50% higher than the device operating frequency. This is followed by a half-wave resonant section of full-width waveguide, and another narrow section to define the end of the resonator. On resonance, the evanescent coupling through the narrow sections is effectively capacitive, and the resonance in the half-wave section is effectively inductive, which means that overall the channel is matched with the main waveguide and power passes through to that channel's detector. Off-resonance, this cancellation does not take place, which keeps out of band radiation from coupling to the detector. The design of a full 54-channel device that fills the single-moded passband of rectangular waveguide is discussed in an earlier publication [START_REF] Bryan | A compact filter-bank waveguide spectrometer for millimeter wavelengths[END_REF] .

Including Detector Effects in the Model

The design and performance modeling of the spectrometer channels was performed with HFSS. This yields a three-port S-matrix model of each spectrometer channel, which we cascaded together into a complete spectrometer model [START_REF] Bryan | A compact filter-bank waveguide spectrometer for millimeter wavelengths[END_REF] . Since it is purely a S-matrix model, this yielded a performance estimate that assumes perfect terminations at all ports of the device. However, even with high optical-efficiency LEKID detectors, the termination at the detector ports of the real device will not be perfect.

To model this imperfect termination effect in the integrated 54-channel device, and verify that it will not significantly impact performance, we took the three-port individual-channel HFSS models, and terminated the detector port in each device with a HFSS simulation of the input impedance of a LEKID detector. This turned the three-port model of a single channel into a two-port model with some loss due to the power absorbed by the detector. In addition, this two-port model also includes the effect of any small reflections caused by imperfect absorption by the LEKID detector.

We developed an approach to track the power absorbed at each individual detector in the integrated device. Consider the ith detector in the model. We cascaded the individual spectrometer channels above that detector into a combined transmission line matrix ABCD subsequent , and cascaded the others, from the first to the ith, into a matrix ABCD be f ore . We then converted the transmission line matrix ABCD subsequent into an S-matrix, which yielded the input impedance looking into these channels. Using this impedance at port 2, and the nominal waveguide impedance at port 1, we then converted ABCD be f ore to its equivalent Smatrix. Absorption in a detector a i , transmission on to the subsequent detectors Fig. 3 Lossless model (left) and measured (right) performance of the WR5 prototype. The two broadly agree, except for an observed decrease in optical efficiency that is consistent with the known impact of waveguide loss. This will be reduced in the science device due to its smaller size, and we may coat the device in superconducting niobium to eliminate waveguide loss in our band. The source power is very low above about 205 GHz, so the apparent increase in coupling in the measured data is likely a measurement artifact.

t ≡ S 21 × conj(S 21 ), or reflection from the entire structure r ≡ S 11 × conj(S 11 ) are the only things that can happen in the first i detectors. So, we simply subtracted the transmission and reflection to leave the absorption, that is,

a i = 1 -r -t.
This quantity a i is actually the combined absorption of all the detectors from the first to the ith. To get a model for an individual detector's absorbed power, we repeated the calculation for the detectors from the first to the (i -1)th to calculate a i-1 . The power absorbed at the ith detector alone is therefore a ia i-1 . We then repeated this calculation for each channel in the spectrometer, as well as sweeping the simulation frequency, to yield a complete model of the integrated spectrometer system and its LEKID detectors. The results are shown in Figure 2. In the end, the overall optical efficiency of the device remains high, and the effect of imperfect detector terminations is not calculated to induce significant out of band coupling or other problems in the spectrometer. This modeling motivates moving forward with building a device with integrated LEKIDs.

WR5 Prototype Test Results

To test that the spectrometer concept will work at the ∼ 200 GHz frequencies that are needed for science observations, we built a five-channel prototype device in WR5 waveguide. The device is nearly identical to the WR10 prototype built initially, but the design is scaled down. To fit a standard UG-387/U flange at each detector port, the entire prototype device needed to be five inches in length, which is significantly longer than the three inches a science device would be. This increased length would significantly increase the impact of any waveguide loss on the prototype system. We measured the performance in the lab with the apparatus described in a previous article [START_REF] Bryan | A compact filter-bank waveguide spectrometer for millimeter wavelengths[END_REF] . A synthesizer generates a single microwave frequency which drives the source port of a WR5 VNA extender. The spectrum is measured by sweeping the microwave frequency, and recording the resulting signal at a diode power detector placed at the desired spectrometer output port. The measured performance of the device, shown in Figure 3 1 Modeled and calculated performance of the WR5 prototype. The measured center frequencies and spectral resolutions R of each channel agree broadly with the model. The measured and modeled optical efficiencies (OE) are shown, as well as the measured optical efficiency scaled up by a calculated value to remove the impact of waveguide loss. This value includes the physical length (5 inches) of the waveguides between the flanges and the device, as well as the R × (λ /2) effective electrical length of the resonator cavities. Since our science device will be smaller and possibly coated with superconducting niobium, waveguide loss will have a significantly smaller impact on the science device performance.

consistent with the model, and the overall reduced optical efficiency observed is consistent with the known impact of waveguide loss in the device. This will be reduced in the science device due to its smaller size, and we may sputter a superconducting niobium coating to the device to eliminate waveguide loss in this band altogether.

Individual Channel Backshorts

Like all filter bank spectrometers, the maximum possible optical efficiency of an individual WSPEC channel is approximately 50%. This is because if both the detector port, and the port that continues to the other spectrometer channels, are perfectly matched (i.e. the best possible case), the power will divide evenly between the two. Creating a filter bank with channels with closer spectral spacing than Nyquist sampling can boost the effective optical efficiency slightly higher. However, drawbacks of this approach include increased channel count and some redundancy between spectral channels.

To go beyond this 50% limit, we are considering placing a resonant backshort behind each individual channel in the spectrometer. As shown in Figure 4, placing an identical copy of the spectrometer channel a quarter-wavelength further down the main transmission line, but terminating this channel in a short instead of an absorber, increases the optical efficiency to nearly unity. Conceptually, this means that on resonance, the power that goes past the detector channel will be reflected from the short in the second channel, and couple back down the original channel to its detector. The distance between the detector channel and the backshort is selected to couple the power down the detector port, which prevents it from coupling back out of the feed horn. The HFSS geometry, and a simulation of an individual channel both with and without an individual backshort, is shown in Figure 4. Fig. 4 Drawing (left) and HFSS simulation (right) of an individual spectrometer channel with a resonant backshort. Adding a second copy of the spectrometer channel further down the main waveguide, but terminating it with a short, reflects nearly all of the in-band power onto the spectrometer channel's detector, well beyond the nominal filter bank limit of 50%.

Conclusions

Focal plane arrays of millimeter-wave spectrometers will soon enable key measurements of cosmic structure, star formation history, and studies of individual high-redshift objects. Complementing other approaches, the straightforward fabrication of these devices and the ability to make room temperature prototypes has enabled us to make significant progress developing WSPEC. In these proceedings, we presented a design and full model that includes both the spectrometer and the LEKIDs, promising measurements in a WR5 prototype that shows the concept works well at 200 GHz, and a first design for a concept that may further increase the optical efficiency beyond the nominal 50% limit for a filter bank spectrometer. The next step is building a prototype device with LEKID detectors, which will be a key step towards building a instrument to bring to a telescope for first light.

Fig. 1

 1 Fig. 1 Illustration of the spectrometer concept. Light from a feed horn (not shown) couples into the main waveguide, and different R ≈200 frequency bands are selected off by each of the five channels shown.

Fig. 2

 2 Fig. 2 Simulated passbands (right) of a full 54-channel device, including the effect of imperfect detector terminations. (The WR10 simulation is shown, but in the low-loss limit only the frequency axis needs scaling to model a device at another band.) The drawing to the left shows an illustration of a two-band WR4 / WR6 focal plane array concept using a dichroic to split the bands. A negative of each individual spectrometer is shown in grey, approximately three inches in total length. The LEKID detector cards are illustrated in green.

  Transmission through 5 inches of WR5 (VDI measured at ambient temperature)

  and Table 1, is broadly

	Measured	Model	Meas. Model Meas. Scaled Model
	Frequency	Frequency	R	R	OE	OE	OE
	156.00 GHz 157.30 GHz 186	182	15%	58%	59%
	175.23 GHz 176.65 GHz 119	120	14%	31%	35%
	175.81 GHz 177.65 GHz 156	193	22%	62%	54%
	176.59 GHz 178.33 GHz 161	207	8%	24%	26%
	203.30 GHz 205.63 GHz 220	154	10%	21%	18%
	Table

Acknowledgements This work is partially supported by the LabEx FOCUS.


R Drummond
email: ross.drummond@eng.ox.ac.uk

G Valmorbida
email: giorgio.valmorbida@l2s.centralesupelec.fr.

S R Duncan

Generalised Absolute Stability using Lyapunov Functions with Relaxed Positivity Conditions

Conditions are given for verifying stability and computing upper bounds on the induced (regional) L2 gain for systems defined by vector fields which are, along with their Jacobian, rational in the states and sector bounded nonlinearities. A class of candidate Lyapunov functions is considered that are polynomial on the states and the nonlinearities and have a polynomial scaled Lurie-Postnikov term. The main result of the paper is a set of conditions that relax the requirement on the candidate Lyapunov function from being sumof-squares with respect to the nonlinearities and the Lurie-Postnikov terms from being non-negative.

Introduction

This paper studies stability and induced L 2 gain for the class of systems described by the feedback structure of Figure 1, i.e. those with vector fields, which are, along with their Jacobian, rational on the state variables and sectorbounded nonlinearities. We refer to the stability analysis of this class of systems as generalised absolute stability. The original problem of absolute stability, formulated by Lurie and Postnikov [START_REF] Lurie | On the theory of stability of control system[END_REF] and summarised in [START_REF] Gusev | Kalman-Popov-Yakubovich lemma and the S-procedure: a historical essay[END_REF], considered the feedback interconnection of a linear system with a generic passive nonlinear operator characterised by sector bounds. The main difference between this classical problem and the one considered here is that the nonlinear system with a rational vector field replaces the linear system in the feedback block.

To analyse such systems requires techniques from both absolute stability theory and polynomial systems. A standard technique for verifying absolute stability is the search for suitable Lyapunov functions, with classical functions being the quadratic function of the circle criterion [START_REF] Khalil | Nonlinear Systems[END_REF] and the Lurie-Postnikov function of the Popov criterion [START_REF] Khalil | Nonlinear Systems[END_REF]. When the nonlinearity is also slope restricted, more advanced candidate Lyapunov functions, which are quadratic on both the states and the nonlinearity and include an integral term, can also be considered [START_REF] Valmorbida | Positivity conditions of Lyapunov functions for systems with slope restricted nonlinearities[END_REF][START_REF] Park | Stability criteria of sector-and sloperestricted Lur'e systems[END_REF][START_REF] Turner | Lyapunov functions and L 2 gain bounds for systems with slope restricted nonlinearities[END_REF]. This is the class of Lyapunov function that is generalised in this paper for the case of rational vector fields.

A generalisation of absolute stability was proposed in [START_REF] Hancock | Generalised absolute stability and sum of squares[END_REF] to extend the class of systems that could be analysed

     ẋ = f (x, φ(y), w) y = Cx z = C z x w z
φ(y) y

Figure 1. The feedback interconnection of the considered class of system. The vector field f (x, φ(y), w) and its Jacobian are assumed to be rational on x and φ(y) and affine in w.

beyond those with the rather restrictive linear feedback structure. A main feature of [START_REF] Hancock | Generalised absolute stability and sum of squares[END_REF] was the use of polynomials to bound the nonlinearities, with the results taking advantage of recent developments in polynomial optimization, namely the sum-of-squares (SOS) relaxation for positive polynomials. The generalisation presented in this paper proposes a refinement based on the introduction of Lyapunov/storage function candidates that contain polynomial terms and scaled integral Lurie-Postnikov terms that are not sum-of-squares. The advances in polynomial optimization also motivated the approximation of vector fields defined by (locally) analytic nonlinear functions by polynomials systems [START_REF] Chesi | Estimating the domain of attraction for non-polynomial systems via LMI optimizations[END_REF][START_REF] Mastellone | Nonlinear stability analysis for non-polynomial systems[END_REF]. The analysis of non-polynomial systems via polynomial optimisation methods was considered in [START_REF] Papachristodoulou | Analysis of non-polynomial systems using the sum of squares decomposition[END_REF], which used a recasting transformation on the system.

Contribution

We propose a method to study systems defined by vector fields that are, along with their Jacobian, rational functions on both the states and sector bounded nonlinearities. The restriction on the Jacobian to be rational leads to tractable numerical solutions of polynomial Lyapunov inequalities. Moreover, we provide examples of non-linearities that are not rational, but which have a rational Jacobian. The results are based on the use of Lyapunov/storage functions that are polynomial on the states and the sector bounded nonlinearities and contain a scaled Lurie-Postnikov term. The conditions for the positivity of the Lyapunov functions extends the results for quadratic functions of [START_REF] Valmorbida | Positivity conditions of Lyapunov functions for systems with slope restricted nonlinearities[END_REF], namely, we do not impose the Lurie-Postnikov term nor the polynomial terms on the sector nonlinearities to be strictly positive. The prescribed analysis also encompasses regional stability and performance, with inclusion conditions on the level sets of the Lyapunov function being used to obtain local properties. More importantly, the conditions presented are expressed in terms of the verification of polynomial inequalities, which are cast as semi-definite programmes, convex optimisation problems that can be efficiently solved for low dimensional systems using software packages such as YALMIP or SOSTOOLS.

A preliminary version of the material in this paper can be found in [START_REF] Drummond | Lyapunov analysis of nonlinear systems with rational vector field and Jacobian[END_REF], where only a global analysis was considered and a simpler version of the Lyapunov function than the one studied here was proposed. Note that since the vector field is defined by rational expressions and not by a linear system, it is not possible to formulate the frequency dependent interpretation of the stability conditions, such as phase equivalent Zames-Falb multipliers [START_REF] Carrasco | Zames-Falb multipliers for absolute stability: From O'Shea's contribution to convex searches[END_REF], which are generally considered to be the least conservative method for the global analysis of absolute stability [START_REF] Carrasco | Zames-Falb multipliers for absolute stability: From O'Shea's contribution to convex searches[END_REF].

The paper is structured as follows. Section I sets up the system to be studied. Section II introduces a set of inequalities for verifying stability and computing upper bounds on the induced L 2 gain from input to output of the system in some local/global domain. Section III details a numerical formulation for the inequalities such that they can be verified using convex programming. Numerical examples illustrate the proposed method in Section IV. The first numerical example shows the existence of a common Lyapunov function for a benchmark saturated system while the second example considers the stability of a nonlinear vehicle model.

Notation: The i th row of matrix C ∈ R p×n is denoted C i and the i th element of vector y is denoted y i . A polynomial R[x] is a function that can be expressed as deg(p,θ) stands for the degree of p on variable θ. The set of non-negative reals is denoted R ≥0 . The set of real symmetric matrices of dimension n is denoted S n and the subset of diagonal matrices is denoted D n . The sets of positive semidefinite symmetric and diagonal matrices are respectively S n ≥0 and D n ≥0 . The set of real positive definite matrices are S n >0 and D n >0 . The vector of dimension n containing zeros is denoted

p(x) = p 1 x + • • • + p m x m . Polynomials of real coefficients of dimension n in variable θ are denoted R n [θ]. We use R[θ] for n = 1. Polynomials on two variables θ 1 and θ 2 are denoted R n [(θ 1 , θ 2 )]. The positive polynomials of di- mension n are denoted R n >0 [(θ 1 , θ 2 )] and the non-negative polynomials are R n ≥0 [(θ 1 , θ 2 )]. The positive sum-of-squares polynomials of dimension n are denoted Σ n [(θ 1 , θ 2 )]. The diagonal sum-of-squares polynomials are Σ diag n [(θ 1 , θ 2 )]. For p ∈ R[θ],
0 n . The L 2 norm of a signal is denoted z 2 = ∞ 0 z T (θ)z(θ)dθ. We define the ρ sublevel set of V by E(V, ρ) = {x ∈ R n | V (x, φ) ≤ ρ} for ρ ≥ 0. The time derivative of a function V (x, φ) : X 0 → R ≥0 is denoted by the inner product V = dV (x,φ) dx ẋ+ dV dφ φ = V (x, φ) = ∇ x V ∇ φ V , ẋ φ .
In a slight abuse of notation, we may drop the argument of functions, most notably with φ(y(x)) = φ and φ i (y i (x)) = φ i .

I. Problem Statement Consider the dynamical system ẋ = D 1 (x, φ) -1 N 1 (x, φ) + B w w (1a) y = Cx (1b) z = C z x (1c) with x ∈ R n , y ∈ Y ⊆ R m , φ(y) : Y → R m , φ ∈ C 1 , output z ∈ R nz , disturbance w ∈ R nw , N 1 ∈ R n [(x, φ)], D 1 ∈ R >0 [(x, φ)].
The nonlinearity φ(y(x)) is defined on Y, which contains the origin and a domain Y 0 satisfying Y 0 ⊆ Y ⊆ R m . Define the set X 0 as X 0 = {x|y(x) ∈ Y 0 }. This system can be interpreted by the feedback form of Figure 1 with f (x, φ(y), w) defined by the RHS of (1a).

We will restrict the considered class of nonlinearities to those which are memoryless, static, decentralised

φ(y) = φ 1 (y 1 ), φ 2 (y 2 ), . . . , φ m (y m ) T , (2a)
sector bounded with the nonlinearities φ i being signpreserving operators on y i satisfying

φ i (y i ) y i ∈ [0, δ i ] ∀y ∈ Y 0 ⊆ Y (2b) and dφ i (y i ) dy i = D 2i (y i , φ i ) -1 N 2i (y i , φ i ), (2c)
where

N 2i ∈ R[(y i , φ i )] and D 2i ∈ R >0 [(y i , φ i )]. Defining D2 (y, φ) ∈ R >0 [(y, φ)] and N2 (y, φ) = diag( N21 , . . . , N2m ) as D2 (y, φ) = i=1, ..., m D 2i (y i , φ i ) (3) N2i (y, φ) = N 2i (y i , φ i ) j=1, ..., m i =j D 2j (y j , φ j ) (4)
such that

D -1 2i (y i , φ i )N 2i (y i , φ i ) = D-1 2 (y, φ) N2i (x, φ) (5)
gives

dφ(y) dy = D2 (y, φ) -1 N2 (y, φ). (6)
The reason for introducing D2 (y, φ) and N2 (y, φ) is to allow the denominator terms to be cancelled out in the stability conditions. The sector boundedness of the nonlinearity can be used to obtain a polynomial inequality, known as the sector inequality, that relates y and φ(y).

Fact 1: If the sector condition of (2b) holds, then

s 1 (y, φ, T (x, φ)) := φ T T (x, φ) ∆y -φ ≥ 0 ∀y ∈ Y 0 (7)
where

T = diag(τ 1 , τ 2 , . . . , τ m ) with τ i ∈ R ≥0 [(x, φ)] for i = 1, 2, . . . , m and ∆ = diag(δ 1 , . . . , δ m ).
The sector inequalities relating φ(y) and its argument are used to verify the positivity of polynomials in y and φ(y) as in the fact below.

p 1 (y, φ)-s 1 (y, φ, T (x, φ)) ≥ 0 ∀x ∈ X 0 (8)
holds then p 1 (y, φ) ≥ 0 ∀x ∈ X 0 . This paper proposes solutions to the following problems: Problem 1 (Stability): Provide conditions that guarantee local stability of the origin of system (1).

Problem 2 (L 2 gain): Obtain (local) upper bounds for the induced L 2 gain from w to z of (1).

II. Main Results

This section gives a set of conditions in terms of inequalities to solve Problems 1 and 2. The analysis is based upon Lyapunov function candidates of the form

V (x, φ) = V 0 (x, φ) + m i=1 yi 0 λ i (σ) φ i (σ) dσ (9)
where V 0 ∈ R[(x, φ)], V 0 (0, 0) = 0, and

λ i ∈ R[σ], i = 1, . . . , m.
The integral term in ( 9) is referred to as the Lurie-Postnikov term. Using ( 7), a polynomial lower bound for V (x, φ) is given by the following lemma:

Lemma 1: Let V 0 ∈ R[(x, φ)] and λ i ∈ R[σ] for i = 1, . . . , m defining (9) be given and let X 0 ⊆ R n , 0 ∈ X 0 with φ satisfying (2). If there exists λi ∈ R ≥0 [σ] for i = 1, . . . , m, such that λ i (σ) ≥ -λi (σ), i = 1, . . . , m, ∀x ∈ X 0 (10a) V (x, φ) = V 0 (x, φ) (10b)
- m i=1 yi(x) 0 δ i σ λi (σ) dσ > 0 ∀x ∈ X 0 ,
where δ i is the sector bound, then V (x, φ) > 0 ∀x ∈ X 0 \{0}. Proof: We use [START_REF] Papachristodoulou | Analysis of non-polynomial systems using the sum of squares decomposition[END_REF] and the sector bounds (2b) to obtain that V (x, φ) is lower bounded by V (x) as follows

V (x, φ) = V 0 (x, φ) + m i=1 yi(x) 0 λ i (σ)φ i (σ)dσ ≥ V 0 (x, φ) - m i=1 yi(x) 0 λi (σ)φ i (σ)dσ = V 0 (x, φ) - m i=1 yi(x) 0 λi (σ)φ i (σ)dσ - m i=1 yi(x) 0 δ i σ λi (σ)dσ + m i=1 yi(x) 0 δ i σ λi (σ)dσ = V 0 (x, φ) - m i=1 yi(x) 0 δ i σ λi (σ)dσ + m i=1 yi(x) 0 λi (σ)(δ i σ -φ i (σ))dσ ≥ V (x, φ). (11)
The following gives a set of inequalities that use the lower bound of Lemma 1 to verify solutions to Problems 1 and 2. Note that the term

p l (x) := m i=1 yi(x) 0 δ i σ λi (σ) dσ (12) in (10b) is a polynomial P l ∈ R[x].
Theorem 1: Suppose there exists polynomials

V 0 ∈ R[(x, φ)], T j = diag(τ j1 , τ j2 , . . . , τ jm ) j = 1, 2 with τ j,• ∈ R ≥0 [(x, φ)], λi ∈ R ≥0 [σ], λ i ∈ R[σ], i = 1, 2, . . . , m
and a scalar η ∈ R >0 such that (10a) and

V (x, φ) -s 1 (x, φ, T 1 (x, φ)) ≥ 0 ∀x ∈ X 0 , (13a)
- ∇ x V ∇ φ V , ẋ φ -Ω(z, w) (13b)
-

s 1 (x, φ, T 2 (x, φ)) > 0 ∀x ∈ X 0 , E(V, ρ) ⊆ X 0 (13c)
hold with a) Ω(z, w) = 0 and w = 0 ; b) Ω(z, w) = -w T w + η -1 z T z; then a) The origin of ( 1) is locally asymptotically stable with domain of attraction E(V, ρ). b) The induced L 2 gain from w to z is upper-bounded by √ η, that is, z 2 w 2 ≤ √ η for all solutions of (1) satisfying

x(0) = 0 and w 2 ≤ ρ 1 2 . Proof: From Fact 2, if (13a) holds, then V (x) > 0 ∀ x ∈ X 0 . Positivity of V (x, φ) > 0 then follows from Lemma 1. Denote by V (x) the time derivative of the Lyapunov function V (x, φ) along the trajectories of ( 1)

such that V (x, φ) = ∇ x V ∇ φ V , ẋ φ . Invoking Fact 2, if (13b) holds, then -V (x, φ) -Ω(z, w) > 0 ∀x ∈ X 0 . (14)
With a) Ω(z, w) = 0, then (14) implies V (x, φ) < 0 ∀x ∈ X 0 . The set inclusion condition (13c) enforces the containment of the sublevel sets of the Lyapunov function E(V, ρ) in the region where the local sector holds. The satisfaction of (13a) and (13b) implies positive invariance and contraction of the sublevel sets E(V, ρ) which define an inner estimate of the region of attraction of the system. b) Ω(z, w) = -w T w + η -1 z T z, then integrating [START_REF] Valmorbida | Region of attraction estimation using invariant sets and rational Lyapunov functions[END_REF] from 0 to t * with x(0) = 0 gives

t * 0 w T (τ )w(τ )dτ > η -1 t * 0 z T (τ )z(τ )dτ + V (x(t * )). For w 2 ≤ ρ 1 2 , we obtain ρ > η -1 t * 0 z T (τ )z(τ )dτ + V (x(t * )) ≥ V (x(t * )), that is x(t * ) ∈ E(V, ρ), thus from (13c) it follows that x(t * ) ∈ X 0 for t * ∈ [0, ∞).
Therefore, with w ≤ ρ 1 2 , (13b) and (13c) hold true and

t * 0 w(τ ) T w(τ )dτ ≥ η -1 t * 0 z(τ ) T z(τ )dτ that is z 2 w 2 ≤ √ η.
Remark 1: Note that Theorem 1 does not require strict positivity of the polynomials λ i , nor of the terms of V 0 related to the variables φ. Positivity of the nonlinear terms of V 0 is relaxed by including the sector information into (13a) while the positivity of λ i is relaxed by using the polynomial lower bound of the Lurie-Postnikov term given in Lemma 1.

III. Numerical Formulation

In this section, we propose numerical relaxations to solve the polynomial inequality constraints of Theorem 1. These relaxations use the SOS formulation for positive polynomials, which is briefly explained in the appendix. See [START_REF] Blekherman | Semidefinite optimization and convex algebraic geometry[END_REF] for a comprehensive analysis.

The set containment condition in (13c) can be cast as a polynomial inequality in the case when the Lurie-Postnikov term is itself a polynomial. Otherwise, we verify (13c) by relying on (10b) and use a polynomial outer approximation of E(V, ρ). This outer approximation is obtained by noting from (10b) that V (x, φ) ≤ V (x, φ), thus E(V, ρ) ⊆ E(V , ρ).

To formulate [START_REF] Blekherman | Semidefinite optimization and convex algebraic geometry[END_REF] as polynomial inequalities, in this section we assume X 0 to be a semi-algebraic set, that is

∃ g i (x) ≥ 0, i = 1, . . . , n p such that X 0 = {x ∈ R n |g i (x) ≥ 0, i = 1, . . . , n p }.
In the proposition below, we denote g(x) = [g 1 (x), . . . , g np (x)].

The following proposition gives a convex formulation of Theorem 1, expressed as an SDP.

Proposition 1:

If there exists V 0 ∈ R[(x, φ)], V 0 (0, 0) = 0, T 1 , T2 ∈ Σ diag m [(x, φ)], λ i (σ) ∈ R[σ], λi (σ) ∈ Σ[σ], satisfying λ i (σ) + λi (σ) ∈ Σ[σ] for i = 1, . . . , m, q j ∈ Σ np [(x, φ)] for j = 1, . . . , 3, η ∈ R >0 , > 0, V (x, φ) -s 1 (y(x), φ, T 1 ) -εx T x -q T 1 (x, φ)g(x, φ) ∈ Σ[(x, φ)], (15a)
- ∇ x V ∇ φ V , D 1 (x, φ)) -1 N 1 (x, φ) + B w w D2 (φ, y) -1 N2 (φ, y) ẏ(x) (15b)
-Ω(w, z)

• D 1 (x, φ) D2 (x, φ) -s 1 (y, φ, T 2 ) -q T 2 (x, φ)g(x, φ) -εx T x ∈ Σ[(x, φ, w)], ρ -V (x, φ) -q T 3 (x, φ)g(x, φ) ∈ Σ[x], (15c)
then a) Ω(w, z) = 0 and w = 0, the origin of ( 1) is locally asymptotically stable with domain of attraction E(V, ρ). b) Ω(z, w) = w T w -η -1 z T (x)z(x) an upper-bound for the induced L 2 gain of ( 1) is given by √ η for all solutions of (1) satisfying x(0) = 0, x(t) ∈ X 0 and w ≤ ρ

Proof:

We show that (15) implies [START_REF] Blekherman | Semidefinite optimization and convex algebraic geometry[END_REF]. From (15a), it follows that V (x) -s 1 (y(x), φ, T 1 ) ≥ εx T x + q T 1 (x, φ)g(x, φ). Since g j (x, φ) ≥ 0 ∀x ∈ X 0 , we obtain (13a).

Similarly, (15b) implies

∇ x V ∇ φ V , D 1 (x, φ) -1 N 1 (x, φ) + B w w D2 (φ, y) -1 N2 (φ, y) ẏ -Ω(w, z) • D 1 (φ, y) D2 (φ, y) ≥ 0 ∀x ∈ X 0 . (16)
Dividing the above by the positive polynomial (D 1 (x, φ) D2 (x, φ)) we obtain [START_REF] Valmorbida | Region of attraction estimation using invariant sets and rational Lyapunov functions[END_REF], which is Fact 2 applied to (13b). Finally, (15c) implies ρ -V (x, φ) ≥ 0 ∀x ∈ X 0 , thus E(V , ρ) ⊆ X 0 and since E(V, ρ) ⊆ E(V , ρ), (13c) holds.

Remark 2: For quadratic formulations V (x, φ) = x T P x, a convex optimization problem that enlarges the estimate of the region of attraction is obtained using the function log( det(P )). For the enlargement of the estimate of E(V, ρ) used here, one can rely on a sequence of solutions to SDPs optimising the ERA as proposed in [START_REF] Valmorbida | Region of attraction estimation using invariant sets and rational Lyapunov functions[END_REF].

IV. Examples

A function satisfying (2b) and (2c) globally (with Y = R) is the tan -1 (y) function. Indeed, in this case, (2b) holds with δ = 1 and the slope constraint (2c) holds since

d tan -1 (y) dy = 1 1 + y 2 ∈ [0, 1). ( 17
)
The function tan -1 (y) can be regarded as a smooth approximation of the saturation function and also appears in applications, as detailed in the example below.

A. Numerical Example 1: Smooth Saturation

Consider the Lurie-Postnikov system ẋ1 ẋ2 = -0.05 1 -10 -0.5

x 1 x 2 + 0 1 tan -1 (y), (18a)
y = 9.9 0.495 x 1 x 2 . ( 18b
)
The above system is adapted from [15, Example 4.3.1] with the difference being that smooth saturation is considered here, with the non-differentiable saturation function approximated by tan -1 (y). In [START_REF] Zaccarian | Modern anti-windup synthesis: control augmentation for actuator saturation[END_REF]Example 4.3.1], the authors show that there does not exist a common quadratic function for both the open-loop system and the linear closed loop system (i.e. the system obtained by replacing tan -1 (y) by y in (18a)). However, using V (x, φ) as in [START_REF] Mastellone | Nonlinear stability analysis for non-polynomial systems[END_REF] with deg(V 0 (x, φ), x) = 4 and the degrees of the other polynomials of Proposition 1 being two, we can prove the global stability of the non-linear closed-loop system. Suitable level sets of V (x, φ) depict the decrease of the function along the trajectories of the system in Figure 2.

B. Numerical Example 2: Vehicle Stability

The second example considers the application of the presented results to a physical system and involves the estimation of the maximal region of attraction of a vehicle model. Consider the following vehicle dynamics model from [START_REF] Masouleh | Non-linear vehicle domain of attraction analysis using sum-of-squares programming[END_REF] m v = -mu 0 ω + F f y (y 1 ) + F ry (y 2 ) (19a)

I ω = -aF f y (y 1 ) -bF ry (y 2 ) (19b) The nonlinearity of the model enters via the tyre forces F f y (y 1 ) and F ry (y 2 ) given in [START_REF] Masouleh | Non-linear vehicle domain of attraction analysis using sum-of-squares programming[END_REF] by

y 1 = v u 0 + aω u 0 (19c) y 2 = v u 0 - bω u 0 (19d)
F f y = -2C f (φ 1 (y 1 ) -K f φ 1 (y 1 ) 3 ) (20a) F ry = -2C r (φ 2 (y 2 ) -K r φ 2 (y 2 ) 3 ) (20b)
where C r = C f = 57300 N rad -1 are the cornering stiffnesses and K f = K r = 4.87 are tyre force parameters. Compared to the "Magic Formula Tire Models" of [17,

Chapter 4], this definition of the tyre forces assumes that the tyre slips do not become large. Under the assumption of small steering angle deflection, the front and rear tyre lateral slip angles are

φ 1 (y 1 ) = tan -1 (y 1 ) , (21a)
φ 2 (y 2 ) = tan -1 (y 2 ) . (21b)
This vehicle model can be expressed as [START_REF] Lurie | On the theory of stability of control system[END_REF] with

D 1 = 1, B w = 0, C z = 0, D2 = (1 + y 2 1 )(1 + y 2 2 ), C = 1 u0 a u0 1 u0 -b u0 , (22a)
N 1 = -u 0 ω 0 (22b) + - 2C f m 2C f K f m -2Cr m 2CrKr m 2aC f I - 2aC f K f I -2bCr I 2bCrKr I     φ 1 (y 1 ) φ 1 (y 1 ) 3 φ 2 (y 2 ) φ 2 (y 2 ) 3     , N2 = (1 + y 2 2 ) 0 0 (1 + y 2 1 )
.

The region of attraction for this model obtained using Proposition 1 is shown in Figure 3. We have used deg(V 0 (x, φ), x) = 8 and with the other polynomial variables of Proposition 1 being of degree 2. Including the nonlinearities of the model directly within the stability analysis of the system in this manner will give a more representative picture of the stability profile of the real system compared to using a linear analysis, as is commonly used [START_REF] Pacejka | Tire and vehicle dynamics[END_REF]Chapter 11].

V. Conclusion

This paper has considered the generalised stability problem with conditions being derived to verify stability and obtain minimal L 2 gain bounds from disturbance to output. The class of nonlinear systems considered is, along with its Jacobian, a rational function of the states and sector bounded nonlinearities. The conditions were derived using a Lyapunov (storage) function that is polynomial with respect to the states and the nonlinearity and has a polynomially scaled integral term. The main result of the paper relaxes strict positivity of the terms of this Lyapunov function so as to reduce the conservatism of the conditions. As a consequence, the Lyapunov functions are not required to be sum-of-squares with respect to the nonlinearity and the scaling terms are also not required to be strictly positive. A convex formulation for the conditions was given such that they could be efficiently verified, as shown in two numerical examples.

A system of practical interest in the form (1) is the finite dimensional approximation of the electrochemical model of a supercapacitor [START_REF] Drummond | Loworder mathematical modelling of electric double layer supercapacitors using spectral methods[END_REF]. In such a system, the (locally) sector bounded nonlinearity ln(1 + ax) appears. We are currently extending the results of this paper to the problem of computing the L 2 gain of the supercapacitor to estimate its resistance. Due to the relatively high dimension (20 states) of the discretised supercapacitor model, a numerical solution of the optimization problem in Proposition 1 currently cannot be obtained. Appendix Definition 1: A multivariate polynomial p(θ) is said to be a sum-of-squares polynomial if there exists a finite set of polynomials h i ∈ R[θ] satisfying

p(θ) = n i=1 h 2 i (θ). ( 23
)
The search for a sum-of-squares decomposition of a polynomial is affine in the coefficients of h i and as such can be expressed as a semi-definite programme (SDP) [START_REF] Blekherman | Semidefinite optimization and convex algebraic geometry[END_REF]. Proposition 2: A polynomial p(θ) ∈ R[θ] of degree 2d is an sum-of-squares polynomial if and only if there exists Q ∈ S ≥0 and a vector of monomials Z(θ) ∈ R M [θ] such that p(θ) = Z(θ) T QZ(θ)

then p(θ) is sum-of-squares.

Fact 2 :

 2 Consider a polynomial p 1 ∈ R[(y, φ)]. If the inequality

2 Figure 2 .w 3 .

 223 Figure 2. Convergence of the state trajectories of (18) (blue). Level sets of the computed function V (x, φ) (black) and lines corresponding to the set {x ∈ R 2 : | tan -1 (y(x))| = 0.8} (red) are also shown.


A S Yigit

Associate Professor. R A G A Ulsoy

Scott

Dynamics of a Radially Rotating Beam With Impact, Part 2

Introduction

A theoretical and computational model using the momen tum balance (coefficient of restitution) method was presented in Part 1 (Yigit et al., 1988a). In this paper experimental and numerical studies on this subject are presented. Experimental results are compared with simulation results in order to assess the validity of the model.

The objective of the experimental work is to investigate the validity of the analytical impact model presented in Part 1, and to determine its applicability for flexible systems. This task is undertaken by considering a radially rotating flexible beam, with a rigid stem, impacting a small semi-cylindrical hard surface. Three types of measurements are made: torque input, rigid body motion, and elastic motion of the flexible beam. The torque input given to the shaft by a de motor was measured through a current measurement. The angular veloc ity of the rigid stem was measured by using an optical in cremental encoder. Electrical resistance strain gages at four different locations along the beam were used to obtain the elastic motion. The measured angular velocity of the rigid shaft and the strains at four different locations were compared with simulations using a torque input approximating the measured torque. Also a set of experiments were performed with no torque input, and the natural frequencies of the non rotating and rotating beam were obtained. High speed video : ecordings were employed to show the existence of multiple Impacts.

Experimental results are compared with simulation results in order to assess the validity of the model. Sensitivity studies are given to determine the applicability of the model for a wider range of system parameters than covered by the experiments.

In Section 2 the experimental design used is presented. The calibration procedure for the instrumentation is explained in Section 3. Experimental and simulation results are given and compared in Section 4. Section 5 is devoted to sensitivity studies. Finally some concluding remarks are given in Section 6.

Experimental Design

The schematic diagram of the experimental setup is shown in Fig. 1. The individual components in the setup are described below.

• Flexible Beam: A slender aluminum bar with the follow ing dimensions and properties was selected: length = .530m, diameter= .0063m, mass density= 2700 kg!m 3 , Young's modulus, E=6.895 1010 Nlm 2 • The beam is can tilevered on a rigid shaft, and was designed so that the natural frequencies are small and at the same time the deflections under the impact are within the elastic limit. Capture of the third mode contribution was to be within the range of the measurements. The specifications given above satisfy these requirements.

• DC Motor: The de motor used for driving the system was an Electro-Craft-Model M-1030 which is an in termediate performance de moving coil servomotor. The damping factor is 4 w-4 Nms and the minimum Motor Tach resonant frequency is 4 kHz. The magnitude of the torque which is delivered by the de motor was controlled using an adjustable voltage and current power supply. Fig. 1 The sketch of the experimental setup is 50kHz. Accuracy (pulse to any other pulse) is ± 1.0 minutes of arc. The encoder produces pulses from two different channels for clockwise and counterclockwise rotations.

• Frequency-to-Voltage Converter: An Analog Devices Model AD650 Frequency-to-Voltage Converter is used to convert the encoder pulse frequency into an analog signal.

The unit is capable of doing frequency-to-voltage and voltage-to-frequency conversions for frequencies up to 1MHz. The typical nonlinearity at 100 kHz is 0.005 per cent, and it can accept input pulses of duration greater than lOOns. Two separate converters had to be used for the two different channels from the encoder because the frequency-to-voltage converter output is always positive. Then these two outputs were summed digitally in the computer.

• Current Transducer: The torque delivered to the rigid shaft by the de motor was measured using a Model CT 1 OOL TT current transducer by Ohio Semitronics which measures the current using the Hall effect (Walden). This transducer provides more accurate and safer measure ments than conventional shunt measurements. The band width is 0 to 5 kHz., linearity is 0.5 percent full scale, and its response time is less than 50 psec. The de signal condi tioner used with the unit was Model CTA101 Ohio Semitronics. Its frequency range is also 0 to 5kHz., and 40 psec is required to reach 90 percent of the steady state value.

• Straingages and Straingage Amplifier: Strains at four locations along the beam are measured using Micro Measurements EA 13-120LZ-120 straingages. The gage factor is 2.095 ±0.5 percent. These straingages have an ex cellent frequency response up to 50kHz, and they have been tested for dynamic measurements [START_REF] Hetenyi | Handbook of Experimental Stress Analysis[END_REF]. A half bridge circuit arrangement was used employing two straingages at a given location on the upper and lower sur faces of the beam. Self-temperature compensation is pro vided. A Model 1020 Micro-Measurements Straingage Amplifier is used for the straingage instrumentation.

Straingages were preferred to accelerometers because of the presence of the impact. An accelerometer which can withstand the impact is not generally capable of measure ing the smaller magnitudes of the vibration before and after the impact. " Analog Filter:The output signals from the frequency-to voltage converters contained high frequency noise in- traduced by the converter itself. In order to eliminate this noise, an analog filter was used with a cutoff frequency of 3kHz. The maximum frequency in the velocity measure ment corresponding to the third mode was about 200 Hz, so this filtering has no significant effect on the measure ment itself. The analog filter was a Khron-Hite Digitally Tuned 3320 Series Low Pass Filter which had an ad justable cut off frequency limits from 0.001 Hz to 99.9 kHz. " Instrumentation Tape Recorder: Data was first recorded on an analog tape recorder to avoid the sampling time restrictions of digitizing. The instrumentation tape recorder was a Model Store 7DS Racal Tape Recorder.

This tape recorder has 7 data channels which have the capability of recording signals of up to 40kHz in frequen cy, depending on the recording speed and mode. In our case, FM recording was used and the speed at which the signals were recorded was 60.0 in/sec, thus achieving a bandwidth of 0 to 20 kHz. Overall system linearity is 0.3 percent, signal to noise ratio is 52 db and the harmonic distortion is less than 1.2 percent at the maximum modulation level. " Computer: The collected data were digitized on a DEC PDP-11123 computer by playing the tape recorder back at a slower speed, thus allowing a much longer actual sam pling rate. The effective sampling rate was 1kHz, cor responding to an actual sampling rate of 62 Hz. Dimensions and structural properties are presented in Table 1. The damping in the rigid shaft joint was calculated as follows: The rigid shaft (without the flexible beam mounted) was rotated with a constant angular velocity, then the motor stopped. The decay of the angular velocity was recorded on a digital storage oscilloscope. The time constant of the system was obtained neglecting the compliance in the motor shaft. The angular velocity of the shaft can be written as:

0 =Ooe-U!r)
where the time constant r can be expressed as:

r=I!c (1) (2)
where I is the inertia of the shaft and cis the viscous damping coefficent. Then knowing the timer for the angular velocity to drop to a value which is 33 percent of the initial value, the viscous damping coefficient can be obtained from equation [START_REF] Hetenyi | Handbook of Experimental Stress Analysis[END_REF].

3 Calibration

Straingage Calibration:

The strain measurement was calibrated statically. The beam was mounted as a cantilever beam and loaded statically by a series of known weights. Measured strains were compared with the calculated strains.

Torque Measurement Calibration: A balance scale was used to calibrate the system. The de motor was turned on with the tip of the beam touching the pan of the scale. Then measuring the force, the applied torque was calculated and compared with the measurements of the current transducer. Very short response time and very low inertia of the com ponents in the system justify the static calibration procedure.

Velocity Measurement Calibration: The optical encoder produces 16,384 pulses per revolution, which means that a fre quency of 16,384Hz would correspond to an angular velocity of 21r[rad/s]. Thus the only calibration necessary was that associated with the frequency-to-voltage converters. This was done using a function generator. The calibration curves for all measurements, and the results of linear regression analyses for these curves are given in [START_REF] Yigit | Dynamics of a Radially Rotating Beam with Im pact: Implications for Robotics[END_REF].

Experimental Results and Comparison

A set of experiments were performed and recorded using the KODAK Ekta-Pro high speed video recording system. Two separate cameras were used to focus on both the entire beam and on the impact region. The camera focused on the entire beam was intended to capture the mode shapes of the beam during its vibration, and the impact point was studied to record multiple impacts clearly. In the first set of recordings only the impact region was studied. As a second trial the screen was divided into two to show both the entire beam and the impact region. In these experiments the recording speed was 1000 frames/second. The mode shapes of the beam can be seen by playing the recordings back at a slower speed (30 frames/second). In order to see more clearly the multiple im pacts the screen was divided into three and the images from the same camera were recorded on one frame thus achieving a recording speed of 3000 frames/second. Multiple impacts were clearly observed in these recordings [START_REF] Yigit | Dynamics of a Radially Rotating Beam with Im pact: Implications for Robotics[END_REF].

4.1 Results Without Torque Input. Some preliminary simulations were run to test accuracy. For the experiment the rigid shaft is held fixed so that a cantilever configuration is ob tained. The tip of the cantilever is given a vertical impulse and the strains at four different locations along the beam are measured. The structural damping present in the flexible beam is seen to be negligible, and it is assumed to be zero for the simulations later. In order to determine the natural frequen cies a discrete Fourier transform (DFT) is applied to the time response, and the first natural frequency was obtained to be 17 Hz. Almost identical results were obtained from the other three locations. Fig. 2(a) shows the free response of the can tilever beam at location 3. This result was obtained using seven modes in the simulation. The fundamental frequency is 18 Hz, and this agrees well with the exact value. Experiments give a fundamental frequency of 17 Hz [Fig. 2(b)]. The small dif ference is probably due to the fact that the ideal cantilever boundary conditions were not achieved in the experimental setup.

Secondly, the beam was considered to fall freely from an in itial angle 0 0 , and impact on an aluminum surface. In the simulation, the coefficient of restitution is specified as a func tion of impact velocity. For aluminum, for impact velocities up to 0.2[m/s] experimental data is available [START_REF] Goldsmith | Impact[END_REF]. For impact velocities higher than 0.2[m/s] the value of 0. 000 1 0 ,---,--,--.---r--.--. .,--.,. --,----,---- e was determined from the following analytical result [START_REF] Johnson | Contact Mechanics[END_REF]:

e= 0.501v -< 1 14> (3)
Simulation and experimental results on angular velocities and strain histories are shown in Figs. 3 and4. Peak strains agree well as can be seen by comparing Figs. [START_REF] Mason | Impact on Beams[END_REF] (a) and 4(b). Due to the Coulomb friction at the joint (not modeled in the simula tion) there is a time delay in the experimental results. This ef fect can be observed by comparing Figs. 3(a) and 3(b). The ef fect of Coulomb friction is seen only at the beginning of the motion where the angular velocity obtained from the simula tion shows a much higher negative peak at the time of impact than the one observed in the experiment. This difference is probably due to the fact that the angular acceleration for this region is beyond the limit of the incremental optical encoder which was used to measure the angular velocity. The average acceleration at the time of impact as calculated from the simulation is 110,000[rad/s 2 ], whereas the maximum accelera tion for the encoder is 100,000[rad/s 2 ]. The negative peaks for the subsequent impacts are captured accurately since they fall ¥•ithin the acceleration limit of the encoder. Discrete Fourier transforms show that the dominant fre quency predicted by the simulation is about 67 Hz, compared to 66 Hz in the experiment [note that the frequency of the rotating beam is higher than the cantilever frequency of 17 Hz, due to the coupling between the elastic and the rigid body mo tions, as discussed in Yigit et a!. (1988b)].

It is interesting to note that both the strain and angular velocity results show the existence of multiple impacts. This is best seen in Fig. 5, which is a magnification of the first con tact/rebound phase. This phenomenon of multiple impacts has also been reported by other researchers (e.g., Goldsmith, The more challenging case of a torque is now examined. The beam impacts an aluminum surface, and the rigid shaft is given a rectangular torque pulse input which approximates the torque input in the corresponding experiments (see Fig. 6). Simulation and ex perimental results on angular velocities and strain histories are shown in Figs. 7 and8. Some differences between the experi ment and simulation are seen in the angular velocity results [(Figs. 7(a) and 7(b)]. The time delay in the experimental results is due to the coulomb friction which is effective during the initial period of the motion when the angular velocity is low. The negative peaks is the simulation differ with those in the experiment because of the limitation on the maximum angular acceleration of the encoder. Some other possible sources for the differences include: (1) using a rectangular pro file for the torque input. As can be seen from Fig. 6 the actual torque input does not drop to zero instantaneously. This con tributes to t he difference in the slopes of the angular velocity just before the impact, (2) neglecting the coulomb friction in the joint: the coulomb friction becomes effective at low angular velocities. It tends to delay the time of the first im pact, and it causes the second impact after bouncing back to occur earlier than predicted by the simulation. Excellent agree ment is found for the peak strains [compare Figs. . .,-3 can conclude that the difference in the dynamic behavior is very minor. It can be concluded that as long as the value of the coefficient of restitution is appropriate for the first impact (which is responsible for the bulk of the energy loss), one can use the same value for the subsequent impacts.

•o 0 � � bD c < � -, 1.. '"'>-
A set of simulations and experiments were also performed for a steel and brass impact surface. In the simulations the value of e for an aluminum-steel and aluminum-brass material pair was obtained using equation (30) of Yigit et a!. (1988a). Torque inputs, which were obtained from the corresponding experiment for the steel and brass surface were used. Again, excellent agreement was found between the simulation and ex perimental results (see [START_REF] Yigit | Dynamics of a Radially Rotating Beam with Im pact: Implications for Robotics[END_REF].

From these comparisons one can conclude that the momen tum balance (coefficient of restitution) method can be used to predict the post impact dynamic behavior of similar systems with good accuracy. It must be noted, however, that the ex periments performed in this study represent a particular con figuration. In order to generalize the results some sensitivity studies were performed. These are presented in the following section.

Sensitivity Studies of the Momentum Balance

Method

In Part 1 (Yigit et a!., 1988a) an analytical investigation was made on a simple one degree of freedom system, and it was found that for most engineering systems the effect of flexibil ity on the coefficient of restitution would be negligible. The good agreement between theory and experiment reported above also supports this conclusion. To give these findings greater generality, some sensitivity studies will now be presented.

There are many possible measures of sensitivity. Here three different sensitivity measures are introduced as follows:

where w(l,t) max are peak elastic deflections calculated for dif ferent values of e. W (l,t) m ax is the arithmetic mean of the w(l,t) max values.

2 Tip displacement sensitivity, s2: Another possible measure for sensitivity is the total displacement w P due to both elastic and rigid body motion. Because of the rigid body mo tion, w P does not have a maximum value, and so replacing w by w P in equation ( 4) is not appropriate. A sensitivity measure s 2 (which is a function of time is introduced by (5)

where w P (l,t) is the displacement corresponding to an ar bitrary value of e, and w P (l,t)e=o.6 is the tip displacement for e = 0.6. Here the value 0.6 is selected arbitrarily as an in termediate value of e in the range of practical values. t1 denotes a time after impact. In the sensitivity plots presented the t1 is taken as the final time in simulations.

3 Global tip displacement sensitivity, s3: As a further measure the time integral of the absolute value of the sensitiv ity s 2 was used. This sensitivity measure is defined as:

r 1 1 Js 2 Jdt Jto s 3 (6) t1 -t0
where t0 is the time of the impact, and t1 is a given time after impact. Values of s1 as functions of e are given in Fig. 11 for two different flexural rigidity values. In both cases elastic peak deflections vary only by about 5 percent for a wide range of coefficient of restitution values in the practical range of 0.4 to 0.8. Figure 12 shows s1 for two different inertia ratios {3 = (I R I I 8). Again in both cases elastic peak deflections change only of the order of 5 percent for e in the range of 0.4 to 0.8.

Results using the s 2 sensitivity measure are presented in Fig. 13 for different values of e and one value of flexural rigidity. In all cases the maximum value is about 10 percent up to the final time in the simulation. Similar small values of s 2 were found for other rigidity values.

The s 3 sensitivity results are presented in Table 2 for a par ticular flexural rigidity (E/ = 7 .50), and different t 1 and e values. The maximum value of s 3 is less than 0.06. These results showing the relative insensitivity to variations in e implies that, even if the flexibility and the system con figuration does affect e, that may not significantly affect the response. Thus, it can be concluded that a coefficient of restitution value which is appropriately chosen from the rigid body data would give sufficiently accurate results, within engineering tolerances, for many situations of practical in terest. By appropriately chosen we mean that it should be selected based on the expected initial impact velocity, and the material pair in question.

Concluding Remarks

In Part 2, we have described experimental and numerical studies on the dynamics of a flexible radially rotating beam with impact.

The momentum balance method has been demonstrated to be capable of accurately predicting the dynamics of systems which consist of both rigid and elastic links undergoing im pact. This has been justified by a set of experiments performed on a specific system.

It was also shown that using a constant coefficient of restitution (provided that this value is appropriate for the first of the multiple impacts) can give almost identical results with those using a velocity dependent coefficient of restitution.

Sensitivity studies show that the model will also work reasonably well over a wide range of parameters covering most engineering situations. It can be concluded that if one is in terested in the post impact dynamic behavior of the system, the momentum balance model and the coefficient of restitu tion concept can be used with reasonable confidence.

  Location of the impact. point, Xp 0.515 m Inertia of the flexible beam, Is 0.0042 kgm2 Total Inertia of the rigid shaft (including motor and encoder shafts), IR 8.610-5kgm2 Inertia ratio, (3 = IR/Is 0.002 Modal damping for the flexible beam 0.0 Joint viscous damping coefficent, c 0.00537 kgm2 / s

Fig. 2

 2 Fig. 2 Free response of the cantilever beam, (a) simulation, (b) experiment

Fig. 3 Fig. 4 Fig. 5 Fig. 6

 3456 Fig. 3 Angular velocity of the rigid shaft, free fall (a) simulation, (b) experiment 0. 00 I 0

4. 2

 2 Results with Torque Input.

  8(a) and 8(b)]. Strains at the other three locations exhibit a very similar behavior. The momentum balance method appears to give very good results in terms of the elastic motion. Another simulation was performed in which the value of the coefficient of restitution was determined based on the velocity of the first impact and held constant throughout the simula tion. Results are presented in Figs. 9 and 10 for the angular � .... .._

Fig. 7 Fig. 8 Fig. 9 Fig. 10

 78910 Fig. 7 Angular velocity of the rigid shaft, (a) simulation, (b) experiment 0.002

Fig. 11 Fig. 12

 1112 Fig. 11 Elastic peak deflection sensitivity, s1

Fig. 13

 13 Fig. 13 Tip displacement sensitivity, s2, for El = 7.50

Table 1 setup

 1 Dimensions and structural properties in the experimental

	Mass per unit length, p	0.0855 kgfm
	Flexural rigidity, EI	5.50 Nm2

Table 2

 2 Global tip displacement sensitivity s3, forE/= 7.5Coefficient ofValues of S3 forRestitution, e tl = 0.10 tl = 0.15 tl = 0.2

				0 20
	0.4	0.0034	0.0033	0.0028
	0.5	0.0059	0.0081	0.0089
	0.6	0.0	0.0	0.0
	0.7	0.0132	0.0335	0.0527
	0.8	0.0076	0.0089	0.0092
	0.9	0.0173	0.0280	0.0378
	1.0	0.0130	0.0104	0.0134

Acknowledgments

The authors would like to thank Mr. Gerald E. Lilly of Eastman-Kodak Company for his help with the high speed video recordings, and Mr. Steve Culp of The University of Michigan, for his assistance with setting up the experiments.


Trianguler le plan projectif réel

Résumé : Nous considérons le calcul de la triangulation d'un ensemble fini de points P = {p 1 , p 2 , . . . , p n } dans le plan projectif. Nous démontrons que la triangulation existe toujours dès lors qu'au moins six points de P sont en position générale, c'est-à-dire que trois d'entre eux ne sont jamais alignés. Nous proposons également un algorithme pour trianguler P 2 si cette condition nécessaire est remplie. À notre connaissance, c'est le premier résultat algorithmique connu pour le plan projectif réel.

Mots-clés : Géométrie algorithmique, triangulation, complexe simplicial, géométrie projective, algorithme

Introduction

The real projective plane P 2 is in one-to-one correspondence with the set of lines of the vector space R 3 . Formally, P 2 is the quotient P 2 = R 3 -{0} / ∼ where the equivalence relation ∼ is defined as follows: for two points p and p ′ of P 2 , p ∼ p ′ if p = λp ′ for some λ ∈ R-{0}.

Triangulations of the real projective plane P 2 have been studied quite well in the past, though mainly from a graph-theoretic perspective. A contraction of and edge e in a map M removes e and identifies its two endpoints, if the graph obtained by this operation is simple. M is irreducible if none of its edges can be contracted. Barnette [START_REF] Barnette | Generating triangulations of the projective plane[END_REF] proved that the real projective plane admits exactly two irreducible triangulations, which are the complete graph K 6 with six vertices and K 4 + K 3 (i.e., the quadrangulation by K 4 with each face subdivided by a single vertex), which are shown in Figure 1. Note that these figures are just graphs, i.e. the horizontal and vertical lines do not imply collinearity of the points. A diagonal flip is an operation which replaces an edge e in the quadrilateral D formed by two faces sharing e with another diagonal of D (see Figure 2). If the resulting graph is not simple, then we do not apply it. Wagner [START_REF] Wagner | Bemerkungen zum Vierfarbenproblem[END_REF] proved that any two triangulations on the plane with the same number of vertices can be transformed into each other by a sequence of diagonal flips, up to isotopy. This result has been extended to the torus [START_REF] Dewdney | Wagner's theorem for the torus graphs[END_REF], the real projective plane and the Klein bottle [START_REF] Negami | Diagonal transformations of triangulations on surfaces[END_REF]. Moreover, Negami has proved that for any closed surface F 2 , there exists a positive integer N (F 2 ) such that any two triangulations G and

G ′ on F 2 with |V (G)| = |V (G ′ )| ≥ N (F 2 )
can be transformed into each other by a sequence of diagonal flips, up to homeomorphism [START_REF] Negami | Diagonal flips in triangulations on surfaces[END_REF]. Mori and Nakamoto [START_REF] Mori | Diagonal flips in hamiltonian triangulations on the projective plane[END_REF] gave a linear upper bound of (8n -26) on the number of diagonal flips needed to transform one triangulation of P 2 into another, up to isotopy. There are many papers concerning with diagonal flips in triangulations, see [START_REF] Negami | Diagonal flips in triangulations on closed surfaces, a survey[END_REF][START_REF] Eppstein | Happy endings for flip graphs[END_REF] for more references. In this paper, we address a different problem, which consists in computing a triangulation of the real projective plane, given a finite point set P = {p 1 , p 2 , . . . , p n } as input.

Definition 1.1 Let us recall background definitions here. More extensive definitions are given for instance in [START_REF] Afra | Topology for Computing[END_REF][START_REF] Henle | A Combinatorial Introduction to Topology[END_REF]. An (abstract) simplicial complex is a set K together with a collection S of subsets of K called (abstract) simplices such that:

1. For all v ∈ K, {v} ∈ S. The sets {v} are called the vertices of K.

If

τ ⊆ σ ∈ S, then τ ∈ S.
Note that the property that σ, σ ′ ∈ K ⇒ σ ∩ σ ′ ≤ σ, σ ′ can be deduced from this.

σ is a k-simplex if the number of its vertices is k + 1. If τ ⊂ σ, τ is called a face of σ.
A triangulation of a topological space X is a simplicial complex K such that the union of its simplices is homeomorphic to X.

All algorithms known to compute a triangulation of a set of points in the Euclidean plane use the orientation of the space as a fundamental prerequisite. The projective plane is not orientable, thus none of these known algorithm can extend to P 2 .

We will always represent P 2 by the sphere model where a point p is same as its diametrically opposite "copy" (as shown in Figure 3(a)). We will refer to this sphere as the projective sphere. A triangulation of the real projective plane P 2 is a simplicial complex such that each face is bounded by a 3-cycle, and each edge can be seen as a greater arc on the projective sphere. We will also sometimes refer to a triangulation of the projective plane as a projective triangulation.

a point p q r p r q z=0 (a) (b) V(p) V(q) V(r) Figure 3: (a)
The sphere model of P 2 . (b) △pqr separated from its "copy" by a distinguishing plane in R 3 .

Stolfi [START_REF] Stolfi | Oriented Projective Geometry: A Framework for Geometric Computations[END_REF] had described a computational model for geometric computations: the oriented projective plane, where a point p and its diametrically opposite "copy" on the projective sphere are treated as two different points. In this model, two diametrically opposite triangles are considered as different, so, the computed triangulations of the oriented projective plane are actually not triangulations of P 2 . Identifying in practice a triangle and its opposite in some data-structure is not straightforward. Let us also mention that the oriented projective model can be pretty costly because it involves the duplication of every point, which can be a serious bottleneck on available memory in practice.

The reader should also note that obvious approaches like triangulating the convex hull of the points in P and their diametrically opposite "copies" (on the projective sphere) separately will not work: it may happen that the resulting structure is not a simplicial complex (see Figure 4 for the most obvious example), so, it is just not a triangulation (see definition 1.1).

Terminology and Notation

We assume that the positions of points and lines are stored as homogeneous coordinates in the real projective plane. Positions of points will be represented by triples (x, y, z) (with z = 0) and

INRIA p q r p r q V(p) V(q) V(r)
Figure 4: The convex hull of duplicated points is not a triangulation. their coordinate vectors will be denoted by small letters like p, q, r, . . .. Positions of lines will also be represented by triples [x, y, z] but their coordinate vectors will be represented by capital letters like L, M, N, . . .. We shall also state beforehand whether a given coordinate vector is that of a line or a point to avoid ambiguity. Point p and line M are incident if and only if the dot product of their coordinate vectors p • M = 0. If p and q are two points then the line L = pq can be computed as the cross product p × q of their coordinate vectors. Similarly the intersections of two lines M and N can be computed as the cross product M × N of their coordinate vectors.

We denote the line in R 3 corresponding to a point p in P 2 by V (p). A plane in R 3 which separates △abc from its diametrically opposite "duplicate copy" on the projective sphere will be referred to as a distinguishing plane for the given triangle (see Figure 3(b)). Note that a distinguishing plane is not unique for a given triangle. Also note that such a plane is defined only for non-degenerate triangles on the real projective plane.

Contents of the Paper

We first prove a necessary condition for the existence of a triangulation of the set P = {p 1 , p 2 , . . . , p n } of P 2 . More precisely, we show that such a triangulation always exists if at least six points in P are in general position, i.e., no three of them are collinear. So if the number of points n in P is very large, the probability of such a set of six points to exist is high, implying that it is almost always possible to triangulate P 2 from a point set.

We design an algorithm for computing a projective triangulation of P if the above condition holds. The efficiency of the algorithm is not our main concern in this paper. The existence of an algorithm for computing a triangulation directly in P 2 is our main goal. As far as we know, this is the first computational result on the real projective plane.

The paper is organized as follows. In section 2, we devise an "in-triangle" test for checking whether a point p lies inside a given △abc. In section 3, we first prove that a triangulation of P 2 always exists if at least six points in the given point set P are in general position. We then describe our algorithm for triangulating P 2 from points in P. Finally, in section 4, we present some open problems and future directions of research in this area.

The Notion of "Interior" in the Real Projective Plane

It is well-known that the real projective plane is a non-orientable surface. However, the notion of "interior" of a closed curve exists because the projective plane with a cell (any figure topologically equivalent to a disk) cut out is topologically equivalent to a Möbius band [START_REF] Henle | A Combinatorial Introduction to Topology[END_REF]. For a given triangle on the projective plane, we observe that its interior can be defined unambiguously if we associate a distinguishing plane with it. The procedure for associating such a plane with any given triangle will be described in Section 3. For now we will assume that we have been given △pqr along with its distinguishing plane in R 3 . We further assume for simplicity that this plane is z = 0 for the given △pqr (as shown in Figure 3(b)). Consider the three lines V (p), V (q) and V (r) in R 3 . These lines give rise to four double cones, three of which are cut by the distinguishing plane. We define the interior of △pqr as the double cone in R 3 which is not cut by its distinguishing plane. Based on the above definition, we define a many-one mapping s : P 2 → R 3 from points in P 2 to points in R 3 as follows:

s(p) = s(x, y, z) =    (1, x z , y z ) z = 0; (0, 1, y x ) z = 0, x = 0; (0, 0, 1) z = 0, x = 0.
Given three points a = (x 0 , y 0 , z 0 ), b = (x 1 , y 1 , z 1 ), c = (x 2 , y 2 , z 2 ) and a point p = (x, y, z), p lies inside △abc if

sign s 0 s 1 s + sign s 1 s 2 s + sign s 2 s 0 s = ±3 (1)
and it lies on the perimeter of △abc if

sign s 0 s 1 s + sign s 1 s 2 s + sign s 2 s 0 s = ±2 (2)
Here s i = s(x i , y i , z i ) for i = 0, 1, 2, and s = s(x, y, z). The function sign(m) returns 1 if m is positive, 0 if m is zero, and -1 if m is negative. The reader should note that similar to the oriented projective model [START_REF] Stolfi | Oriented Projective Geometry: A Framework for Geometric Computations[END_REF], there is no notion of interior when all a, b, c and p are at infinity. We now consider the case when the distinguishing plane in R 3 is αx + βy + γz = 0, where α, β, γ are arbitrary constants. We use a linear transformation matrix M for transforming the given plane into the plane z = 0 according to the equation M • p ′ = p, where orientation is preserved. This transformation takes the coordinate vector p of a point to the vector p ′ . Now the s-mapping of equation ( 1) can be used for the "in-triangle" test with the new coordinate vectors, as described above. For the case when γ = 0, we have

M =      0 -(β 2 +γ 2 ) √ (β 2 +γ 2 )(α 2 +β 2 +γ 2 ) α √ (α 2 +β 2 +γ 2 ) γ √ β 2 +γ 2 αβ √ (β 2 +γ 2 )(α 2 +β 2 +γ 2 ) β √ (α 2 +β 2 +γ 2 ) -β √ β 2 +γ 2 αγ √ (β 2 +γ 2 )(α 2 +β 2 +γ 2 ) γ √ (α 2 +β 2 +γ 2 )      (3)
For the case when γ = 0, β = 0, we have

M =     β √ α 2 +β 2 0 α √ α 2 +β 2 -α √ α 2 +β 2 0 β √ α 2 +β 2 0 -1 0     (4)
Finally, we have the case when γ = β = 0. In this case, we simply make the X-axis the new Y ′ -axis, the Y -axis the new Z ′ -axis, and the Z axis the new X ′ -axis. So our tranformation matrix M is as follows:

M =   0 0 1 1 0 0 0 1 0   (5)
Note that all the transformation matrices given by equations ( 2), ( 3) and (4) are orthogonal matrices, i.e., M -1 = M T .

Computing the Projective Triangulation of a Point Set

We now proceed to discuss our algorithm for triangulating the real projective plane given a point set P = {p 1 , p 2 , ..., p n } as input. We number the points in this section for diagramatic clarity. We first prove the following simple result for point sets in P 2 : Proof: It is easy to see that the lemma holds if all points in P are collinear. So we may safely assume that this is not the case. We prove the above lemma by the method of contradiction.

Assuming that no set of (n -1) points are collinear. Consider a set of four points {1, 2, 3, 4} as shown in Figure 5(a). Since among every set of four points, at least three are collinear, so we assume that 1, 2 and 3 are collinear. Now consider a fifth point 5 instead of 3, and assume that it does not lie on the line L = 12. From the given condition, it must lie on the line M = 14. But now no three points among the set {2, 3, 4, 5} are collinear, a contradiction! Corollary 3.1 If (n -1) points are not collinear in the given point set P, then there exists a set of four points in P, no three of which are collinear.

We call such a set of four points a K 4 -quadrangulation. Corollary 3.1 states that every point set P in which no (n -1) points are collinear contains a K 4 -quadrangulation. We now make the following important observation that such a point set can be used to construct a triangulation of the projective plane (see Figure 5(b,c)). We have the following lemma: Lemma 3.2 A K 4 -quadrangulation can be used to construct a projective triangulation.

Proof: Consider the points of the K 4 -quadrangulation on the projective sphere and construct the lines (great circles) {12, 13, 14, 24, 23, 34} (see Figure 5(b,c)). The intersection of these six lines define three more points {p, q, r}. We call these points pseudo-points because these may or may not be points in P. It is now easy to see that the resulting triangulation is a simplicial complex and is isomorphic to the projective triangulation shown in Figure 1(a). The reader should observe that every triangle in the above triangulation has precisely two copies on the projective sphere which are diametrically opposite (see Figure 5(b,c)). So it now becomes possible to associate a distinguishing plane with each triangle in the above triangulation unambiguously. For every △abc in the projective triangulation, we can take the plane through the center O of the projective sphere and parallel to the plane passing through the end-points a, b, c of one copy of △abc. Given a query point u, we can now determine the triangle inside which it lies. We will use this fact quite extensively in our algorithm. The procedure described above is incomplete in the sense that we triangulate the real projective plane with the help of some pseudo-points. We now give a necessary condition for computing a projective triangulation from a point set P. The reader should observe that every triangle in Figure 5(b,c) is incident with exactly one pseudo-point. We will refer to the set of triangles incident to one pseudo-point as a region. Note that any two regions have the same set of vertices. For constructing a projective triangulation from P we will initially take help of pseudo-points, but we will go on deleting them as their use is over. We now present the following lemma: are in different regions of the projective triangulation formed by S, then it is possible to triangulate the projective plane using these six points, unless (n -2) points in P are collinear.

Proof: We give a constructive proof of the above statement. We first construct a projective triangulation with the set {1, 2, 3, 4} (as described above). Suppose points 5 and 6 lie in the regions associated with the pseudo-points p and q respectively (see Figure 6(a,b)). We now add point 5 and make it adjacent to the vertices of its bounding region, deleting the pseudo-point p and the edges it was incident with. The newly added edges are shown by dashed lines. The pseudo-points have also been kept for better understanding. We now add point 6 and delete the corresponding pseudo-point q and the edges it was incident with. Now we intend to delete the pseudo-point r and construct a valid projective triangulation using only points in P. Here we make the important observation that either the edge 12 or 34 can be flipped. To see this, note that if flipping of neither of these edges was possible, then 6 must lie to the "left" (as shown in Figure 6(a)) or "right" (as shown in Figure 6(b)) of both the lines 52 and 54, in which case flipping of edges would induce crossings. (Note that we refer to a point being on the "left" or "right" of a line only locally with respect to front half of the projective sphere.) However, the edge 24 lies in between these two lines and 6 cannot lie to its left (resp. right). Thus, our claim holds. Suppose the edge 12 can be flipped. We then construct a valid projective triangulation by flipping 12, deleting the pseudo-point r and adding the edge 12 in that region. Observe that the projective triangulation constructed is isomorphic to that shown in Figure 1(b). In the event that flipping of neither 12 nor 34 is possible, all four points 5, 2, 4, 6 must be collinear. Since such a flip is also not possible with any other point in P, they must all lie on the line 524, implying that (n -2) of the points in P are collinear.

We will refer to such a K 4 -quadrangulation which has two points of P in different regions as a canonical set. We now have almost all the basic tools required for triangulating the real projective plane from a point set P. All that we need to characterize is the existence of a canonical set. So far we have not used anywhere the assumption that at least six points in P are in general position. It turns out that there always exists a canonical set in P in this case. We have the following lemma: Lemma 3.4 If at least six points in P are in general position, then there exists a canonical set.

Proof: We prove this lemma by the method of contradiction. We assume that the lemma does not hold, so for every K 4 -quadrangulation in P, all other points of P are in the same region. Consider a K 4 -quadrangulation {1, 2, 3, 4} in P. Suppose we add two more points 5 and 6, and they lie in the same region (as shown in Figure 6(c)). Now consider the K 4 -quadrangulation formed by {4, 6, 3, 5}. If this is to satisfy the property that all points in P lie in exactly one of its regions, then it is easy to see that 2 must lie on or to the right of the line 54. But now 2 and 6 lie in different regions of the K 4 -quadrangulation {4, 5, 1, 3}, a contradiction! So we now have a procedure for triangulating the real projective plane given a point set P with at least six points in general position. We summarize our results in the following theorem: Theorem 3.1 Given a point set P = {p 1 , p 2 , . . . , p n } with at least six points in general position, it is always possible to construct a projective triangulation.

We now present our algorithm which outputs a triangulation of the real projective plane given a point set P = {p 1 , p 2 , . . . , p n } with at least six points in general position.

1. Find a set S = {1, 2, 3, 4, 5, 6} of six points such that no three points in S are collinear.

2. Construct a projective triangulation with the set S. Associate distinguishing planes with every triangle of the triangulation.

3. for all points p ∈ P\S do

4.

Identify the triangle △abc in which p lies.

5.

Make p adjacent to the vertices a, b and c. Make the distinguishing plane of △apb, △bpc, and △cpa the same as that for △abc.

6. end for 7. return(triangulation of P 2 ).

There are two possible approaches for finding the set S in step 1. In the first approach, we arbitrarily choose a starting point q and initialize our set S = {q}. For any point p ∈ P\S, we add p in S if p is not collinear with any two points in S. We stop when S contains six points. It may happen that we are not able to find such a set S of six points if we start with any random starting point q. So we iterate over all points in P for choosing the starting point. This approach has a worst-case time complexity of O(n 2 ). A slightly better approach can be adopted for performing step 1, which works in O(n) time if we assume that the minimum line cover of the point set P is greater than 4. In this approach, we first choose any two points 1 and 2. Let the line defined by them be L. We delete all other points in P on L. We now choose two more points 3 and 4. Let the line defined by them be M . We delete all other points in P on M . We also delete all other points on N (the line defined by 1 and 3) and T (the line defined by 2 and 4). Now choose two more points 5 and 6. We now have the required set S = {1, 2, 3, 4, 5, 6}. It is easy to see that this approach takes O(n) time if the minimum line cover of P is greater than 4. The above two approaches work reasonably well for most point sets. However, for certain point sets, it may happen that both these approaches fail to find such a set S. We are currently unaware of an optimal method for finding such a set which works in all cases. We believe that some approach similar to that used for solving the "ordinary line" problem can be adopted for finding the same (see for instance [START_REF] Kelly | On the number of ordinary lines determined by n points[END_REF][START_REF] Mukhopadhyay | On the ordinary line problem in computational geometry[END_REF][START_REF] Devillers | Finding an ordinary conic and an ordinary hyperplane[END_REF]). After having found such a set S, we find a canonical set within S by a procedure similar to that described in Lemma 3.4.

Once we have a canonical set, constructing the projective triangulation in step 2 takes O(1) time. We store the triangulation in a DCEL so that addition and deletion of edges and vertices takes O(1) time. The loop in steps 3-6 runs once for every point p ∈ P\S. Inside the loop, we use our "in-triangle" test (as described in Section 2) for testing whether a point lies inside a given triangle. We use a procedure similar to that described by Devillers et al. in [START_REF] Devillers | Walking in a triangulation[END_REF] for identifying △abc inside which p lies. We first choose any arbitrary vertex t of the current projective triangulation. We then identify the triangle whose interior is intersected by the line L = tp. This test is performed by checking for all edges E of all triangles sharing vertex t whether the intersection of L and the line described by E lies inside the given triangle (see Figure 7(a)). After having identified the starting triangle, we move to its neighbor sharing the edge E. In this way, we "walk" in the triangulation along the line L. We stop when p lies inside the current triangle. Although this method of "walking" in a triangulation has a worst-case time complexity of O(n), it is reasonably fast for most practical purposes. So the loop takes a total of O(n 2 ) steps. Thus, our algorithm computes a projective triangulation from a given point set P in O(n 2 ) steps.

As mentioned in the introduction, the complexity of the algorithm is not our main concern in the present paper. Still, note that our algorithm is incremental, which is an important property in practice. O(n 2 ) is a standard worst-case complexity for incremental algorithms computing triangulations in the Euclidean plane. After step 2, instead of inserting the points incrementally, we could do the following1 : for each point, find the triangular face of the initial triangulation containing it. Then, in each of these faces, triangulate the set of points using the usual affine method. This can be done since the convex hulls of subsets of points in a triangular face of the initial triangulation can be defined with the help of distinguishing planes. This yields an optimal O(n log n) worst-case time (non-incremental) algorithm.

Conclusion and Open Problems

It woud be interesting to check whether the metric on P 2 allows to define a triangulation of the projective plane that would extend the notion of Delaunay triangulation, which is well-known in the Euclidean setting. Then, extending the randomized incremental insertion with a hierarchical data-structure such as [START_REF] Devillers | The Delaunay hierarchy[END_REF] to the projective case, if possible, would lead to an incremental algorithm with better theoretical (and practical, too) complexity.

Also, problems like the Minimum Weight Triangulation [START_REF] Mulzer | Minimum weight triangulation is NP-hard[END_REF], Minmax Length Triangulation [START_REF] Edelsbrunner | A quadratic time algorithm for the minmax length triangulation[END_REF], etc., may have meaning even on the real projective plane. The Minimum Weight Triangulation problem was neither known to be NP-Hard nor solvable in polynomial time for a long time [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF]. This open problem was recently solved and was shown to be NP-Hard by Mulzer and Rote [START_REF] Mulzer | Minimum weight triangulation is NP-hard[END_REF]. The Minmax Length Triangulation problem asks about minimizing the maximum edge length in a triangulation of a point set P. This problem was shown to be solvable in time O(n 2 ) by Edelsbrunner and Tan [START_REF] Edelsbrunner | A quadratic time algorithm for the minmax length triangulation[END_REF]. It would be interesting to analyze the complexity of these problems on the real projective plane P 2 .

Figure 1 :

 1 Figure 1: The two irreducible triangulations of P 2 .

Figure 2 :

 2 Figure 2: A diagonal flip.

Figure 5 :

 5 Figure 5: (a) Every point set P has a K 4 -quadrangulation, unless (n -1) points are collinear. (b) and (c) A K 4 -quadrangulation is sufficient for triangulating the real projective plane.

Figure 6 :Lemma 3 . 3

 633 Figure 6: (a) and (b) Symmetric cases for constructing a projective triangulation. (c) A canonical set always exists when six points in P are in general position.Lemma 3.3 If there exists a set of six points (say, {1, 2, 3, 4, 5, 6}) in a given point set P such that four of them (say, S = {1, 2, 3, 4}) form a K 4 -quadrangulation and the other two (say, {5, 6}) are in different regions of the projective triangulation formed by S, then it is possible to triangulate the projective plane using these six points, unless (n -2) points in P are collinear.

Figure 7 :

 7 Figure 7: (a) "Walking" in a projective triangulation. (b) Two copies of the edge pq, only one is cut by the distinguishing plane.

RR n 6296

as suggested by an anonymous reviewer INRIA

Acknowledgements

The authors would like to thank Olivier Devillers for his valuable suggestions and helpful discussions.