Jean-Philippe Condomines
email: jean-philippe.condomines@recherche.enac.fr.

Cédric Seren
email: cedric.seren@onera.fr.

Gautier Hattenberger
email: gautier.hattenberger@recherche.enac.fr.

Nonlinear state estimation using an invariant unscented Kalman filter

In this paper, we proposed a novel approach for nonlinear state estimation, named π-IUKF (Invariant Unscented Kalman Filter), which is based on both invariant filter estimation and UKF theoretical principles. Several research works on nonlinear invariant observers have been led and provide a geometrical-based constructive method for designing filters dedicated to nonlinear state estimation problems while preserving the physical properties and systems symmetries. The general invariant observer guarantees a straightforward form of the nonlinear estimation error dynamics whose properties are remarkable. The developed π-IUKF estimator suggests a systematic approach to determine all the symmetry-preserving correction terms, associated with a nonlinear state-space representation used for prediction, without requiring any linearization of the differential equations. The exploitation of the UKF principles within the invariant framework has required the definition of a compatibility condition on the observation equations. As a first result, the estimated covariance matrices of the π-IUKF converge to constant values due to the symmetry-preserving property provided by the nonlinear invariant estimation theory. The designed π-IUKF method has been successfully applied to some relevant practical problems such as the estimation of Attitude and Heading for aerial vehicles using low-cost AH reference systems (i.e., inertial/magnetic sensors characterized by low performances).

Nomenclature e i

Orthonormal basis vector, i ∈ [ [ 1; 3 ]] ω Instantaneous angular velocity vector, rad/s ω 0

Reference angular velocity vector in R Vertical speed component (body frame) A Constant gravity vector in the North-East-Down (NED) coordinate system i.e., A = ge 3 , m/s 2 a Specific acceleration vector, m/s 2 V Velocity vector of the center of gravity (Earth frame), m/s

I. Introduction

The attitude estimation problem for nonlinear dynamic system is an important research topic and is a major concern in the aerospace engineering community. [START_REF] Crassidis | A survey of Nonlinear Attitude Estimation Methods[END_REF] The Extended Kalman Filter (EKF), and its numerous variants e.g., the multiplicative EKF, is the most widely used signal processing methodology. Its principle relies on a first-order linearization of the nonlinear model equations so that the traditional linear Kalman Filter (KF) can be applied. However, it has two well known drawbacks: [START_REF] Xiong | Performance evaluation of UKF-based nonlinear filtering[END_REF] (1) the first-order linearization can introduce large errors in the mean and covariance of the state vector and (2) the calculation of Jacobian matrices is nontrivial for many applications. To overcome this problem, the Unscented KF 3 (UKF) is an efficient derivative free estimation algorithm for determining approximate solutions to discretetime nonlinear optimal filtering problems. It has been successfully applied to various practical problems and it has been shown that it outperforms the EKF in many cases. [START_REF] Van Der Merwe | Sigma-point Kalman filters for probabilistic inference in dynamic state-space models[END_REF][START_REF] Sarkka | On Unscented Kalman Filtering for State Estimation of Continous-Time Nonlinear Systems[END_REF][START_REF] Xu | Stochastic Stability of the Continuous-Time Unscented Kalman Filter[END_REF] Recently, several research works on nonlinear invariant observers have been led and provide a geometricalbased constructive method for designing filters dedicated to the state estimation problem of dynamical systems while preserving the physical properties and systems symmetries. Such an approach can be very tedious and non-systematic for complex dynamical systems. That's why, a more general EKF-based method dedicated to gains calculation (and which exploits a re-linearization of the invariant state estimation error dynamics) has been recently developed in. [START_REF] Bonnabel | Invariant Extended Kalman Filter:Theory and Application to a Velocity-Aided Attitude Estimation Problem[END_REF][START_REF] Barczyk | Invariant Extended Kalman Filter Design for a Magnetometer-plus-GPS Aided Inertial Navigation System[END_REF][START_REF] Martin | Generalized Multiplicative Extended Kalman Filter for Aided Attitude and Heading Reference System[END_REF] Unfortunately, the Invariant EKF (IEKF) is based on a kind of second-order approximations of the invariant estimation error. Once again, this technique can be very complex to implement due to this re-linearization. The contribution of this paper is to propose a novel and original approach for nonlinear state estimation, named π-IUKF (Invariant Unscented Kalman Filter), which is based on both invariant filter estimation and UKF theoretical principles. The developed π-IUKF estimator suggests a systematic approach to determine all the symmetry-preserving correction terms, associated with a given nonlinear state-space representation used for state prediction, without requiring any linearization of the dynamical systems equations. It is named "π-IUKF" because a compatibility condition, which makes use of a specific nonlinear function denoted "π", on the observation equations has been defined and integrated to the standard UKF algorithm. The paper is organized as follows: in section II, kinematics and nonlinear evolution equations are presented and the π-Invariant Unscented Kalman Filter is detailed in section III. To conclude, section IV shows some results obtained in simulation and validate the whole approach.

II. Dynamical system modeling

In order to tackle a wide range of applications, various implementations of flight dynamics models, in terms of assumptions and numerical techniques, therefore exist. Mini-Unmanned Aerial Vehicles (mini-UAVs) are commonly represented using a 6 Degrees of Freedom (DoF) quaternion-based kinematic model (3 DoF correspond to the translational motion (u, v, w) and the 3 remaining DoF to the rotational motion (φ, θ, ψ)). This choice is due to:

1. the high dynamical behavior of mini-UAVs which causes problems when using the typical Euler angles representation;

2. the difficulty for such systems to derive a reliable and representative aerodynamic modeling from scratch (see Fig. 1) (which could be included within a full nonlinear flight dynamics model) since mini-UAVs must remain low-cost experimental means (no CFD calculations, no wind-tunnel tests, . . . ).

A. Kinematics and Nonlinear Evolution Equations

Assuming a flat non-rotating Earth, the flying rigid body motion can be mathematically described by:

     q = 1 2 q * ω V = A + q * a * q -1 (1)
In Eq.( 1), symbol * corresponds to the quaternion product. The first expression of Eq.( 1) describes the kinematic relation of the body of the considered mini-UAV and the second one corresponds to the Newton's law. [START_REF] Farrell | The global positioning system and inertial navigation[END_REF] The instantaneous attitude of the flying mini-UAV can be generalized and deduced using the aforementioned standard quaternionial form which provides at the same time:

• a global parametrization;

• avoids the mathematical singularities inherent to Euler angles;

• and is convenient for calculations and simulations.

For more details about formulas used on quaternion in this paper, see Appendix I.

B. Comments on the Measurement and Sensor Reliability

Recent progresses in the miniaturization of the inertial instruments i.e., Inertial Navigation System (INS) and Global Positioning System (GPS) hardware have led to the design of small, low-cost integrated navigation systems. Under conditions where GPS is available, this one provides a position accuracy of approximately 3 to 10 meters. In our case, where the position is unavailable or intermittent such as in urban or indoor environments, a major issue for autonomous maneuvering and obstacles avoidance for mini-UAVs is to estimate vehicle orientation with respect to an inertial reference frame. Thus navigation quality is limited by inertial sensor performance: given the size, power and cost constraints of miniature systems. To recover navigation accuracy using low-cost sensors such as the well-known non-aided Attitude and Heading Reference Systems (AHRS), it is necessary to use, if possible, additional instruments (e.g., magnetometers, which are used to increase the heading accuracy) and/or nonlinear estimation algorithm to improve the global mini-UAVs handling qualities. The nonlinear state estimation makes use of 3 triaxial sensors which deliver 9 scalar measurement signals:

• 3 of them are associated with 3 magnetometers which provide a local measurement of Earth's magnetic field, which is known constant a and expressed in the body-fixed frame such that the vector y B = q -1 * B * q where B = [B x , B y , B z ] can be considered as an output of the observation equations;

• 3 gyroscopes allow to obtain a measurement of the instantaneous angular velocity vector denoted by

ω m ∈ R 3 s.t ω m = [ω mx , ω my , ω mz ] T ;
• and finally, 3 accelerometers give a measurement of the specific acceleration denoted by a m ∈ R 3 s.t a m = [a mx , a my , a mz ] T .

All these measurements are obviously corrupted by noises whose effects will be taken into account in §III-B.

In this non-aided AHRS, no velocity and position measurements are available (no GPS, no Pitot sensors). Thus, to keep the whole nonlinear state representation observable given these available informations, the assumption that the linear acceleration V remains small is made, i.e., V = 0. As a result, the specific acceleration vector, expressed in the body-fixed frame, can be approximated by a = -y A = -q -1 * A * q where A is the local Earth's gravity vector. This one can be compared, as an another output vector, with its corresponding measurement a m . Thereby, we consider the following system M s given by:

M s            q = 1 2 q * ω (evolution) y =   y A y B   =   q -1 * A * q q -1 * B * q   (observation)
a The magnetic field can be determined from a world magnetic model such as B is equal to [0.5156 0.0570 0.8549] at the local flight coordinates 43.617-43 • -37 (N) and 1.450-1 • -27 (E)

All the sensors embedded are low-cost and so have imperfections. The major error sources in the navigation system are due to: -the inertial instrument imperfection; -the possible incorrect navigation system initialization (e.g., on magnetometers); -and the imperfections in the local Earth's gravity vector used in computations. For both instruments, the largest errors are usually a bias instability (expressed respectively in (deg/hr) for gyros and µg for the accelerometers biases) and a wrong scale-factor value. Obviously, smaller will be the inertial sensor errors and better will be the quality of the instruments (improved accuracy of the navigation results), but higher will be the cost of the system. Fig. 2 from [START_REF] Barbour | Inertial Sensor Technology Trends[END_REF] shows INS relative strapdown system cost "projections" as a function of the inertial instruments technology and performance. Note that colored areas are function of the quantity of sensor systems produced which influences their production cost. Thus, only MEMS (Micro-Electro-Mechanical Systems) gyros and accelerometers (performing at around 100 • /hr and 100mg bias stabilities respectively) are currently suitable for a potential use as a low-cost AHRS.

Figure 2. Strapdown INS cost as a function of instrument technology

Based on these technical elements and the observation relations which composed M s , it is possible to choose how these imperfections can be modeled since some degrees of freedom exist. A first-order observability analysis -which can be led analytically by linearizing M s -(see [START_REF] Salan | Filtering Algoritms and Avionics Systems for Unmanned Aerial Vehicles[END_REF] for more details) shows that 6 additional unknown constants can be estimated without introducing inobservability. Thus, an additive constant bias vector ω b is basically considered on the angular velocity vector measurement ω m . Then, for mini-UAV safety (magnetic measurements used to estimate the mini-UAV's heading only and accelerometric imperfection estimated without relying on the disturbed magnetic measurement), only two more constant positive scaling factors a s and b s are introduced in the estimation scheme. The latter quantity b s will be used to adjust and preserve the unit norm on vector y B . All these sensor imperfections are modeled as gaussian random walks which can be physically interpreted as slowly varying parameters. The M s modeling now becomes:

M + s                              q = 1 2 q * (ω m -ω b ) ωb = 0 (evolution) ȧs = 0 ḃs = 0   y A y B   =   a s q -1 * A * q b s q -1 * B * q   (observation) (2)
where ω m is seen as an imperfect and noisy but known measured input and a m as an available imperfect and noisy measured output like B. The nonlinear state space representation corresponding to M + s can be described in a compact form such as: ẋ = f (x, u) and y = h(x, u) where:

x = [q T , ω T b , a s , b s ] T , u = [ω m ] and y = [y T
A , y T C ] are the state, input and output vectors respectively.

III. π-Invariant Unscented Kalman Filter

A. Invariant observer and system symmetries This section presents briefly the main theoretical principles of some research works dealing with dynamical system symmetries and invariant observer. [START_REF] Bonnabel | Non-Linear Symmetry-Preserving Observers on Lie Groups[END_REF] Let's consider the following nonlinear state-space representation defined by:

   dx dt = f (x, u) y = h(x, u) (3)
where the state (resp. input)(resp. output) vector belongs to an open set

X ∈ R n (resp. U ∈ R m ) (resp. Y ∈ R p , p ≤ n).
Then, we consider the transformation group defined on X × U × Y such that :

G × (X × U × Y) → (X × U × Y) (g, x, u, y) → (ϕ g (x), ψ g (u), ρ g (y)) = (X, U, Y )
where G is a Lie-group b (i.e., a finite-dimensional differentiable manifold in which multiplication and inversion group operations are also differentiable) and (ϕ g , ψ g , ρ g ) are three diffeomorphisms c parametrized by g ∈ G.

Practically, the objective is to define 3 specific coordinates transformations whose actions on state, input and output variables respectively are distinct such that the evolution (dx/dt=f (x, u)) and observation (y=h(x, u)) equations are left unchanged i.e., Ẋ = f (X, U ) and Y = h(X, U ). Thus, a dynamical system ẋ = f (x, u) will be said

G-invariant if ∃(ϕ g , ψ g ) g∈G , ∀(g, x, u) ∈ G × (X × U), f (ϕ g (x), ψ g (u)) = Dϕ g (x) • f (x, u) and G-equivariant if moreover ∃(ρ g ) g∈G acting on Y, ∀(g, x, u) ∈ G × (X × U), h(ϕ g (x), ψ g (u)) = ρ g (h(x, u)
). These last definitions mean that all state and output equations remain explicitly identical. Based on such G-invariant/equivariant dynamical systems it is possible to design invariant observers (whose properties are often remarkable) by copying the predetermined group actions on the estimated state and output spaces. By doing this, it is possible to replace the usual linear output error on which a lot of observers are based (for instance KF), and which does not preserve systems symmetries, with a nonlinear invariant one. In our inertial navigation application, this invariance feature corresponds to some physical aspects of the mini-UAV's movement such as the independent choice of the reference frame (Galilean invariances through transformations composed with rotations and translations). So, we are looking for frame transformations that leave the system identical. Considering the explicit M + s modeling given in Eq.( 2) and the Lie-group G = H 1 × R 5 (where H 1 is the unit quaternions manifold) acting on the whole state space of M + s , the following variable transformations give to our considered dynamical system the G-invariant and G-equivariant properties s.t. ∀g = (q T 0 ω T 0 a 0 c 0 ) ∈ G :

1) ϕ g x = q T ω T b a s b s T =      q * q 0 q -1 0 * ω b * q 0 + ω 0 a s • a 0 b s • b 0      2) ψ g (ω = (ω m )) = q -1 0 * ω m * q 0 + ω 0 3) ρ g y = y T A y T B T = a 0 • q -1 0 * y A * q 0 b 0 • q -1 0 * y B * q 0
These latter are equivalent to time-constant rotations and translations in both earth-and body-fixed frames. Moreover, let's consider Q = (q * q 0 ) , Ω b = (q -1 0 * ω b * q 0 + ω 0 ) and Ω m = (q -1 0 * ω m * q 0 + ω 0 ). The system M + s is indeed invariant by transformation group since:

Q = • (q * q 0 )= ( 1 2 q * (ω m -ω b )) * q 0 = 1 2 (q * (ω m -ω b )) * q 0 = 1 2 (q * ω m -q * ω b ) * q 0 = 1 2 (q * ω m * q 0 -q * ω b * q 0 ) = 1 2 (q * q 0 * q -1 0 * ω m * q 0 -q * q 0 * q -1 0 * ω b * q 0 ) = 1 2 (q * q 0 ) * ((q -1 0 * ω m * q 0 + ω 0 ) -(q -1 0 * ω b * q 0 + ω 0 )) = 1 2 Q * ((q -1 0 * ω m * q 0 + ω 0 ) -(q -1 0 * ω b * q 0 + ω 0 )) = 1 2 Q * (Ω m -Ω b ) Ωb = • (q -1 0 * ω b * q 0 + ω)= q -1 0 * ωb * q 0 = 0 • (a s • a 0 ) = ȧs • a 0 = 0 • (b s • b 0 ) = ḃs • b 0 = 0
whereas the output (observation) is G-compatible since:

ρ g (y) = a 0 • q -1 0 * y A * q 0 b 0 • q -1 0 * y b * q 0 = a 0 • q -1 0 * (a s • q -1 * A * q) * q 0 b 0 • q -1 0 * (b s • q -1 * B * q) * q 0 =        a s • a 0 • q -1 0 * q -1 * A * Q (q * q 0 ) b s • b 0 • (q -1 0 * q -1 ) =(q * q0) -1 i.e., Q -1 * B * Q (q * q 0 )        = a s • a 0 • Q -1 * A * Q b s • b 0 • Q -1 * B * Q = h ϕ g (x), ψ g (u)
Then, if we consider the particular case where : g = x -1 = [(q -1 ) T (-q * ω b * q -1 ) T a -1 s b -1 s ] T , we obtain:

ϕ x -1 (x) =      q * q -1 q * ω b * q -1 -q * ω b * q -1 a s • a -1 s b s • b -1 s      =      1 0 1 1      = e
which is equal to the neutral element e of the local transformation ϕ g . Indeed, each component of the vector (1 0 1 1) T corresponds to the neutral element of:

• the quaternion Hamilton product (1 * );

• the vector translation operation (0 +);

• and the real scalar product (1 •).

This vector forms an invariant and defines x -1 ∈ G as the invariant Lie-group parametrization to choose. This result is known as the normalization equations solution and allows to build an invariant output error denoted by E. In the followings, we will use for commodity reasons the notation:

γ(x) =      q -1 -q * ω b * q -1 a -1 s b -1 s      (4)
The general form of a symmetry-preserving invariant observer can be defined such as :

ẋ = prediction G-invariant f (x, u) + invariant correction term n i=1 K i ψ x-1 , E := h e = ϕ x-1 (x -1 ), ψ x-1 (u) -ρ x-1 (y) w i (x) (5)
In Eq.( 5), x corresponds to the estimated state vector and is an element of R n . We recognize in Eq.( 5) the typical mathematical expression of a nonlinear state estimator in which a model-based prediction whose dynamics is G-invariant and its output G-equivariant is summed with an invariant correction term. The idea is to build an invariant additive correction term base on the nonlinear gains K i (which depend on the invariances (e, ψ x-1 , ρ x-1 ) in a nonlinear way) which keeps the whole estimated state vector dynamics ẋ G-invariant such as for the dynamical system representation f . In Eq.( 5) the standard analytical expression of the invariant output error E is given. Its explict mathematical expression can be calculated directly such that:

E = h e = ϕ x-1 (x), ψ x-1 (u) -ρ x-1 (y) = A -â-1 s q * y A * q-1 B -b-1 s q * y B * q-1 = E A E B
Moreover, the w i (x) vectors form an invariant frame i.e., for each x ∈ X , (w 1 (x), • • • , w n (x)) project the invariant correction term on each component of f (x, u) (i.e., the tangent state space). In our case, considering the canonical basis of R 3 , (e i ) i∈[ [ 1;3 ] ] , the expression of these vectors which define an invariant reference frame (also called invariant natural basis in the tangent state space) is obtained by the 8 vector equations with respect to variable w(q, ωb , âs , bs ) s.t.:

Jacobian matrix of the geometrical

transformation ϕg (x)      Dϕ γ(x)      q ω b a s b s           .w(x) =                  0 e i 0 0 0 0 0                  i∈[ [ 1;3 ] ] ,                  0 0 0 0 e i 0 0                  i∈[ [ 1;3 ] ] ,                  0 0 0 0 0 0 0 1 0                  and                  0 0 0 0 0 0 0 0 1                  (6)
By inverting the Jacobian matrix Dϕ γ(x)

     q ω b a s b s     
, it follows that the invariant vector field w i (x) is given by:

     w q i (x) w ω b i (x) w as (x) w bs (x)      =                  0 e i * q 0 0 0 0 0                  i∈[ [ 1;3 ] ] ,                  0 0 0 0 q -1 * e i * q 0
Then, we have:

∀i ∈ [[ 1; n ]], K i (ψ x-1 (u), E = 0) = 0 (i.

e., no correction if no output error !)

It follows that the 1 st order Taylor series of each nonlinear gain is reduced to:

∀i ∈ [[ 1; n ]], K i (ψ x-1 (u), E) = Ki (ψ x-1 (u), E) • E
So that a linear approximation of Eq.( 5) can be obtained easily s.t.:

ẋ = f (x, u) + n i=1 Ki (ψ x-1 (u), E) • E • w i (x) (8)
Combining all these results, Eq.( 8) can be detailed as follows:

OM s                          q = 1 2 q * (ω m -ωb ) + 3 i=1 ( Kq Ai • E A + Kq Bi • E B )e i * q ωb = 3 i=1 q-1 * ( Kω b Ai • E A + Kω b Bi • E B ) * q ȧs = âs ( Kas Ai • E A + Kas Bi • E B ) ḃs = bs ( Kbs Ai • E A + Kbs Bi • E B )
where the quantities Kq Ai , Kq Bi , Kω b Ai , Kω b Bi , Kas Ai , Kas Bi , Kbs Ai and Kbs Bi are arbitrary 1×3 gain matrices whose choice will determine the convergence properties of the estimation scheme. A direct and analytical design can be done thanks to the linearization of the invariant state error denoted η(x, x) and defined such as:

η(x, x) = ϕ γ(x) (x) -ϕ γ(x) (x) =      q * q -1 q * ( ŵb -w b ) * q-1 a s /â s b s / bs      ∼ x -1 x.
Such an approach can be very tedious and non-systematic for more complex dynamical systems than the ones represented by a fully kinematic model. That's why recently, a more general EKF based method dedicated to gains calculation (exploiting a re-linearization of the invariant state estimation error dynamics) has been developed in. [START_REF] Bonnabel | Invariant Extended Kalman Filter:Theory and Application to a Velocity-Aided Attitude Estimation Problem[END_REF][START_REF] Barczyk | Invariant Extended Kalman Filter Design for a Magnetometer-plus-GPS Aided Inertial Navigation System[END_REF][START_REF] Martin | Generalized Multiplicative Extended Kalman Filter for Aided Attitude and Heading Reference System[END_REF] In that way, the Invariant EKF can be viewed as a "second-order" approximation of the invariant estimation error. Once again, this technique can be very complex to implement due to this re-linearization. As a consequence a novel and original UKF-based approach, named π-IUKF, has been defined. It addresses the approximation issue of the invariant EKF without requiring any linearization of the dynamical systems equations. It is called "pi-"IUKF because a compatibility condition, pronounced "pi", on the measurement equations has been defined and integrated to the standard UKF algorithm. This new way for calculating the matricial gains has motived our research works.

B. π-IUKF algorithm

In this section, we propose a modification of the UKF 3 algorithm for nonlinear systems possessing symmetries.

As mentioned in the previous section, this algorithm presents the advantage to compute an approximated solution of a nonlinear optimal filtering problem in discrete-time without requiring any linearization of the nonlinear systems equations. A method, named Sigma-Points Approach 14 (SPA), is used for calculating the gaussian approximation of an unknown random variable x k|k at time kT on which is applied a nonlinear transformation such as follows:

   x k|k ∼ N (x k|k , P k|k ) s.t x k+1 = f (x k , u k , k) + w k y k|k = h(x k|k , k) + v k (9)
where x ∈ R n , y ∈ R p and w, v are independent white gaussian noises. The aims of the basic SPA is to form a fixed number of deterministically chosen sigma-points χ k|k , which capture the "true" mean xk|k and covariance P k|k of the original distribution x k|k . This set of points must represent accurately the first and second order moments. These sigma-points are then propagated through the nonlinear functions Eq.( 9) providing a cloud of evolving points. The mean xk|k+1 and estimated covariance matrix P k|k+1 of the transformed points are then computed based on their statistics. Therefore, the UKF is only a straightforward extension of the SPA for computing approximate solutions to the estimation problem posed whose mathematical form is presented by:

x k+1 = f (x k , u k ) + w k (10)
z k = h(x k , u k ) + v k (11)
Concerning both evolution and observation noises (w and v), we assume that they enter into the system in an "invariant way" such as in. [START_REF] Bonnabel | Invariant Extended Kalman Filter:Theory and Application to a Velocity-Aided Attitude Estimation Problem[END_REF] Noises can be considered as additive exogeneous invariant inputs and outputs, which do not change systems symmetries, if their insertion in the modeling equations relies on Ginvariant/equivariant mathematical terms. That's why some specific operations (such as the right quaternion multiplications • * q) appear in the discrete differential equations at the beginning of §IV. Note that in reference, [START_REF] Bonnabel | Invariant Extended Kalman Filter:Theory and Application to a Velocity-Aided Attitude Estimation Problem[END_REF] the authors make an important assumption on the local diffeomorphism ρ g which is motivated by the fact that this latter application is not generally linear. This can raise problems in considering the measurement noises within the G-equivariant observation model. The exploitation of the UKF principles within the nonlinear invariant estimation framework required the definition of a particular condition on the observation equations, named compatibility condition, which corresponds to the application of a nonlinear function, denoted π, to the predicted measurement vector. The resulting operation can be interpreted as the construction of a set of new predicted measurements whose expressions can be projected in the natural base w q i (x k+1|k ) (see §III-A). Indeed, we have:

π(h(x k+1|k , u k )) = 3 i=1 ρ γ (x k+1|k ) (h(x k+1|k , u k ))w q i (x k+1|k )e i (12)
where:

γ(x k+1|k ) =       q-1 k+1|k -q k+1|k * ωb k+1|k * q-1 k+1|k â-1 s k+1|k b-1 s k+1|k      
Therefore, Eq.( 12) can be detailed for M + s as it follows :

π(h(x k+1|k , u k )) =       âs k+1|k • 3 i=1 (A * qk+1|k )e i bs k+1 • 3 i=1 (B * qk+1|k )e i       (13)
This compatibility condition ensures that the measurement prediction step will be led in an "invariant way" so that the Square-Root UKF [START_REF] Van Der Merwe | Sigma-point Kalman Filters for Nonlinear Estimation and Sensor-Fusion: Applications to Integrated Navigation[END_REF][START_REF] Van Der Merwe | The Square-Root Unscented Kalman Filter for State and Parameter-Estimation[END_REF] (SR-UKF) gains now depend on the invariances. As a result, the π-IUKF relies on the same ideas as the usual UKF but instead of using a linear correction term based on a linear output error, it uses a geometrically adapted correction term based on an invariant output error. Likewise, the covariance matrix is not updated from an initial linear state error (x 0 -x0 ), but from an intrinsic initial invariant state error η(x 0 , x 0 ). The stochastic differential equations corresponding to the π-IUKF for nonlinear systems possessing symmetries are given by the modifications of the initialization, observation and correction steps of the standard SR-UKF. Note that we use the square-root form of the UKF method because it avoids to perform a direct matrix square-root operation at each time kT . The pseudo-code of the π-IUKF algorithm is detailed in Appendix II. For commodity reasons, the addition of a given vector u to each column of a matrix M is noted M ± u. This technique provides an increasing computational efficiency as well as a numerical stability. The SR-UKF calculations are based on three linear algebra computation techniques:

-the QR decomposition (denoted: qr{•}); -the Cholesky factor updating technique (cholupdate 16 within MATLAB TM ); -and the least squares 17 (MATLAB's '\' operator).

IV. Simulation results

We now illustrate the performances reached by the developed π-IUKF method on the basis of simulated data provided by the ONERA d lab. As it was aforementioned, the driving perturbation noises are introduced within the G-invariant system of Eq.( 10) such as follows:

             q k+1 = 1 2 q k * (ω m k -ω b k ) + W q k w q k * q k w b k+1 = q k -1 * W w k w w k * q k a s k+1 = a s k W a k w a k b s k+1 = b s k W b k w b k
and the output noises enter the observation Eq.( 11) such as:

z A k+1 z B k+1 = a s k+1 q -1 k+1 * (A + V A k+1 v A k+1 ) * q k+1 b s k+1 q -1 k+1 * (B + V B k+1 v B k+1 ) * q k+1
In the last expressions, W q k , W w k , W a k and W b k (resp. V A k and V B k ) refer to the process (resp. measurement) noise covariance matrices which are similar to weighting terms. Followingly, we define two block diagonal matrices, denoted W k and V k , by concatenating these covariance matrices. The experimental data are sampled with a frequency equal to 50Hz which characterized the inertial measurement unit and the magnetometers. The measurements are corrupted by gaussian white noises whose standard deviations are set to: σ gyro = 0.2 • /s, σ accelero = 0.2 g and σ magneto = 500 nT. In order to validate our π-IUKF filter, we have also introduced a biases vector on ω m s.t. ω b = [0.1 rad/s 0.05 rad/s 0.02 rad/s] T . The two positive scalar factor values are set to a s = 1.2 and b s = 0.9 respectively.

A. Acceleration disturbance ( V = 0)

First of all, we present the acceleration V in Fig. 3(a) for which the assumption made previously ( V = 0) may be wrong. Indeed, we can observe that this acceleration goes down and up strongly at time 20s < t < 35s. Thus, in Fig. 3(b), biases, as well as scaling factors, are a little bit influenced by this inconsistency since these ones are estimated while system is evolving. However, it appears that both biases and scaling factors converge quickly to their respective true constant values.

B. Attitude estimation

In Fig. 4(a), we firstly illustrate a typical realization of the gaussian white noise stochastic processes with zero mean which perturbate the nonlinear state estimation problem to solve. Then, the real (red solid line) and estimated (blue dashed line) attitude angles (φ, θ, ψ) are plotted on the same subfigures respectively in order to be compared. These first estimation results have been obtained without introducing an initial condition error on the state vector. We can observe on the last graph that the estimation errors are small in spite of strong acceleration phases for which the assumption V = 0 is not true. Moreover, we would like to see the behavior (the local domain of attraction, even if we have no proof of convergence) with an initial condition error on the estimated state vector. In Fig. 4(b), initial condition errors are now introduced on the whole state vector. For instance, the quaternion is initialized s.t. q = [0.8 0.05 0.05 0.01] T instead of q = [1 0 0 0] T which is not an unit norm quaternion. The real (red solid line) and estimated (blue dashed line) Euler angles are represented and a zoom on the first seconds is depicted. We can observe the good convergence behavior of our developed π-IUKF algorithm, in terms of accuracy and settling time, in presence of process and measurement noises, sensors imperfections and erroneous initial conditions. C. Estimated covariance matrices P (t) of the π-IUKF vs. standard UKF Fig. 5 shows the estimation error covariance matrix P (t) calculated by the standard SR-UKF method (on the left) and by the proposed π-IUKF technique. As explained in the theoretical development part, the new π-IUKF invariant observer must have the same invariance properties as, for instance, the IEKF developed in. [START_REF] Bonnabel | Invariant Extended Kalman Filter:Theory and Application to a Velocity-Aided Attitude Estimation Problem[END_REF] If the π-IUKF estimation error covariance matrix is constant then the whole approach will be validated. The results depicted in Fig. 5 show that P (t) (which is initialized with a given value P (0)) converge, as expected, to constant values over the whole flight trajectory due to the symmetry-preserving property provided by the nonlinear estimation theory based on invariant filter. Meanwhile, the SR-UKF estimation error covariance matrices do not convergence to constant values. Obviously, they also depend on the followed trajectory which is illustrated by the peaks. This is indicative of a smaller convergence domain for this estimation technique in comparison with the new π-IUKF method. A last cross-comparison of the results show that in some steady-state interval, the values of the standard SR-UKF and π-IUKF covariance matrices are equal. This invariance property is remarkable and can be potentially exploited to characterize some flight parameters in terms of estimation acuracy which is constant for all trajectories. Indeed, these constant covariances can be viewed as a very useful characterization of the attitude angles in a complete and future robust control law design.

V. Conclusions and Perspectives

The new nonlinear state estimation approach named π-IUKF has been presented in this paper. The first results obtained in simulation appear as promising. The hybridation of the standard SR-UKF with the invariant observer theory represents a significant research topic to improve the handling qualities of mini-UAVs. Indeed, the computational cost of such an approach seems suitable given the low-cost sensor technology used so far and whose micro-processors performances increase year after year. The perspectives of these research works are numerous and various. In the future, a complete validation of this method must be led by treating the case of real data which will be gathered from real flight experimentations. Moreover, as the final objective concerns the insertion of mini-UAVs in complex environments (such as urban areas or indoor flights), future works will deal with the robustness of this developed methodology when one or several measurement signals become unavailable. In that case, some adaptations of the π-IUKF algorithm must be made in order to recover satisfying performances. These latter works aim at providing to the estimation scheme adaptive features.

  Nonlinear state estimation using an invariant unscented Kalman filter Jean-Philippe Condomines * Ecole Nationale de l'Aviation Civile -BP 54005, Toulouse Cedex 4, 31055, France Cédric Seren † ONERA -The French Aerospace Lab -BP 74025, Toulouse Cedex 4, 31055, France Gautier Hattenberger ‡ Ecole Nationale de l'Aviation Civile -BP 54005, Toulouse Cedex 4, 31055, France

Figure 1 .

 1 Figure 1. Quadrirotor (Blender, 350mm diameter and 300g weight) and aircraft (Solar Storm, 500mm span and 400g weight) min-UAVs from the ENAC (Ecole Nationale de l'Aviation Civile) University

Figure 3 .

 3 Figure 3. Biases and scaling factors locally influenced by V = 0

  Results obtained with state initial error

Figure 4 .

 4 Figure 4. Attitude estimation results -(φ, θ, ψ)

Figure 5 .

 5 Figure 5. Estimated covariance matrices P (t): SR-UKF (left) vs. π-IUKF (right)

b e.g., SE 3 (R) (rotation matrices and translation vectors of dimension 3) is a Lie-group. c Given a Lie-group, it is called diffeomorphism if it is bijective and its inverse is differentiable.

Acknowledgments

The autors of this paper thank Thierry Le Moing from the ONERA, Toulouse, France and Prof. Erwan Salaun from the Georgia Institute of Technology, USA, for their valuable helps about nonlinear estimation theory. This research is financed by the european project "D3CoS" (see: http://www.d3cos.eu/).

Appendix I: quaternions and rotations

An unit quaternion provide a convenient mathematical notation and a global parametrization for representing orientation and rotation of a rigid body in three dimensions. Indeed, for any unit quaternion q = q 0 + q = cos θ 2 + usin θ 2 and for any vector p ∈ R 3 the action of the operator

is associated to a rotation matrix R q ∈ S0(3) which is a rotation of the coordinate frame about axis u = q ||q|| through an angle θ.

But how can a quaternion, which lives in R 4 , operate on a vector, which lives in R 3 ?

Note that a vector p ∈ R 3 can be viewed as a pure quaternion whose real part is zero 0 p . Thus, when the real part is a scalar denoted p 0 ∈ R the quaternion p is given as :

We can utilize the inner product and cross product of two vector in R 3 to write the above (non commutative) quaternion product * in more concise form:

, and (p * q) -1 = q -1 * p -1

We have to useful formulas

Appendix II: π-IUKF algorithm

The pseudo-code for the π-UKF with additive noise:

• parameters definition:

the scalar λ : is a parameter set to capture the moments of the given distribution. √ P k|k : is shorthand notation for a matrix S k|k such that P k|k = S k|k S T k|k .

W (m) , W (c) : are weights associated to the mean and covariance distributions s.t :

(m,c) = 1. W k , V k : are the predifined covariance matrices of the process and observation noises.

• Initialization:

• For kT = 1, . . . , ∞ :

1. Calculate sigma-points:

2. Time-update equations:

S k+1|k = qr W

3. Measurement-update equations:

k+1|k -ẑk+1|k


B Picinbono
email: bernard.picinbono@lss.supelec.fr

Simulations of some Doubly Stochastic Poisson Point Processes

Keywords: Point processes, lifetime, counting

Computer simulations of point processes are important either to verify the results of certain theoretical calculations that can be very awkward at times, or to obtain practical results when these calculations become almost impossible. One of the most common methods for the simulation of nonstationary Poisson processes is random thinning. Its extension when the intensity becomes random (doubly stochastic Poisson processes) depends on the structure of this intensity. If the random density takes only discrete values, which is a common situation in many physical problems where quantum mechanics introduces discrete states, it is shown that the thinning method can be applied without error. We study in particular the case of binary density and we present the kind of theoretical calculations that then become possible. The results of various experiments realized with data obtained by simulation show a fairly good agreement with the theoretical calculations.

Introduction

Point processes (PP) play an important role in various domains of Science, Technology and Economics. In the one-dimensional case, which is the only one considered in the following pages, they are convenient to describe events that arrive randomly in time either at a microscopic scale (emission of photons or electrons) or at a macroscopic scale (traffic problems). There exists abundant literature for their presentation as, for example [START_REF] Cox | Point Processes[END_REF][START_REF] Snyder | Random Point Processes in Time and Space[END_REF][START_REF] Lowen | Fractal-Based Point Processes[END_REF]. Explicit calculations of mathematical properties of PPs are often very complicated to the effect that, in a number of cases, the only possible approach for their study is by using simulation procedures, which is one of the main purposes of this paper.

They exist various papers or books devoted to simulation of Poisson processes.

In the stationary case the problem is quite straightforward because the distances between successive points are IID positive random variables (RV) with an exponential distribution. The simulation of a stationary Poisson process is then reduced to the generation of a sequence of IID positive exponential RVs, which is a very simple problem.

The question becomes more difficult in the case on nonstationary (or nonhomogeneous) Poisson processes. A nonstationary Poisson process is entirely defined by its density λ(t), (λ(t) ≥ 0), sometimes also called intensity or rate function. It is often useful to introduce the cumulative rate function Λ(t) = t 0 λ(θ)dθ, which is the mean value of the number of points of the process in the interval [0, t]. If the inverse function Λ -1 (s) of Λ(t) can be calculated explicitly, this provides a simple simulation method because the nonstationary Poisson process can be obtained by an inversion procedure from a stationary Poisson process. The use of this method, however, is often impossible because of the lack of knowledge of Λ -1 (s).

In this case, the most convenient method for the simulation of a nonstationary Poisson process is by using the thinning procedure. As said in [START_REF] Devroye | Non-Uniform Random Variate Generation[END_REF], p. 251, "the major work on simulation of nonhomogeneous Poisson processes is [START_REF] Lewis | Simulation of nonhomogeneous Poisson processes by thinning[END_REF]". Until recently various authors have presented some improvements concerning simulations of nonstationary Poisson processes, especially for applications in problems related to insurance and risks [START_REF] Harrod | Numerical methods for realizing nonstationary Poisson processes with piecewise-constant instantaneous-rate functions[END_REF][START_REF] Burnecki | Simulation of risk Processes[END_REF] and also in biology [START_REF] Shinomoto | Modeling spiking behavior of neurons with time-dependent Poisson processes[END_REF].

In this paper we are interested in simulation problems of doubly stochastic Poisson processes (DSPP). This expression refers to nonstationary Poisson processes in which the density λ(t) becomes a positive random signal λ(t, ω),

where ω stands for randomness. The first and most important use of such PPs was for Statistical Optics [START_REF] Saleh | Photoelectron Statistics[END_REF][START_REF] Lowen | Fractal-Based Point Processes[END_REF][START_REF] Picinbono | Some properties of point processes in statistical optics[END_REF], which was the origin of various studies in optical communications and information theory, some of them being summarized in [START_REF] Snyder | Random Point Processes in Time and Space[END_REF].

DSPPs, however, are not limited to physics, and also appear in numerous other domains such as, for example, management, economics, finance and traffic problems [START_REF] Burnecki | Simulation of risk Processes[END_REF][START_REF] Centanni | A Monte Carlo approach to filtering for a class of marked doubly stochastic Poisson processes[END_REF][START_REF] Lowen | Estimation and simulation of fractal stochastic point processes[END_REF].

Since DSPPs are extensions of nonstationary Poisson processes, it seems natural to extend for their simulation the thinning method introduced in the original and basic paper [START_REF] Lewis | Simulation of nonhomogeneous Poisson processes by thinning[END_REF]. This was explicitly announced in the last sentence of this paper saying "It is also possible to extend the method of thinning to simulation of doubly stochastic or conditioned Poisson processes. This will be discussed elsewhere." At least to our knowledge, this intention was never realized but realized some years late in [START_REF] Ogata | On Lewis simulation method for point processes[END_REF]. Another attempt in this direction was presented in [START_REF] Apanasovich | Computer simulation of photoevents[END_REF], but the results are not convincing and the theoretical presentation is rather obscure.

We shall below indicate other more recent approaches of this problem.

The paper is organized as follows. In Section 2 we recall the principles of the thinning method introduced in [START_REF] Lewis | Simulation of nonhomogeneous Poisson processes by thinning[END_REF] for the simulation of nonstationary Poisson processes. We indicate how this procedure can be extended for DSPPs and we present the problems appearing for this extension due to the fact that the density of a DSPP is random. In Section 3 we present a statistical model of random density for which these problems can be solved without approximation. This model is constructed with the help of a random process called random telegraph signal (RTS) widely used in many other areas of physics and signal processing. In Section 4 we calculate some statistical properties of the integral of the RTS that are useful for the following mathematical analysis. Various computer simulations confirm the validity of these calculations. In Section 5 we show that the model of random density with discrete values allows us to use the thinning method for the computer simulations of a DSPP and various experimental measurements on these simulations are presented in the last section. We show especially that measurements of moments of the random number of points in some time intervals or of PDFs of the lifetimes of such PPs yield results very similar to those obtained from the calculations.

Principles and problems of the thinning method

The principle of the thinning method is derived from the following result.

Consider a Poisson stationary point process P P S of density µ and let λ(t) be a bounded nonnegative function satisfying λ(t) < µ. Let T i be the random points (or time instants) of P P S noted in increasing order (T i < T i+1 ). To each T i we associate a thinning coefficient R i which is a Bernoulli RV taking only the values 0 or 1. We assume that these RVs are independent of each other and also independent of P P S . Furthermore we suppose that, for a given realization of P P S in which the T i s take the values T i = t i , the probability that R i = 1 is equal to λ(t i )/µ. Finally we suppose that if R i = 0 the point t i is erased and that it is not erased if R i = 1. The non erased points of P P S constitute a new point process P and it results from all these assumptions that P is a nonstationary Poisson process of density λ(t) [START_REF] Lewis | Simulation of nonhomogeneous Poisson processes by thinning[END_REF].

From this result there are various ways to obtain P from computer simulations. The simplest one is to adapt to this problem the techniques already used for the analysis of dead time effects in PPs [START_REF] Picinbono | Output dead-time in point processes[END_REF] which also leads to a problem where an initial PP P I is transformed into another one P by thinning. In this last case, however, the statistical properties of the coefficients R i are defined from the physical structure of the dead time effects, while in our problem they are defined from the previous assumptions.

Suppose now that P is a DSPP defined by a random density λ(t, ω). This means that for each value of ω, P is a nonstationary Poisson process of density λ(t, ω). We can then apply the thinning method indicated above to generate a trajectory of a DSPP from that of a pure stationary Poisson process P P S of density µ. The validity of this procedure was shown by [START_REF] Ogata | On Lewis simulation method for point processes[END_REF].

Two problems however immediately appear. The first one comes from the principle of this method. Indeed we have seen that one of its basic assumptions is that the density λ(t) of the PP P to be simulated must be bounded. When this density becomes random, [λ(t, ω)], it has in general no reason to be bounded, which makes impossible to satisfy the conditions of validity of the thinning procedure. For example if λ(t, ω) is the square of a normal signal, it is positive but not bounded. This is the reason why we shall study in the following statistical models of random density ensuring that it is bounded.

The second problem concerns the calculation of the coefficients R i used in the thinning procedure. We have seen above that the probability that R i = 1 is equal to λ(t i )/µ. For a given value of ω, or for a given trajectory of the random density of the PP to be simulated, we must then calculate λ(t i , ω).

The instants t i are the points in the time axis of a realization of a stationary

Poisson process P P S of density µ and in simulation procedure they are obtained from a realization of a sequence of IID positive exponential RVs. Consider now the calculation of λ(t i , ω). If the analytic expression of this function is known, which is often the case when λ(t) is not random, this can easily be done. On the other hand when λ(t) becomes random, the analytic expression of each trajectory is very rarely known and it is in general described only from some of its statistical properties.

In simulation experiments we must work with a density also obtained by simulation and this implies that in most of the cases this density is in practice only known at certain time instants τ i coming from the simulation procedure.

These instants have no relation with the t i s coming from the simulation of the PP P P S . Then, starting from times t i we must estimate λ(t i , ω), while we have at our disposal the values λ(τ i , ω). This requires the use of some approximations introducing a systematic bias in the procedure. Similar problems also appear with other simulations methods and some are discussed in (Giesecke et al., 2011). As seen below, this disadvantage completely disappears when using random densities taking only two possible values, which is a limit case of a piecewise-constant density used in [START_REF] Harrod | Numerical methods for realizing nonstationary Poisson processes with piecewise-constant instantaneous-rate functions[END_REF].

As we shall see later, it is possible to overcome without approximation the two previous difficulties of the thinning method for the simulation of DSPPs by using a binary density. Furthermore such a density is a rather good approximation of various physical situations. For example signals taking only two possible values are widely used in digital communications and information theory and in optical digital communications a point process with a binary intensity corresponds to various physical devices.

It is however interesting to note that DSPP are recently widely used in mathematical analysis of financial markets playing an essential role in the evolution of the world economic situation. A long list of papers appears in the reference list of [START_REF] Giesecke | Monte Carlo algorithms for default timing problems[END_REF]. As in statistical physics many studies require the use of simulation procedures and there are various attempts in order to simulate Poisson processes with random intensity. The two papers cited above present the most complete overview of the results obtained in this field.

Looking at this literature, the first point that clearly appears is that discrete valued intensities are not at all appropriate in this field for the description of random intensities of PPs. The fluctuations of the financial parameters studied are more conveniently described by diffusion or Brownian motion processes.

Furthermore the purpose of the simulation are rarely the properties of the simulated PP and more often some particular quantities especially interesting in financial problems such as loss or gain or evaluation of the risks. This is also the case of papers in the domain of insurance business. Even with this difference of perspective, the two problems appearing in the simulation of DSPPs and presented above are well noted and an attempt for their solution is introduced in the two papers of Giesecke cited in references. By using very abstract mathematical procedures of measure and martingale theory, Giesecke introduces a specific model of diffusion for the random intensity. It becomes then possible to use the thinning method for the simulation of DSPPs even for unbounded density, which is the case of diffusion processes, and to suppress the bias in the estimation of the sampling of the random intensity. The algorithms deduced from the theoretical analysis are appropriate for the domain of the paper but seem without interest in statistical physics, in such a way that the is almost no overlapping between these methods and those discussed below.

Doubly stochastic Poisson processes with binary density

Let us begin by showing the origin of the theoretical difficulties that appear in almost all the calculations with DSPPs, which justifies the necessity of a simulation approach. Consider for example the case of counting analysis which is one of the most common ways to study PPs. Let N be the random number of points of a PP in some interval which can arbitrarily be [0, T ]. If the PP is a Poisson process of density λ(t), the probability p n that N = n, where n is an arbitrary nonnegative integer, is p n = exp(-m)m n /n! where m = T 0 λ(θ)dθ. In the case of a DSPP, m becomes random, because of the randomness of λ(t, ω), in such a way that

p n = E exp(-M) M n n! with M = T 0 λ(θ, ω)dθ, (1)
where E means the expected value with respect to ω specifying the randomness of the density of the PP. As said on p. 171 of [START_REF] Saleh | Photoelectron Statistics[END_REF], where this problem is discussed, "To find the statistics of an integral over a stochastic process is, in general, a formidable task". This is effectively the origin of the difficulty of the calculations of the probabilities p n because of the lack of knowledge of the statistics of the random variable M appearing in (1).

Various attempts have been made to solve this problem. However the number of cases where an explicit solution can be found remains very limited. The first approach in the domain of statistical optics was presented in [START_REF] Bédard | Photon counting statistics of Gaussian light[END_REF] and an overview of some other cases is given in [START_REF] Saleh | Photoelectron Statistics[END_REF]. This requires in general a Karhunen Loeve expansion for which it is necessary to solve an integral equation and explicit analytical expressions are rarely obtained. Another simplification is often used in optics in the so-called long coherence time case. This means that the variations of λ(t, ω) are very small in intervals of duration of the order of T in such a way that in (1) the integral over [0, T ] can be replaced by λ(0, ω)T , which in a number of cases makes it possible to calculate of the probabilities p n .

Let us now present an example of random density in which the previous difficulties concerning the conditions of simulation by thinning or the mathematical calculations disappear. This density is defined by a general expression of the form

λ(t, ω) = λ 0 [1 + αB(t, ω)], (2)
where B(t, ω) is a random perfect clipped signal taking only the values ±1, the changes of signs appearing at the points of a PP P B assumed to be stationary.

It is first clear that, whatever the nature of P B , if |α| < 1, then the random density λ(t, ω) is positive and bounded, which ensures that the first condition of validity of the thinning method is satisfied.

Let us show that this is also the case for the second condition. Let S i denote the ordered sequence of the random points of P B and s i a particular realization of this sequence. It results from the definition of B(t) that in the interval s i < t < s i+1 it remains constant and equal to B(s + i ), the value of B(t) obtained immediately after the change of sign appearing at s i . As a result, if

B(s + 0
) is known we have B(t) = B(s + 0 )(-1) i for any t satisfying s i < t < s i+1 . From an algorithmic point of view we see that in order to calculate B(t) we must associate to any value of t the integer i(t) such that s

i(t) < t < s [i(t)+1] .
Applying the previous relation we obtain B(t) = B(s + 0 )(-1) i(t) , valid for any t. There are various ways to obtain i(t), the simplest one coming from the fact that i(t) is the only integer such that the product [t -s

i(t) ][t -s [i(t)+1] ] is
negative. This exact knowledge of the values of B(t) for any value of t allows one to use the thinning algorithm, according to the previous discussion.

It remains to show that in some cases this model introduces the possibility to make analytical calculations which in general are almost impossible. This is especially the case when the PP P B is a stationary Poisson process defined by its density ν.

In this specific case the signal B(t) introduced above becomes the random telegraph signal (RTS), widely used in random signal theory. It is introduced on p. 334 of [START_REF] Picinbono | Random Signals and Systems[END_REF] and studied as a particular case of ordered signals in [START_REF] Picinbono | Polyspectra of ordered signals[END_REF]. The main property is that all its moments of order higher than two have a simple mathematical expression, in such a way that, as shown in the next section, the calculation of various statistical properties of the integral appearing in (2) becomes possible. It is worth pointing out that the RTS is not only a simple mathematical extension of Poisson processes, but also a rather good statistical model to represent physical phenomena in which a system makes random transitions between two states, which is a very common situation [START_REF] Efros | Random telegraph signal in the photoluminescence intensity of a single quantum dot[END_REF]. It is also widely used in binary communications, which is the origin of its name, and densities such as (2), where B(t) is a RTS, can be used in a description of binary optical communications.

Statistics of the integral of the RTS

Let B(t) be the RTS of a stationary Poisson process P B of density ν and Z the integral

Z(t, T ) △ = t+T t B(θ)dθ. (3)
Because of the stationarity of P B the statistical properties of Z do not depend on t. The moment of order k of Z is given by

M Z k = ... D E[B(θ 1 )B(θ 2 )...B(θ k )]dθ 1 dθ 2 ...dθ k , ( 4
)
where D is a k-dimensional cube of edge length T . The expression of the high order moments of B(t) appearing in (4) is well known (see [START_REF] Picinbono | Polyspectra of ordered signals[END_REF].

The moments of odd order are zero and for even order we have

E[B(θ 1 )B(θ 2 )...B(θ 2k )] = exp[-2ν( θ2 -θ1 + θ4 -θ3 + ... + θ2k -θ2k-1 )], (5)
where { θi } is the ordered permutation of {θ i }, satisfying θj < θj+1 , ∀j, 1 ≤ j < 2k -1. After what are sometimes rather complicated calculations it is possible to deduce from ( 5) the values of the moments M Z k defined by (4). For example the moments of order 2 and 4 are given by

M Z 2 = (1/ν)[T -(1/2ν)(1 -e -2νT )] (6)
and

M Z 4 = (3/2ν 4 )[3 -4νT + 2(νT ) 2 -(3 + 2νT )e -2νT ]. (7)
Let us now see how the integral (3) can be calculated in a computer simulation experiment. Let ξ i be the distances between successive points in the Poisson process P B . As indicated above the instants s i of this PP are then defined by the recursion s i = s i-1 + ξ i . It results from the structure of the RTS that for any arbitrary positive integer N

Z(s i , s i+N -s i ) = s i+N s i B(θ)dθ = B(s + i ) N k=1 (-1) k-1 ξ i+k , (8)
where B(s + i ) is defined above. Then the integral of the RTS over an interval between two points of P B is reduced to a sum of N terms, where N -1 is the number of points of the PP between s i and s i+N . Consider now the calculation of Z(s i , T ). It is clear that s i + T is almost surely not a point of P B . There exist however simple algorithms that make it possible to calculate the number N i of points of P B appearing in the interval [s i , s i + T ] and they are deduced from the inequalities s i+N i < s i + T < s i+N i +1 . As a result we have

Z(s i , T ) = Z(s i , s i+N i -s i ) + s i +T s i+N i B(θ)dθ. (9)
Since s i+N i is the last point of P B before s i +T , there is no change of sign of B(θ)

in the interval defining the last integral of (9). Then it is equal to B(s s + i +N i )∆θ, where ∆θ = s i +T -s i+N i . As the ξ i s are the intervals between successive points of P B defined recursively by s i+1 = s i +ξ i+1 , we have s i+N i = s i +ξ i+1 +...+ξ i+N i , which allows us to calculate Z(s i , T ) uniquely from the sequence of M values ξ i , which are M realizations of IID positive exponential RVS. This the starting point of the statistical analysis that follows in which M is of the order of 10 6 .

In Table 1 we present results concerning the first four moments of Z(t, T ).

Their theoretical values are given above by ( 6) and ( 7) and furthermore M Z 1 = M Z 3 = 0. These moments are calculated and measured for two values of T and several values of the density ν of P B The theoretical and experimental results are noted M Z k and MZ k respectively. This table shows a very good agreement between the theoretical and experimental values. The precision of the method decreases when νT increases which comes from the structure of the algorithm of simulation. Indeed, as said above, for each i we have to calculate the numbers N i s defined above which are RVs of mean value νT . In order to reduce the duration of the calculation the N i greater than a certain limit are deleted. This clearly introduces a bias in the measurements which increases with νT . It can obviously be reduced but at the price of longer calculations.

Table 1. Moments of the Integral of the RTS. Let us now present some results concerning measurements of probability density functions (PDF). Our purpose is to determine by computer simulated experiments the PDF of the RV Z(t, T ). This PDF does not depend on t because of the assumption that P B is stationary. We can a priori expect that the distribution of Z contains a discrete part. Indeed if there is no point of P B in the interval [t, t+ T ] the possible values of the RV Z are ±T because there is no change of sign of the function B(θ) in the interval of integration. Since P B is a Poisson process, the probability that the interval [t, t + T ] does not contain points of P B is exp(-νT ). Then we have Pr(Z = ±T ) = exp(-νT ). This probability can also be obtained from computer simulation by using the same kind of experimental approach as previously and the results appear in Table 2. For two values of T we present the theoretical value of Pr(Z = ±T ) and the results of four independent experiments. The agreement between theory and experiment is quite good.

ν M Z 1 MZ 1 M Z 2 MZ 2 M Z 3 MZ 3 M Z 4 MZ 4 T = 1 0.
On the other hand there is no simple theoretical expression of the continuous part of the PDF p Z (z) of Z given by (3). This PDF can be measured and the results appear on Fig. 1 obtained for ν = T = 1. The experimental values of the PDF are located with a rather good precision on the curve -0.1z 2 +0.54, which does not have a simple explanation. The first step is to generate a sequence of N ordered points t i of the stationary Poisson process P P S of density µ. This is a straightforward standard task. Indeed the intervals between the successive points of this process are IID positive RVs with the exponential distribution µ exp(-µx) easily obtained from various computer programs.

The second step is quite similar and yields the sequence of ordered points s j of another stationary Poisson process P B of density ν. These points are the instants at which the RTS B(t) changes its value from ±1 to ∓1.

The third step is to associate to each t i of P P S the point s j(i) of P B such that the product (t i -s j(i) )(t i -s j(i)+1 ) is negative. This can be done in two ways. In the first, indicated above, for each t i we calculate the N products (t i -s j )(t is j+1 ) and j(i) is the only value of j for which this product is negative. This procedure however requires very long calculation times because there are of the order of N 2 products to calculate which yields in our experiments an order of 10 12 multiplications. With standard computers this can require more than several hundred hours of calculations. It is then more appropriate to introduce a recursive procedure which is much less time consuming. Its principle is as follows.

Suppose that s j(i) is known and consider the difference d n (i + 1) = t i+1s j(i)+n . It is clear that j(i + 1) = j(i) + n if and only if d n (i + 1) > 0 and d n+1 (i + 1) < 0. This yields that

j(i + 1) = j(i) + ∞ n=1 nF n (i), (10)
with

F n (i) = u[d n (i + 1)]v[d n+1 (i + 1)], (11)
where u(x) is the unit step function equal to 1 if x > 0 and to -1 otherwise, and v(x) = u(-x). The advantage of (10) seems weak, because the presence of a series of an infinite number of terms. It appears however that this number is in practice very small. In particular when the densities of the PPs P P S and P B satisfy ν << µ, which is a very common situation, a short calculation shows that there is practically only one term corresponding to n = 1 in the series of (10), and for each set of experiments we have introduced a procedure ensuring that replacing the series by a sum of a finite number of terms does not introduce an error in the calculation of (10).

The fourth step is to calculate the value of the thinning coefficient R i defined in the first paragraph of Section 2. It is a Bernoulli RV and the probability that the point t i be not erased, event corresponding to

R i = 1, is λ(t i )/µ.
This requires the knowledge of λ(t i ). It is the value of the density (2) for t = t i which is known when B(t i ) is known. For this we use the expression i) . Using the values of B(t i ), and then of λ(t i ), associated with each point t i of P P S , we can apply the thinning procedure and the non erased points from this procedure constitute a trajectory of the simulated DSPP to be used for the following experiments.

B(t) = B(s + 0 )(-1) i(t) indicated in Section 3. It results from the definition of s j(i) that B(t i ) = B(s + 0 )(-1) j(

Presentation of experimental results

Since it is very rarely possible to carry out theoretical calculations of statistical properties of DSPPs, our first task is to introduce computer experiments in which the results can be compared to those obtained from a mathematical analysis. This is especially the case for the study of the moments of the random number of points counted in a given time interval (counting experiments).

As indicated in various papers, it is necessary, except for Poisson processes, to make a clear distinction between relaxed and triggered counting experiments.

In the relaxed case the beginnings of the counting time intervals have no relation with the PP analyzed, while in the triggered case they depend on this PP, the simplest example being the one where the counting intervals are open by a point of the PP. Since the theoretical analysis is much more complicated for the triggered case, we shall only present results concerning relaxed counting experiments, which are also in general simpler to realize than the triggered ones.

Let us note N the number of random points of a stationary DSPP in an interval [t, t + T ]. It is a positive RV and the main purpose of counting experiments is to determine some properties of its probability distribution. Its theoretical calculation is in general very difficult, even for the particular model of random density introduced in this paper. On the other hand, we shall now verify that the results obtained in the previous sections make the theoretical analysis of the moments of the RV N much simpler.

The first four moments of a Poisson RV N of mean m are given by

E(N) = m, (12)
E(N 2 ) = m + m 2 , (13)
E(N 3 ) = m + 3m 2 + m 3 , (14)
E(N 4 ) = m + 7m 2 + 6m 3 + m 4 . ( 15
)
If N is now the number of points of a DSPP recorded in the arbitrary time interval introduced above (relaxed counting), its moments are simply obtained by replacing the quantities m k of these relations by E(M k ), where M is given by (1) in which 0 and T are replaced by t and t + T respectively. For the calculation of these expected values we start from (2) which introduces values of the integral of the RTS calculated above. By using the notations previously introduced we obtain

E(M) = λ 0 T, (16)
E(M 2 ) = λ 2 0 T 2 + α 2 M Z 2 , (17)
E(M 3 ) = λ 3 0 T 3 + 3α 2 T M Z 2 , (18)
E(M 4 ) = λ 4 0 T 4 + 6α 2 T 2 M Z 2 + α 4 M Z 4 , (19)
where M Z 2 and M Z 4 are given by ( 6) and ( 7) respectively. Replacing m k by E(M k ) in Eqs ( 12) to (15) yields the theoretical values m k = E(N k ) of the moments of order k of the random number of points of the DSPP considered.

These moments can also be measured by a procedure described in [START_REF] Picinbono | Algorithms for point processes analysis[END_REF] applied to a computer simulation of this PP obtained by the procedure introduced above. The experiment is realized with the parameters µ = 1, λ 0 = 1/2 and T = 2 in such a way that the mean value of the random number of points is 1. Furthermore various values of the parameter α defining the variance of the fluctuations of the random density of the PP are taken in consideration. For α = 0 the DSPP becomes a pure Poisson process. The results appear in Table 3. The first two lines correspond to the value α = 0, which means that the DSPP is a Poisson process. The first four moments are 1, 2, 5, and 15. The experimental values obtained are almost the same. For the other values taken by the parameter α there is also quite a good correspondence between the theoretical and experimental values of the moments. This means that the moment analysis of the relaxed number of points of the PP in a time interval correctly exhibits the properties of a DSPP defined by the random density (2).

By using the procedure described in [START_REF] Picinbono | Algorithms for point processes analysis[END_REF] it is also possible to measure the PDFs of the time intervals between successive points (lifetimes).

Unfortunately, even for this particular case of DSPP the theoretical calcula- that since λ 0 = 1/2 and since the mean value of B(t) is zero, then the density of the Poisson process obtained by this particular mode of thinning is also equal to 1/2, which implies that the PDF of the lifetime is exponential with the density p(x) = 0.5 exp(-0.5x). This clearly appears on the figure where the points are the results of experimental measurements while the continuous curve corresponds to this theoretical DDP. On the other hand, as indicated above, there is no theoretical result concerning the PDF of the lifetime when α = 0 and the figure presents experimental results obtained for a DSPP defined by α = 0.9. The only fact that exhibits this PDF is that there are two kinds of exponential behaviour corresponding approximately to the two values taken by the random density λ(t) given by (2). We shall now conduct a more precise analysis of this point in a case where some calculations become possible.

This case, already introduced above, is often used in the interpretation of experimental results of statistical optics and corresponds to the so-called long coherence time approximation. This means that the fluctuations of B(t) in (2) are sufficiently slow so that in integrals such as (3), B(t) remains constant and becomes a RV, which greatly simplifies the calculation of the statistical properties of this integral, which, without this assumption, is the major difficulty of the theoretical analysis. Looking at (2) we see that this approximation can be expressed in terms of the densities µ and ν of the PPs P P S and P B respectively by the relation µ >> ν. In the previous experiments these values were µ = 1 and ν = 0.2, which do not satisfy this inequality. The last experiments of this paper, on the other hand, are realized with the values µ = 1 and ν = 10 -3 , which does correspond much better to the previous inequality. From this approximation we can deduce that the random density is now given by (2) where B(t) is simply a RV taking the values ±1 with the same probabilities 1/2.

By using this approximation we can use a method introduced on p. 345 of [START_REF] Picinbono | Random Signals and Systems[END_REF] in order to calculate the PDF of two distinct time intervals appearing in the study of PPs. The residual lifetime is the random distance between an arbitrary point A and the first point of the process posterior to A, while the lifetime or order one is the distance between two successive points of the PP. In general, it is almost impossible to calculate the PDFs of these RVs noted p RL (x) and p L (x) respectively. But under the assumption introduced just above and by using the methods introduced in [START_REF] Picinbono | Random Signals and Systems[END_REF] we obtain p RL (x) = (1/2)[λ 1 exp(-λ 1 x) + λ 2 exp(-λ 2 x)], (20)

p L (x) = λ 2 1 exp(-λ 1 x) + λ 2 2 exp(-λ 2 x), (21)
where

λ 1 = λ 0 (1 + α), λ 2 = λ 0 (1 -α). ( 22
)
The results of experimental measurements of the PDFs of simulated data are presented in Fig. 3. For these experiments we have used 510 6 samples of exponential RVs for the generation of the PP P P S of density µ = 1. The density ν of the Poisson process generating the RTS B(t) was 10 -3 and finally α = 0.9.

The theoretical values give by ( 20) and ( 21) are represented by continuous curves while the results of PDF measurements appear as points in the figure.

The agreement between the theoretical analysis and experimental results is rather good. There is a small difference between them and we have verified that it results from the insufficient approximation of the calculation. The results are better if we take ν < 10 -3 , but the duration of the experiment also increases and obtaining greater precision would require more memory than we can obtain in our experiment. Nevertheless the main point that clearly appears from this experiment is the difference between the statistics of the lifetime and the residual lifetime, as predicted by the theoretical analysis.

-1 -0.8 -0.6 -0.

  tions of these PDFs are almost impossible and initially we shall only present the results of computer experiments. They appear in Fig.2. In this figure the points represent experimental measurements of the PDFs of the lifetimes obtained for two values of α. For α = 0 the DSPP becomes a Poisson Process, as already verified in other experiments carried out previously. It results from (2)

Figure 1 :

 1 Figure 1: Probability density functions p z (z). Points: experiment, continuous curves: -0.1z 2 + 0.54.

Figure 2 :Figure 3 :

 23 Figure 2: Probability density functions p(x) of lifetimes of a DSPP for various values of a. Points: experiment, continuous curves: theory.

Table 2 .

 2 Probability that Z of (3) is equal to ±T

	T = 1	T = 2
	Th.	0.3679	0.1353
	Exp.	0.3676	0.1353
		0.3682	0.1352
		0.3678	0.1358
		0.3681	0.1355
	5 Algorithms for the simulation of a DSPP
	with binary density

Following the definitions introduced in Section 3 we shall now present the sequences of algorithms used for the simulation of a DSPP with a binary density given by (2) where B(t) is the RTS defined from a stationary Poisson process P B .

Table 3 .

 3 First Four Moments of the Number of Points

		m 1	m 2	m 3	m 4
	α = 0
	Th.	1	2	5	15
	Exp.	0.99	1.99	4.99	15.04
	α = 0.5
	Th.	1.0000 2.1948 6.1687 21.0844
	Exp.	1.001	2.196	6.168	21.083
	α = 0.8
	Th.	1.0000 2.4987 7.9919 30.7596
	Exp.	1.004	2.504	7.995	30.735
	α = 0.9
	Th.	1.0000 2.6311 8.7867 35.0470
	Exp.	1.007	2.645	8.809	35.036

Unité mixte de recherche (UMR 8506) du Centre national de la recherche scientifique (CNRS), de l' École supérieure d'électricité (Suplec) et de l'Université de Paris-Sud 11 (UPS).


Mnacho Echenim

Nicolas Peltier

Sophie Tourret

❛❞❞❡❞

A Superposition Strategy for Abductive Reasoning in Ground Equational Logic

♦ t❤❛t ❛•• ♣r✐♠❡ ✐♠♣•✐❝❛t❡s ❝❛♥ ❜❡ ♦❜t❛✐♥❡❞✱ ❜✉t ♦♥ t❤❡ ♦t❤❡r ❤❛♥❞✱ t❤❡② s❤♦✉•❞ ❜❡ ♠♦r❡ ❢♦❝✉s❡❞ t❤❛♥ t❤❡ ❜•✐♥❞ ❛♣♣•✐❝❛t✐♦♥ ♦❢ t❤❡ r❡s♦•✉t✐♦♥ r✉•❡✱ ✐♥ ♦r❞❡r t♦ ②✐❡•❞ ❛ ♠♦r❡ ❡✣❝✐❡♥t ❛♥❞ ❤♦♣❡❢✉••② s❝❛•❛❜•❡ ❛•❣♦r✐t❤♠✳ ❲❡ ❞❡✈✐s❡ s✉❝❤ ❛ str❛t❡❣② ✐♥ t❤❡ ♣r❡s❡♥t ♣❛♣❡r✳ ✸ |❡♣r❡s❡♥t❛t✐♦♥ ♦❢ ❈•❛✉s❡s ▼♦❞✉•♦ ❊q✉❛•✐t② ■♥ ♣r♦♣♦s✐t✐♦♥❛• •♦❣✐❝✱ ❞❡t❡❝t✐♥❣ r❡❞✉♥❞❛♥t ❝•❛✉s❡s ✐s ❛♥ ❡❛s② t❛s❦✱ s✐♥❝❡ ❛ ✭♥♦♥✲ t❛✉t♦•♦❣✐❝❛•✮ ❝•❛✉s❡ C ✐s ❛ •♦❣✐❝❛• ❝♦♥s❡q✉❡♥❝❡ ♦❢ D ✐✛ D ✐s ❛ s✉❜❝•❛✉s❡ ♦❢ C✳ ❚❤✉s C ✐s r❡❞✉♥❞❛♥t ✐♥ ❛ ❝•❛✉s❡ s❡t ✐✛ t❤❡r❡ ❡①✐sts ❛ ❝•❛✉s❡ D ∈ S s✉❝❤ t❤❛t D ⊆ C✳ ❚❤❡ •✐t❡r❛•s ❝❛♥ ❜❡ t♦t❛••② ♦r❞❡r❡❞ t♦ ❤❛♥❞•❡ ❝♦♠♠✉t❛t✐✈✐t② ❛♥❞ t❤❡ ❝•❛✉s❡ s❡t S ❝❛♥ ❜❡ r❡♣r❡s❡♥t❡❞ ❛s ❛ tr✐❡ ✭❛ tr❡❡✲❜❛s❡❞ ❞❛t❛✲str✉❝t✉r❡ ❝♦♠♠♦♥•② ✉s❡❞ t♦ r❡♣r❡s❡♥t str✐♥❣s ❬✹(✮✱ s♦ t❤❛t ✐♥❝•✉s✐♦♥ ❝❛♥ ❜❡ t❡st❡❞ ❡✣❝✐❡♥t•②✳ ❍♦✇❡✈❡r✱ ✐♥ ❣r♦✉♥❞ ❡q✉❛t✐♦♥❛• •♦❣✐❝✱ t❤❡ ❛❜♦✈❡ ♣r♦♣❡rt② ❞♦❡s ♥♦t ❤♦•❞ ❛♥②♠♦r❡✿ ❢♦r ✐♥st❛♥❝❡ t❤❡ ❝•❛✉s❡ a ≃ b∨b ≃ c ✐s ❛ •♦❣✐❝❛• ❝♦♥s❡q✉❡♥❝❡ ♦❢ a ≃ c ❜✉t ♦❜✈✐♦✉s•② a ≃ c ✐s ♥♦t ❛ s✉❜❝•❛✉s❡ ♦❢ a ≃ b ∨ b ≃ c✳ ❚❤✉s t❡st✐♥❣ ❝•❛✉s❡ ✐♥❝•✉s✐♦♥ ✐s ♥♦ •♦♥❣❡r s✉✣❝✐❡♥t ❛♥❞ r❡♣r❡s❡♥t✐♥❣ ❝•❛✉s❡ s❡ts ❛s tr✐❡s ✇♦✉•❞ ②✐❡•❞ ♠❛♥② ✉♥❞❡s✐r❡❞ r❡❞✉♥❞❛♥❝✐❡s✿ ❢♦r ✐♥st❛♥❝❡ t❤❡ ❝•❛✉s❡s a ≃ b ∨ b ≃ c ❛♥❞ a ≃ b ∨ a ≃ c ✇♦✉•❞ ❜❡ ❜♦t❤ st♦r❡❞✱ ❛•t❤♦✉❣❤ t❤❡② ❛r❡ ❡q✉✐✈❛•❡♥t✳ ❖✉r ✜rst t❛s❦ ✐s t❤✉s t♦ ❞❡✈✐s❡ ❛ ♥❡✇ r❡❞✉♥❞❛♥❝② ❝r✐t❡r✐♦♥ t❤❛t ❣❡♥❡r❛•✐③❡s s✉❜s✉♠♣t✐♦♥✱ t♦❣❡t❤❡r ✇✐t❤ ❛ ♥❡✇ ✇❛② ♦❢ r❡♣r❡s❡♥t✐♥❣ ❝•❛✉s❡s✱ t❤❛t t❛❦❡s ✐♥t♦ ❛❝❝♦✉♥t t❤❡ s♣❡❝✐❛• ♣r♦♣❡rt✐❡s ♦❢ t❤❡ ❡q✉❛•✐t② ♣r❡❞✐❝❛t❡✳ ❚♦ t❤✐s ♣✉r♣♦s❡ ✇❡ s❤♦✇ ❤♦✇ t♦ ♥♦r♠❛•✐③❡ ❣r♦✉♥❞ ❝•❛✉s❡s ❛❝❝♦r❞✐♥❣ t♦ ❛ t♦t❛• ♦r❞❡r✐♥❣ ♦♥ ❝♦♥st❛♥t s②♠❜♦•s✱ ❛♥❞ ✇❡ ✐♥tr♦❞✉❝❡ ❛ ♥❡✇ ♥♦t✐♦♥ ♦❢ ♣r♦•❡❝t✐♦♥✳ ▲❡t ≺ ❜❡ ❛ ✜①❡❞✱ t♦t❛• ♦r❞❡r✐♥❣ ♦♥ ❝♦♥st❛♥t s②♠❜♦•s ❛♥❞ C ❜❡ ❛ ❝•❛✉s❡✳ ❚❤❡ C✲r❡♣r❡s❡♥t❛t✐✈❡s ♦❢ ❛ ❝♦♥st❛♥t a✱ ❛ •✐t❡r❛• ✷ l ❛♥❞ ❛ ❝•❛✉s❡ D ❛r❡ r❡s♣❡❝t✐✈❡•② ❞❡✜♥❡❞ ❛s ❢♦••♦✇s✿ a ⇂C = min ≺ {b | b ≃ a |= C}; l ⇂C = a ⇂C ⊲⊳ b ⇂C ✇❤❡♥ l = a ⊲⊳ b; ❛♥❞ D ⇂C = {l ⇂C | l ∈ D}✳ ❚❤❡ •✐t❡r❛• l ⇂C ✭r❡s♣✳ t❤❡ ❝•❛✉s❡ D ⇂C ✮ ✐s ❝❛••❡❞ t❤❡ ♣r♦•❡❝t✐♦♥ ♦❢ l ✭r❡s♣✳ D✮ ♦♥ C✳ ■t ❝❛♥ ❜❡ ✈❡r✐✜❡❞ t❤❛t ❡✈❡r② ❝•❛✉s❡ C ✐s ❡q✉✐✈❛•❡♥t t♦ t❤❡ ❝•❛✉s❡✿ C ↓ ❞❡❢ = a∈C,a =a ⇂C a ≃ a ⇂C ∨ a≃b∈C a ⇂C ≃ b ⇂C ✳ ❚❤❡ ❝•❛✉s❡ C ↓ ✐s t❤❡ ♥♦r♠❛• ❢♦r♠ ♦❢ C ❛♥❞ ✐t ✐s ✉♥✐q✉❡✳ ❈•❛✉s❡s ✐♥ ♥♦r♠❛• ❢♦r♠ ❛r❡ st♦r❡❞ ✐♥ ❛ tr❡❡ ❞❛t❛✲str✉❝t✉r❡✱ ❝❛••❡❞ ❛ ❝•❛✉s❛• tr❡❡✱ s♣❡❝✐✜❝❛••② t❛✐•♦r❡❞ t♦ st♦r❡ s❡ts ♦❢ •✐t❡r❛•s ✇❤✐•❡ t❛❦✐♥❣ ✐♥t♦ ❛❝✲ ❝♦✉♥t t❤❡ ✉s✉❛• ♣r♦♣❡rt✐❡s ♦❢ t❤❡ ❡q✉❛•✐t② ♣r❡❞✐❝❛t❡✳ ❆s ✐♥ tr✐❡s✱ t❤❡ ❡❞❣❡s ♦❢ t❤❡ tr❡❡ ❛r❡ •❛❜❡•❡❞ ❜② •✐t❡r❛•s ❛♥❞ t❤❡ •❡❛✈❡s ❛r❡ •❛❜❡•❡❞ ❡✐t❤❡r ❜② ✭r❡♣r❡s❡♥t✐♥❣ t❤❡ ❡♠♣t② ❝•❛✉s❡✮ ♦r ❜② ∅ ✭❢❛✐•✉r❡ ♥♦❞❡✮✳ ❚❤❡ s❡t ♦❢ ❝•❛✉s❡s C(T ) r❡♣r❡s❡♥t❡❞ ❜② ❛ tr❡❡ T ✐s t❤❡ ❞✐s•✉♥❝t✐♦♥ ♦❢ •✐t❡r❛•s •❛❜❡•✐♥❣ ❛ ♣❛t❤ ❢r♦♠ t❤❡ r♦♦t t♦ t❤♦s❡ •❡❛✈❡s •❛❜❡•❡❞ ❜② ✳ ❋❛✐•✉r❡ ♥♦❞❡s ❛r❡ ✉s❡❢✉• ♠❛✐♥•② t♦ r❡♣r❡s❡♥t ❡♠♣t② s❡ts ✕ ✐♥ ❢❛❝t t❤❡② ❝❛♥ ❛•✇❛②s ❜❡ ❡•✐♠✐♥❛t❡❞ ❜② str❛✐❣❤t❢♦r✇❛r❞ s✐♠♣•✐✜❝❛t✐♦♥ r✉•❡s✱ ❡①❝❡♣t ✷ ❲❤❡♥ ✇❡ ✇r✐t❡ l = a ⊲⊳ b✱ t❤❡ s②♠❜♦• ⊲⊳ st❛♥❞s ❢♦r ≃ ♦r ≃✳ ✐❢ t❤❡ r♦♦t ✐ts❡•❢ ✐s •❛❜❡•❡❞ ❜② ∅✳ ❆ ❢❡✇ ❝♦♥str❛✐♥ts ❛r❡ ✐♠♣♦s❡❞ ♦♥ t❤❡ •✐t❡r❛•s •❛❜❡•✐♥❣ t❤❡ ❡❞❣❡s ♦❢ ❛ ❝•❛✉s❛• tr❡❡ ✕ t❤❡ ♠❛✐♥ ♦♥❡ ❜❡✐♥❣ t❤❛t ♥❡❣❛t✐✈❡ •✐t❡r❛•s ♦❝❝✉r ❜❡❢♦r❡ ♣♦s✐t✐✈❡ ♦♥❡s ✐♥ ❛•• ❜r❛♥❝❤❡s ♦❢ t❤❡ tr❡❡ ✕ ✐♥ ♦r❞❡r t♦ ❣✉❛r❛♥t❡❡ t❤❛t t❤❡ ❝•❛✉s❡s r❡♣r❡s❡♥t❡❞ ❜② ❛ ❝•❛✉s❛• tr❡❡ ❛r❡ ✐♥ ♥♦r♠❛• ❢♦r♠✱ ❛♥❞ t❤❛t t❤✐s ❞❛t❛✲ str✉❝t✉r❡ ❝❛♥ ❜❡ ✉s❡❞ t♦ t❡st r❡❞✉♥❞❛♥❝② ❡✣❝✐❡♥t•②✳ ❚❤❡ ❝♦♥❞✐t✐♦♥s ♣❡r♠✐t t♦ ❞❡s✐❣♥ ❛♥ ❡✣❝✐❡♥t t❡st ✇❤❡t❤❡r ❛ ❣✐✈❡♥ ❝•❛✉s❡ ✐s r❡❞✉♥❞❛♥t ✇✳r✳t✳ t❤♦s❡ st♦r❡❞ ✐♥ t❤❡ tr❡❡ ✭❢♦r✇❛r❞ s✉❜s✉♠♣t✐♦♥✮ ❛♥❞ ❝♦♥✈❡rs❡•② t♦ ❞❡•❡t❡ ❢r♦♠ t❤❡ tr❡❡ ❛•• ❝•❛✉s❡s t❤❛t ❛r❡ r❡❞✉♥❞❛♥t ✇✳r✳t✳ ❛ ♥❡✇•② ❣❡♥❡r❛t❡❞ ❝•❛✉s❡ ✭❜❛❝❦✇❛r❞ s✉❜s✉♠♣t✐♦♥✮✳ ❚❤❡♦r❡♠ ✶✳ ▲❡t C ❛♥❞ D ❜❡ t✇♦ ❝•❛✉s❡s✳ ❆ss✉♠❡ t❤❛t C ✐s ♥♦t ❛ t❛✉t♦•♦❣②✳ D |= C ✐✛ ❡✈❡r② ♥❡❣❛t✐✈❡ •✐t❡r❛• ✐♥ D ⇂C ✐s ❛ ❝♦♥tr❛❞✐❝t✐♦♥ ❛♥❞ ❡✈❡r② ♣♦s✐t✐✈❡ •✐t❡r❛• ✐♥ D ⇂C ✐s ❛•s♦ ✐♥ C ⇂C ✳ ❲❡ ✇✐•• ✉s❡ ❛ s•✐❣❤t•② ♠♦r❡ r❡str✐❝t✐✈❡ ✈❡rs✐♦♥ ♦❢ t❤✐s ❝r✐t❡r✐♦♥ ❢♦r r❡❞✉♥❞❛♥❝② ❡•✐♠✐♥❛t✐♦♥✱ ❜② ✈✐❡✇✐♥❣ t❤❡ ❝•❛✉s❡s ❛s ♠✉•t✐s❡ts ❛♥❞ ❜② ✐♠♣♦s✐♥❣ t❤❛t t❤❡ ♣♦s✐t✐✈❡ •✐t❡r❛•s ✐♥ D ⇂C ❛r❡ ♠❛♣♣❡❞ t♦ ♣❛✐r✇✐s❡ ❞✐st✐♥❝t •✐t❡r❛•s ✐♥ C ⇂C ✳ ❚❤✐s ❛❞❞✐t✐♦♥❛• r❡str✐❝t✐♦♥ ✐s ♥❡❝❡ss❛r② t♦ ♣r❡✈❡♥t t❤❡ ❢❛❝t♦rs ♦❢ ❛ ❝•❛✉s❡ ❢r♦♠ ❜❡✐♥❣ r❡❞✉♥❞❛♥t ✇✳r✳t✳ t❤❡ ✐♥✐t✐❛• ❝•❛✉s❡✳ ❋♦r ✐♥st❛♥❝❡✱ t❤❡ ❝•❛✉s❡ a ≃ b ∨ a ≃ a ′ ✇✐•• ♥♦t ❜❡ r❡❞✉♥❞❛♥t ✇✳r✳t✳ a ≃ b ∨ a ′ ≃ b ∨ a ≃ a ′ ✱ ❛•t❤♦✉❣❤ a ≃ b ∨ a ′ ≃ b ∨ a ≃ a ′ |= a ≃ b ∨ a ≃ a ′ ✳ ❆•❣♦r✐t❤♠ ✶ ✐s❊♥t❛✐•❡❞(C, T ) |❡q✉✐r❡✿ C ✐s ❛ ❝•❛✉s❡ ✐♥ ♥♦r♠❛• ❢♦r♠ ❛♥❞ T ✐s ❛♥ RN ✲❝•❛✉s❛• tr❡❡✳ ❊♥s✉r❡✿ ✐s❊♥t❛✐•❡❞(C, T ) = true ⇔ ∃D ∈ C(T ), D |= C ✐❢ T = t❤❡♥ r❡t✉r♥ true ❡♥❞ ✐❢ ✐❢ C = t❤❡♥ r❡t✉r♥ false ❡♥❞ ✐❢ l1 ← min < {l ∈ C} ❢♦r ❛•• (l, T ′ ) ∈ T s✉❝❤ t❤❛t l ≥ l1 ❞♦ ✐❢ l1 = a ≃ b✱ ✇✐t❤ a ≻ b t❤❡♥ ✐❢ l = l1 t❤❡♥ ✐❢ ✐s❊♥t❛✐•❡❞(C\l1, T ′ ) t❤❡♥ r❡t✉r♥ true ❡♥❞ ✐❢ ❡•s❡ ✐❢ ¬(l = a ≃ c)✱ ✇✐t❤ a ≻ c t❤❡♥ ✐❢ ✐s❊♥t❛✐•❡❞(C\l1, (T ′ | ↑ l )[a := b]) t❤❡♥ r❡t✉r♥ true ❡♥❞ ✐❢ ❡♥❞ ✐❢ ❡•s❡ ✐❢ l ∈ C t❤❡♥ ✐❢ ✐s❊♥t❛✐•❡❞(C\l, T ′ ) t❤❡♥ r❡t✉r♥ true ❡♥❞ ✐❢ ❡♥❞ ✐❢ ❡♥❞ ❢♦r r❡t✉r♥ false ❚❤❡ ✜rst ❛•❣♦r✐t❤♠ ✐s ✐♥✈♦❦❡❞ ♦♥ ❛ ❝•❛✉s❡ C ❛♥❞ ❛ tr❡❡ T ✱ ❛♥❞ r❡t✉r♥s true ✐❢ ❛♥❞ ♦♥•② ✐❢ t❤❡r❡ ❡①✐sts ❛ ❝•❛✉s❡ D ✐♥ C(T ) s✉❝❤ t❤❛t C ✐s r❡❞✉♥❞❛♥t ✇✳r✳t✳ D✳ ❚♦ t❡st t❤✐s ❡♥t❛✐•♠❡♥t✱ t❤❡ ❛•❣♦r✐t❤♠ ♣❡r❢♦r♠s ❛ ❞❡♣t❤✲✜rst tr❛✈❡rs❛• ♦❢ T ❛♥❞ ❛tt❡♠♣ts t♦ ♣r♦•❡❝t ❡✈❡r② ❡♥❝♦✉♥t❡r❡❞ •✐t❡r❛• ♦♥ C✳ ■❢ ❛ •✐t❡r❛• ❝❛♥♥♦t ❜❡ ♣r♦•❡❝t❡❞✱ t❤❡ ❡①♣•♦r❛t✐♦♥ ♦❢ t❤❡ s✉❜tr❡❡ ❛ss♦❝✐❛t❡❞ t♦ t❤✐s •✐t❡r❛• ✐s ✉s❡•❡ss✱ s♦ t❤❡ ❛•❣♦r✐t❤♠ s✇✐t❝❤❡s t♦ t❤❡ ❢♦••♦✇✐♥❣ •✐t❡r❛•✳ ❆s s♦♦♥ ❛s ❛ ❝•❛✉s❡ ❡♥t❛✐•✐♥❣ C ✐s ❢♦✉♥❞✱ t❤❡ tr❛✈❡rs❛• ❤❛•ts ❛♥❞ true ✐s r❡t✉r♥❡❞✳ ❚❤❡ s❡❝♦♥❞ ❛•❣♦r✐t❤♠ ❞❡•❡t❡s ❢r♦♠ ❛ tr❡❡ T ❛•• ❝•❛✉s❡s ❡♥t❛✐•❡❞ ❜② C✳ ■t ♣❡r❢♦r♠s ❛ ❞❡♣t❤✲✜rst tr❛✈❡rs❛• ♦❢ T ❛♥❞ ❛tt❡♠♣ts t♦ ♣r♦•❡❝t C ♦♥ ❡✈❡r② ❝•❛✉s❡ ✐♥ C(T )✱ ❞❡•❡t✐♥❣ t❤♦s❡ ♦♥ ✇❤✐❝❤ s✉❝❤ ❛ ♣r♦•❡❝t✐♦♥ s✉❝❝❡❡❞s✳ ❆s s♦♦♥ ❛s ❛ ♣r♦•❡❝t✐♦♥ ✐s ✐❞❡♥t✐✜❡❞ ❛s ✐♠♣♦ss✐❜•❡✱ t❤❡ ❡①♣•♦r❛t✐♦♥ ♦❢ t❤❡ ❛ss♦❝✐❛t❡❞ s✉❜tr❡❡ ❤❛•ts ❛♥❞ t❤❡ ❛•❣♦r✐t❤♠ ♠♦✈❡s ♦♥ t♦ t❤❡ ♥❡①t ❝•❛✉s❡✳ ❲❤❡♥ ❡✈❡r② •✐t❡r❛• ✐♥ C ❤❛s ❜❡❡♥ ♣r♦•❡❝t❡❞✱ ❛•• t❤❡ ❝•❛✉s❡s r❡♣r❡s❡♥t❡❞ ✐♥ t❤❡ ❝✉rr❡♥t s✉❜tr❡❡ ❛r❡ ❡♥t❛✐•❡❞ ❜② C✱ ❛♥❞ ❛r❡ t❤❡r❡❢♦r❡ ❞❡•❡t❡❞✳ ❆•❣♦r✐t❤♠ ✷ ♣r✉♥❡❊♥t❛✐•❡❞(C, T ) |❡q✉✐r❡✿ C ✐s ❛ ❝•❛✉s❡ ✐♥ r❡•❛①❡❞ ♥♦r♠❛• ❢♦r♠✱ T ✐s ❛♥ N ✲❝•❛✉s❛• tr❡❡ ❛♥❞ ✐s❊♥t❛✐•❡❞(C, T ) = false ❊♥s✉r❡✿ ∀D ∈ C(T )✱ C |= D ✐❢ C = t❤❡♥ T ← ∅ ❡①✐t ❡♥❞ ✐❢ ✐❢ T = t❤❡♥ ❡①✐t ❡♥❞ ✐❢ l1 ← min < {li ∈ C} ❢♦r ❛•• (l, T ′ ) ∈ T s✉❝❤ t❤❛t l ≤ l1 ❞♦ ✐❢ l1 = l t❤❡♥ ♣r✉♥❡❊♥t❛✐•❡❞(C\l1, T ′ ) ❡•s❡ ✐❢ l = a ≃ b t❤❡♥ ♣r✉♥❡❊♥t❛✐•❡❞(C, T ′ ) ❡•s❡ ✐❢ l = a ≃ b✱ ✇✐t❤ a ≻ b ❛♥❞ ∄c✱ l1 = a ≃ c✱ ✇✐t❤ a ≻ c t❤❡♥ ♣r✉♥❡❊♥t❛✐•❡❞(C[a := b], T ′ ) ❡♥❞ ✐❢ ❡♥❞ ✐❢ ❡♥❞ ❢♦r ✹ •❡♥❡r❛t✐♦♥ ♦❢ ■♠♣•✐❝❛t❡s ❚❤❡ ❝•❛✉s❛• tr❡❡s ❞❡s❝r✐❜❡❞ ❛❜♦✈❡ ❛r❡ ✉s❡❞ ❞✉r✐♥❣ t❤❡ ❣❡♥❡r❛t✐♦♥ ♦❢ ✐♠♣•✐❝❛t❡s t♦ st♦r❡ t❤❡ ♥♦♥✲r❡❞✉♥❞❛♥t ❝•❛✉s❡s ❣❡♥❡r❛t❡❞ ❜② t❤❡ ❢♦••♦✇✐♥❣ ✐♥❢❡r❡♥❝❡ r✉•❡s✿ a ≃ b ∨ C a ′ ≃ c ∨ D a ≃ a ′ ∨ b ≃ c ∨ C ∨ D a ≃ b ∨ a ′ ≃ b ′ ∨ C a ≃ b ∨ a ≃ a ′ ∨ b ≃ b ′ ∨ C n i=1 (a i ≃ b i ) ∨ P 1 c ≃ d ∨ P 2 n i=1 (a i ≃ c ∨ d ≃ b i ) ∨ P 1 ∨ P 2 ❚❤❡ ♣r❡♠✐s❡s ❛r❡ ❛ss✉♠❡❞ t♦ ❜❡ ✐♥ ♥♦r♠❛• ❢♦r♠ ❛♥❞ t❤❡ ❝♦♥❝•✉s✐♦♥ ✐s ♥♦r♠❛•✲ ✐③❡❞ ❜❡❢♦r❡ ❜❡✐♥❣ st♦r❡❞✳ ❚❤❡ ✜rst r✉•❡ ✐s s✐♠✐•❛r t♦ t❤❡ ✉s✉❛• ♣❛r❛♠♦❞✉•❛t✐♦♥ r✉•❡✱ ❡①❝❡♣t t❤❛t t❤❡ ✉♥✐✜❝❛t✐♦♥ ❜❡t✇❡❡♥ t❤❡ t❡r♠s a ❛♥❞ a ′ ✐s ♦♠✐tt❡❞ ❛♥❞ r❡♣•❛❝❡❞ ❜② t❤❡ ❛❞❞✐t✐♦♥ ♦❢ t❤❡ •✐t❡r❛• a ≃ a ′ ❡♥s✉r✐♥❣ t❤❛t t❤❡s❡ t❡r♠s ❛r❡ s❡♠❛♥t✐❝❛••② ❡q✉❛•✳ ❙✐♠✐•❛r•②✱ t❤❡ s❡❝♦♥❞ r✉•❡ ❢❛❝t♦r✐③❡s t❤❡ •✐t❡r❛•s a ≃ b ❛♥❞ a ′ ≃ b ′ ❛♥❞ ❛❞❞s t❤❡ ❝♦♥❞✐t✐♦♥s a ≃ a ′ ❛♥❞ b ≃ b ′ ✳ ❚❤❡ •❛st r✉•❡ ❝♦rr❡s♣♦♥❞s t♦ ❛♥ ❛♣♣•✐❝❛t✐♦♥ ♦❢ t❤❡ ❢❛❝t♦r✐③❛t✐♦♥ r✉•❡ ♦♥ t❤❡ ♥❡❣❛t✐✈❡ •✐t❡r❛•s a i ≃ b i ✱ ❢♦•✲ •♦✇❡❞ ❜② ❛ ♣❛r❛♠♦❞✉•❛t✐♦♥ st❡♣ ✇❤✐❝❤ r❡♠♦✈❡s t❤❡s❡ •✐t❡r❛•s✱ ✇❤✐•❡ ❛❞❞✐♥❣ t❤❡ ❝♦♥❞✐t✐♦♥s ❡♥s✉r✐♥❣ t❤❛t a i = c ❛♥❞ b i = d✳ ◆♦t❡ t❤❛t ♥♦ ♦r❞❡r✐♥❣ ❝♦♥❞✐t✐♦♥ ✐s ✐♠♣♦s❡❞✳ ■t ✐s ❡❛s② t♦ s❡❡ t❤❛t t❤❡ ❛❞❞✐t✐♦♥ ♦❢ ♦r❞❡r✐♥❣ r❡str✐❝t✐♦♥s ✭❛s ✐♥ t❤❡ ✉s✉❛• ✈❡rs✐♦♥ ♦❢ t❤❡ r✉•❡s✮ ♠❛❦❡s t❤❡ ❝❛•❝✉•✉s ✐♥❝♦♠♣•❡t❡✳ ❋♦r ✐♥st❛♥❝❡✱ ❛ss✉♠❡ t❤❛t S = {a ≃ b ∨ c ≃ d, a ≃ b}✱ ✇✐t❤ a ≺ b ≺ c ≺ d✳ ❚❤❡♥ S |= c ≃ d✱ ❜✉t a ≃ b ✐s ♥♦t ♠❛①✐♠❛• ✐♥ a ≃ b ∨ c ≃ d✱ t❤✉s ✐t ✐s ✐♠♣♦ss✐❜•❡ t♦ ✐♥❢❡r t❤❡ ❝•❛✉s❡ c ≃ d ✐❢ ♣❛r❛♠♦❞✉•❛t✐♦♥ ✐s r❡str✐❝t❡❞ t♦ ♠❛①✐♠❛• •✐t❡r❛•s✳ ❙✐♠✐•❛r•②✱ ✐❢ S ′ = {a ≃ b, c ≃ b}✱ t❤❡♥ S ′ |= a ≃ c✱ ❜✉t a ≃ c ❝❛♥ ♦♥•② ❜❡ ❞❡r✐✈❡❞ ✐❢ t❤❡ ♣❛r❛♠♦❞✉•❛t✐♦♥ r✉•❡ ✐s ❛♣♣•✐❡❞ ♦♥ t❤❡ ♥♦♥✲♠❛①✐♠❛• t❡r♠ a ✐♥ a ≃ b ✭r❡♣•❛❝✐♥❣ a ❜② b ✐♥ a ≃ b ②✐❡•❞s ❛ ❝♦♥tr❛❞✐❝t✐♦♥ ❛♥❞ t❤❡ ❝♦♥❞✐t✐♦♥ a ≃ c ✐s ❛❞❞❡❞ ✐♥ t❤❡ ❝•❛✉s❡✮ ♦r ❢r♦♠ t❤❡ ♥♦♥✲♠❛①✐♠❛• t❡r♠ b ✐♥ c ≃ b ✭r❡♣•❛❝✐♥❣ b ❜② c ✐♥ a ≃ b ②✐❡•❞s a ≃ c✮✳ ❙♦♠❡ ♦❢ t❤❡ •✐t❡r❛•s ✐♥ t❤❡ ❝♦♥❝•✉s✐♦♥s✱ ♠♦r❡ ♣r❡❝✐s❡•② t❤❡ ❝♦♥❞✐t✐♦♥s a ≃ a ′ ✱ b ≃ b ′ ✱ a i ≃ c ❛♥❞ d ≃ b i ❝❛♥ ❜❡ ✏❢r♦③❡♥✑ ✐♥ t❤❡ s❡♥s❡ t❤❛t ♥♦ ❢✉rt❤❡r ✐♥❢❡r❡♥❝❡ ✐s ❛••♦✇❡❞ ✇✐t❤✐♥ t❤❡♠ ✭t❤❡s❡ •✐t❡r❛•s ✇✐•• ❡✈❡♥t✉❛••② r❡♠❛✐♥ ✕ ❛❢t❡r ♥♦r♠❛•✐③❛t✐♦♥ ✕ ✐♥ t❤❡ ❝♦♥s✐❞❡r❡❞ ♣r✐♠❡ ✐♠♣•✐❝❛t❡✮✳ ❆•t❤♦✉❣❤ t❤✐s r❡♠❛r❦ ❝❛♥ ❞✐s♠✐ss ♠❛♥② ✐♥❢❡r❡♥❝❡s✱ ✐ts ♣r❛❝t✐❝❛• ✐♥t❡r❡st r❡♠❛✐♥s ✉♥❝•❡❛r✱ s✐♥❝❡ t❤❡ ❢r♦③❡♥ •✐t❡r❛•s ❤❛✈❡ t♦ ❜❡ ❝♦♥s✐❞❡r❡❞ ❛♣❛rt ✇❤❡♥ ❛♣♣•②✐♥❣ t❤❡ r❡❞✉♥❞❛♥❝② ❞❡t❡❝t✐♦♥✱ ✇❤✐❝❤ ♠❛② ♣r❡✈❡♥t t❤❡ r❡♠♦✈❛• ♦❢ ♥✉♠❡r♦✉s ❝•❛✉s❡s✳ ❆ s❡t ♦❢ ❝•❛✉s❡s S ✐s s❛t✉r❛t❡❞ ✐✛ ❢♦r ❡✈❡r② ❝•❛✉s❡ C t❤❛t ❝❛♥ ❜❡ ❞❡r✐✈❡❞ ❢r♦♠ S ✉s✐♥❣ t❤❡s❡ r✉•❡s✱ t❤❡r❡ ❡①✐sts ❛ ❝•❛✉s❡ C ′ ∈ S s✉❝❤ t❤❛t C ✐s r❡❞✉♥❞❛♥t ✇✳r✳t✳ C ′ ✳ ❚❤❡ ❢♦••♦✇✐♥❣ t❤❡♦r❡♠ st❛t❡s t❤❡ ❝♦♠♣•❡t❡♥❡ss ♦❢ t❤✐s ❛♣♣r♦❛❝❤✿ ❚❤❡♦r❡♠ ✷✳ ▲❡t S ❜❡ ❛ s❡t ♦❢ ❣r♦✉♥❞ ✢❛t ❝•❛✉s❡s ✐♥ ♥♦r♠❛• ❢♦r♠✳ ■❢ S ✐s s❛t✲ ✉r❛t❡❞ ❛♥❞ S |= C t❤❡♥ t❤❡r❡ ❡①✐sts ❛ ❝•❛✉s❡ C ′ ∈ S s✉❝❤ t❤❛t C ′ |= C✳ ❆ s❡❡♠✐♥❣•② ♥❛t✉r❛• ✐❞❡❛ ✇♦✉•❞ ❜❡ t♦ r❡str✐❝t t❤❡ •❛st r✉•❡ t♦ t❤❡ ❝❛s❡ n = 1 ✭✐✳❡✳ t♦ r❡♠♦✈❡ t❤❡ ♣r❡•✐♠✐♥❛r② ✐♠♣•✐❝✐t ❢❛❝t♦r✐③❛t✐♦♥ st❡♣ ✇❤❡♥ ♣❛r❛♠♦❞✉•❛t✐♥❣ ✐♥t♦ ♥❡❣❛t✐✈❡ •✐t❡r❛•s✮✳ ❍♦✇❡✈❡r✱ t❤✐s ♠❛❦❡s t❤❡ ❝❛•❝✉•✉s ✐♥❝♦♠♣•❡t❡✳ ❋♦r ✐♥st❛♥❝❡✱ ❝♦♥s✐❞❡r t❤❡ s❡t✿ S = {c ≃ a ∨ d ≃ a, a ≃ b ∨ c ≃ b ∨ d ≃ b}✳ ■t ✐s ❝•❡❛r t❤❛t S |= c ≃ b ∨ d ≃ b ❛♥❞ t❤❛t S ✐s s❛t✉r❛t❡❞ ✐❢ t❤❡ ♥❡❣❛t✐✈❡ ♣❛r❛♠♦❞✉•❛t✐♦♥ r✉•❡ ✐s ❛♣♣•✐❡❞ ♦♥ ❛ ✉♥✐q✉❡ •✐t❡r❛• ♦♥•②✳ ✺ ❈♦♥❝•✉s✐♦♥ ❲❡ ❤❛✈❡ ❞❡✈✐s❡❞ ❛ ❝❛•❝✉•✉s ❢♦r ❣❡♥❡r❛t✐♥❣ ♣r✐♠❡ ✐♠♣•✐❝❛t❡s ♦❢ ❝•❛✉s❡ s❡ts ❞❡✜♥❡❞ ♦✈❡r ❡q✉❛t✐♦♥s ❛♥❞ ❞✐s❡q✉❛t✐♦♥s ❜❡t✇❡❡♥ ❝♦♥st❛♥ts✳ ❙✉❝❤ ❛♥ ❛•❣♦r✐t❤♠ ✐s ♠✉❝❤ ♠♦r❡ r❡str✐❝t❡❞ ❛♥❞ t❤✉s ✭♣r♦❜❛❜•②✮ ♠♦r❡ ❡✣❝✐❡♥t t❤❛♥ t❤❡ ♥❛✐✈❡ ❛♣♣r♦❛❝❤ ❝♦♥✲ s✐st✐♥❣ ✐♥ ❛♣♣•②✐♥❣ t❤❡ r❡s♦•✉t✐♦♥ ❝❛•❝✉•✉s ♦♥ ❡q✉❛•✐t② ❛①✐♦♠s✳ ■♥ ♣❛rt✐❝✉•❛r✱ t❤❡ ❝•❛✉s❡s ❛r❡ t❛❦❡♥ ♠♦❞✉•♦ ❡q✉✐✈❛•❡♥❝❡✱ ❛♥❞ ❡✣❝✐❡♥t

  t❤❛♥ ❤❛✈✐♥❣ t♦ ❡①♣r❡ss t❤❡♠ ❡①♣•✐❝✐t•② ❛s ❛①✐♦♠s✳ ❆ ♥❛t✉r❛• ❝❛♥❞✐❞❛t❡ ✐s t❤❡ s✉✲ ♣❡r♣♦s✐t✐♦♥ ❝❛•❝✉•✉s✱ ❜✉t ✐t ✐s ❡❛s② t♦ s❡❡ t❤❛t ✐t ✐s t♦♦ r❡str✐❝t✐✈❡ ❢♦r ♦✉r ♣✉r♣♦s❡✿ ❢♦r ❡①❛♠♣•❡✱ ❛ ♣r✐♠❡ ✐♠♣•✐❝❛t❡ ♦❢ t❤❡ ❝•❛✉s❡s a ≃ b ❛♥❞ c ≃ d✱ ✐s c ≃ a ∨ d ≃ b✱ ②❡t t❤❡r❡ ✐s ♥♦ ✇❛② ♦❢ ❣❡♥❡r❛t✐♥❣ t❤✐s ❝•❛✉s❡ ❢r♦♠ t❤❡ t✇♦ ♣r❡♠✐s❡s ❜② s✉♣❡r✲ ♣♦s✐t✐♦♥✱ ❡✈❡♥ ✐❢ ♦r❞❡r✐♥❣ ❝♦♥❞✐t✐♦♥s ❛r❡ r❡•❛①❡❞✳ ■t ✐s t❤✉s ♥❛t✉r❛• t♦ s❡❛r❝❤ ❢♦r ♥❡✇ s✉♣❡r♣♦s✐t✐♦♥✲❜❛s❡❞ str❛t❡❣✐❡s ❢♦r ❣❡♥❡r❛t✐♥❣ ✐♠♣•✐❝❛t❡s✳ ❖♥ ♦♥❡ ❤❛♥❞✱s✉❝❤ str❛t❡❣✐❡s ♠✉st ❜❡ ♠♦r❡ r❡•❛①❡❞ t❤❛♥ t❤❡ st❛♥❞❛r❞ s✉♣❡r♣♦s✐t✐♦♥ ✐♥❢❡r❡♥❝❡ r✉•❡s✱ s

  ❞❛t❛✲str✉❝t✉r❡s ❛r❡ ✉s❡❞ t♦ r❡♣r❡s❡♥t ❝•❛✉s❡ s❡ts✳ ❆•❣♦r✐t❤♠s ❛r❡ ♣r♦✈✐❞❡❞ ❢♦r ✉♣❞❛t✐♥❣ s✉❝❤ ❞❛t❛✲str✉❝t✉r❡s ❛♥❞ ❞❡t❡❝t✐♥❣ r❡❞✉♥❞❛♥❝②✳ ❚❤❡ ❣❡♥❡r❛t✐♦♥ ♦❢ t❤❡ ✐♠♣•✐❝❛t❡s ✐s ♣❡r❢♦r♠❡❞ ❜② ❛♣♣•②✐♥❣ ❛ ❝♦♥str❛✐♥❡❞ ♣❛r❛♠♦❞✉•❛t✐♦♥ r✉•❡✱ ✇❤❡r❡ ❡q✉❛t✐♦♥s ♣❡r♠✐tt✐♥❣ t❤❡ ❛♣♣•✐❝❛t✐♦♥ ♦❢ t❤❡ tr❛♥s✐t✐✈✐t② ❛①✐♦♠ ❛r❡ ❛••♦✇❡❞ t♦ ❜❡ ❛ss❡rt❡❞ ✐♥st❡❛❞ ♦❢ ❜❡✲ ✐♥❣ ♣r♦✈❡❞✳ ❲❤✐•❡ s✉❝❤ ❛ r✉•❡ ✐s ❢❛r ♠♦r❡ r❡str✐❝t❡❞ t❤❛♥ t❤❡ r❡s♦•✉t✐♦♥ r✉•❡✱ ✐ts ❞❡✜♥✐t✐♦♥ ♣r♦❜❛❜•② •❡❛✈❡s r♦♦♠ ❢♦r ✐♠♣r♦✈❡♠❡♥ts✿ ✐♥ ♣❛rt✐❝✉•❛r✱ ♥♦ ♦r❞❡r✐♥❣ r❡✲ str✐❝t✐♦♥ ✐s ❝♦♥s✐❞❡r❡❞ ✐♥ t❤❡ ❝✉rr❡♥t ✈❡rs✐♦♥✳ ❚❤✐s ✐s t♦ s♦♠❡ ❡①t❡♥t ✉♥❛✈♦✐❞❛❜•❡ ✐❢ ♦♥❡ ✇❛♥ts t♦ ❞❡r✐✈❡ ❛•• ✐♠♣•✐❝❛t❡s ✸ ✱ ❜✉t✱ st✐••✱ s♦♠❡ ♣❛rt✐❛• ♦r❞❡r✐♥❣ ❝♦♥❞✐t✐♦♥s ❝❛♥ ♣r♦❜❛❜•② ❜❡ ❡♥❢♦r❝❡❞ ✇❤✐•❡ r❡t❛✐♥✐♥❣ ❝♦♠♣•❡t❡♥❡ss✱ ❛♥❞ ✇❡ ❛r❡ ❝✉rr❡♥t•② ✸ ❋♦r ✐♥st❛♥❝❡✱ ❣✐✈❡♥ t❤❡ ❝•❛✉s❡s a ≃ b ∨ c ≃ d ❛♥❞ c ≃ d✱ ♦♥❡ s❤♦✉•❞ ❜❡ ❛❜•❡ t♦ ♣r♦✈❡ t❤❛t a ≃ b ✐s ❛♥ ✐♠♣•✐❝❛t❡✱ ❤❡♥❝❡ t♦ ♣❡r❢♦r♠ ✐♥❢❡r❡♥❝❡s ♦♥ c ≃ d ❡✈❡♥ ✐❢ a, b ≻ c, d✳