Qinghua Zhang
email: qinghua.zhang@inria.fr

Lennart Ljung
email: ljung@isy.liu.se

From Structurally Independent Local LTI Models to LPV Model

Keywords: System identification, LPV model, coherent local linear models

The local approach to linear parameter varying (LPV) system identification consists in interpolating individually estimated local linear time invariant (LTI) models corresponding to fixed values of the scheduling variable. It is shown in this paper that, without any global structural assumption of the considered LPV system, individually estimated local state-space LTI models do not contain su cient information for determining similarity transformations making them coherent. It is possible to estimate these similarity transformations from input-output data under appropriate excitation conditions.

Introduction

Linear parameter varying (LPV) models provide an effective approach to handling nonlinear control systems [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF][START_REF] Mohammadpour | Control of Linear Parameter Varying Systems with Applications[END_REF][START_REF] Lopes Dos Santos | Linear Parameter-Varying System Identification[END_REF][START_REF] Sename | Robust Control and Linear Parameter Varying Approaches[END_REF]. Some successful methods for LPV system identification have been reported recently [START_REF] Van Wingerden | Subspace identification of bilinear and LPV systems for open-and closed-loop data[END_REF][START_REF] Mercere | Continuous-time linear parameter-varying identification of a cross: A local approach[END_REF][START_REF] Lopes Dos Santos | An LPV modeling and identification approach to leakage detection in high pressure natural gas transportation networks[END_REF][START_REF] Tóth | Instrumental variable scheme for closed-loop LPV model identification[END_REF][START_REF] Zhao | Prediction error method for identification of LPV models[END_REF][START_REF] Piga | Identification of LPV models from noise-corrupted scheduling parameter and output observations[END_REF]. In the local approach to LPV system identification, interpolation is essential to establishing global models from a collection of locally estimated linear time invariant (LTI) models [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF][START_REF] De Caigny | Interpolation-based modeling of MIMO LPV systems[END_REF][START_REF] De Caigny | Interpolated modeling of LPV systems[END_REF]. As LTI state-space models can be estimated in an arbitrary state basis, it is necessary to use a coherent collection of local models for the purpose of interpolation. This paper is focused on the problem of making local state-space models coherent, without treating the interpolation step. Only state-space models are considered in this paper, as local model coherence is not relevant for other models. For shorter expressions, the words "statespace" will be omitted from terms like "local state-space model" and "LTI state-space model". In practice, interpolation is based on a finite set of local LTI models, each corresponding to a specific value of the scheduling variable p(t), the main discussion of this paper is thus about the case where p(t) evolves within a finite set, but its motivation is indeed with the perspective of interpolation for continuous values of p(t). It seems natural to transform all the local LTI models to some canonical form in order to make them coherent.

The main purpose of this paper is to point out the fact that, in the local approach to LPV system identification, structurally independent local LTI models themselves do not contain su cient information to determine similarity transformations making them coherent. However, locally estimated LTI models can be made coherent by making use of the information contained in some inputoutput data sequences across all the working points, notably with an algorithm initially introduced in the framework of piecewise linear hybrid systems [START_REF] Verdult | Subspace identification of piecewise linear systems[END_REF][START_REF] Zhang | LPV System Common State Basis Estimation from Independent Local LTI Models[END_REF]. Preliminary results of this work have been presented in [START_REF] Zhang | LPV System Common State Basis Estimation from Independent Local LTI Models[END_REF], which are completed in the present paper with a rigorous proof of the main result.

2 Problem statement Let u(t) 2 R q and y(t) 2 R s be respectively the input and the output at discrete time instant t = 0, 1, 2, . . . , p(t) be the scheduling variable evolving within a compact set ⇧. An LPV system is described by the state-space model

x(t + 1) = A(p(t))x(t) + B(p(t))u(t) + w(t) (1a)
y(t) = C(p(t))x(t) + D(p(t))u(t) + v(t) (1b)
where

x(t)
Q(p(t)), R(p(t)).
Based on the fact that the LPV system (1) becomes an LTI system when the scheduling variable p(t) is maintained at a fixed value, the following definitions aim at establishing a link between LPV and LTI models.

Consider a set of m LTI models indexed by the integer i:

x(t + 1) = A i x(t) + B i u(t) + w(t) (2a) y(t) = C i x(t) + D i u(t) + v(t) (2b)
characterized by matrices A i , B i , C i , D i of appropriate sizes, and noise covariance matrices

Q i and R i . The notation i , (A i , B i , C i , D i , Q i , R i ) (3)
will be used to denote the matrices characterizing the ith local LTI model [START_REF] De Caigny | Interpolation-based modeling of MIMO LPV systems[END_REF], or the LTI model itself by abuse of notation. The set of LTI models will be denoted by

⌃ = { i : i = 1, 2, . . . , m}. (4)
Definition 1 A set of local LTI models

⌃ ⇤ = {(A ⇤ i , B ⇤ i , C ⇤ i , D ⇤ i , Q ⇤ i , R ⇤ i ) : i = 1, 2, . . . , m} (5)
is called a multi-snapshot of the LPV system (1) for

p(t) 2 P = {p 1 , . . . , p m } ⇢ ⇧, (6)
if

A ⇤ i = A(p i ), B ⇤ i = B(p i ), C ⇤ i = C(p i ), D ⇤ i = D(p i ), Q ⇤ i = Q(p i ), R ⇤ i = R(p i ). (7)
2 In the local approach to LPV system identification [14, 2, 1], a set of locally estimated LTI models are interpolated to obtain a global model. As such local LTI models are typically estimated up to di↵erent and arbitrary similarity transformations, they do not constitute a multisnapshot of the underlying LPV system in the sense of Definition 1. It is thus important to make the local models "coherent" before the interpolation step. What does mean a "coherent" set of local LTI models with the perspective of their interpolation? The following auxiliary definition will be helpful.

Definition 2

The input-output behavior of a set of local LTI models ⌃ = { 1 , . . . , m } is the input-output behavior of the multi-model switching system consisting of the same set of LTI models, such that the i-th LTI model i is active when p(t) = p i 2 P = {p 1 , . . . , p m }, and at every transition between two LTI models, the initial state of the new active model is equal to the final state of the previous active LTI model.

2 For example, following this definition, the input-output behavior of a multi-snapshot ⌃ ⇤ (see Definition 1) of an LPV system (1) is identical to input-output behavior the LPV system when p(t) evolves within the restricted set P.

If the interpolation of a set of local LTI models is expected to describe correctly an LPV system for all sequences of p(t) within ⇧, it should also be true in the particular case where p(t) evolves within the restricted set P = {p 1 , . . . , p m } ⇢ ⇧, including when p(t) switches between di↵erent values within P. It means that a "coherent" set of local LTI models should have the same inputoutput behavior as the underlying LPV system when p(t) 2 P, in the sense of Definition 2. This requirement will be satisfied by Definition 4 through Property 1.

Definition 3 Two LTI models , (A, B, C, D, Q, R) and ˜ , ( Ã, B, C, D, Q, R) are related by a similarity transformation characterized by an invertible matrix T 2 R n⇥n and denoted as

T ! ˜ , (8)
if the matrices characterizing the two LTI models satisfy

à = T AT 1 , B = T B, C = CT 1 , (9a) D = D, Q = T QT T , R = R. (9b) 2 Definition 4 A set of LTI models ⌃ = {˜ i : i = 1, 2, . . . , m}
is said coherent with another set

⌃ = { i : i = 1, 2, . . . , m},
if there exits an invertible transformation matrix T 2 R n⇥n , common to the local models, such that i T ! ˜ i for all i = 1, . . . , m. [START_REF] Piga | Identification of LPV models from noise-corrupted scheduling parameter and output observations[END_REF] This relationship between ⌃ and ⌃ is then denoted by

⌃ T =) ⌃. ( 11
)
The local LTI models ˜ i are simply said coherent when the reference model set ⌃ is obvious. 2 The relevance of this definition is justified by the following property, as discussed before Definition 3.

Property 1 If a set of local models ⌃ is coherent with the multi-snapshot ⌃ ⇤ of the LPV system (1), then ⌃ has the same input-output behavior (in the sense of Definition 2) as that of the LPV system (1) with the scheduling variable p(t) restricted to the finite set {p 1 , . . . , p m } and with appropriate initial states.

2 The proof of this property is trivial: under the assumed conditions, the LPV system (1) and its multi-snapshot ⌃ ⇤ have the same input-output behavior, and the states of ⌃ ⇤ and ⌃ are related by the transformation matrix T . Property 1 is a necessary condition that a relevant definition of local model coherence should satisfy. It does not exclude other possible definitions, notably those based on p-dependent state transformation matrices T (p). Such considerations, related to LPV system equivalent state-space representations as investigated in [START_REF] Kulcsar | On the similarity state transformation for linear parameter-varying systems[END_REF], would be out of the scope of this technical communiqué. In practice, local LTI models are estimated from a finite data sample subject to random uncertainties, thus the definition of coherent local models is understood in an approximative sense. If the estimation of each local model is consistent, then Definition 4 can also be understood for the limiting models when the data size for each local model estimation tends to infinity. If some global structural assumptions of the matrix functions A(p), B(p), etc. were assumed, then they could be used to make estimated local LTI models coherent. This paper is focused on structurally independent local LTI models, as defined below. Definition 5 A set of local LTI models i are structurally dependent if their parametrizations are such that fixing the matrices A i , B i , C i , D i , Q i , R i reduces the degrees of freedom of the matrices A j , B j , C j , D j , Q j , R j , for j 6 = i; otherwise they are structurally independent. 2 For instance, if a set of local LTI models is parametrized such that all the matrices A i share an equal entry at the same position, say A i (1, 1) = A j (1, 1) for all i, j = 1, 2, . . . , m, then the local models are structurally dependent. A less trivial counterexample will be given in Section 3.2 with equation (27).

Making local LTI models coherent

Apply di↵erent similarity transformations Ti to a set of local LTI models ⌃ ⇤ , resulting in a transformed set of models ⌃, which is not coherent with ⌃ ⇤ . Assume that neither Ti nor ⌃ ⇤ is known, is it possible to derive, solely from ⌃, a model set ⌃ coherent with ⌃ ⇤ ? 3.1 Structurally independent local LTI models Proposition 1 Given a set of structurally independent LTI models

⌃ ⇤ = { ⇤ i : i = 1, 2, . . . , m}, (12)
assume that a set of arbitrary unknown invertible matrices

{ Ti 2 R n⇥n : i = 1, 2, . . . , m} (13)
transform the LTI models

⇤ i Ti ! ˆ i , (14)
resulting in a set of transformed models

⌃ = {ˆ i : i = 1, 2, . . . , m}, (15)
then it is impossible to determine invertible matrices Ti 2 R n⇥n (not necessarily equal to Ti ) for i = 1, . . . , m, solely based on the set of transformed models ⌃, such that the further transformed set of models

⌃ = {˜ i : i = 1, 2, . . . , m} (16)
through

ˆ i Ti ! ˜ i ( 17
)
is coherent with the original set ⌃ ⇤ , or more precisely,

⌃ ⇤ T =) ⌃ ( 18
)
for some unknown T 2 R n⇥n . 2

Proof. The proof proceeds by contradiction, starting by assuming the contrary of Proposition 1, in which the word "impossible" is replaced by "possible", or more compactly, "solely based on ⌃, the original model set ⌃ ⇤ can be determined up to an unknown transformation matrix T common to all local LTI models". As the matrix T 2 R n⇥n can be any invertible matrix, what is determined from ⌃, is a set of (infinitely many) sets of LTI models, which will be denoted by S( ⌃) in what follows. The contrary of Proposition 1 then means

S( ⌃) = {⌃ : ⌃ ⇤ T =)⌃, T 2 R n⇥n , det(T ) 6 = 0}. ( 19
)
Within this set, the particular element ⌃ ⇤ corresponds to the identity transformation matrix T , hence

⌃ ⇤ 2 S( ⌃). ( 20
)
Consider a set of m invertible and pairwise distinct matrices S i 2 R n⇥n , i = 1, 2, . . . , m and a model set

⌃ 0 = { 0 i : i = 1, 2, . . . , m} (21)
such that

0 i Si ! ⇤ i ( 22
)
for i = 1, 2, . . . , m. The two transformations ( 22) and ( 14) then imply

0 i T 0 i ! ˆ i ( 23
)
with T 0 i = Ti S i , for i = 1, 2, . . . , m. This result means that the model set ⌃ originating from ⌃ ⇤ through transformations ( 14) could also have originated from ⌃ 0 through transformations (23). Therefore, following the same reasoning as at the beginning of this proof, that has led to ⌃ ⇤ 2 S( ⌃), it yields ⌃ 0 2 S( ⌃). This result and [START_REF] Hoyt | Gain scheduling for geometrically nonlinear flexible space structures[END_REF] imply that ⌃ ⇤ T =)⌃ 0 for some invertible T . In terms of the LTI models composing ⌃ ⇤ and ⌃ 0 , it means

⇤ i T ! 0 i (24)
with a common T for all i = 1, 2, . . . , m. This result is clearly in contradiction with (22), where S i are pairwise distinct matrices. This contradiction invalidates the assumption made at the beginning of this proof, hence proving Proposition 1. 2

Examples based on global structural assumptions

The result of Proposition 1 may seem in contradiction with some known methods for making local LTI models coherent. In fact, these methods assume (implicitly) some particular structure of the matrix-valued functions A(p), B(p) etc., therefore, the local LTI models are not structurally independent in the sense of Definition 5. To better clarify the situations, some examples are recalled below, by pointing out their particular global structural assumptions.

Coherent LTI models based on canonical forms

In order to make local LTI models coherent, a natural idea is to find similarity transformations of the local LTI models leading to some LTI canonical state-space form, or to estimate local LTI models directly in such a form [START_REF] Hoyt | Gain scheduling for geometrically nonlinear flexible space structures[END_REF][START_REF] Steinbuch | Experimental modelling and LPV control of a motion system[END_REF][START_REF] Lovera | Identification for gain-scheduling: a balanced subspace approach[END_REF][START_REF] De Caigny | Interpolation-based modeling of MIMO LPV systems[END_REF]. This practice assumes that the local LTI models are coherent when they are all transformed into the same LTI canonical form 1 .

For the sake of presentation simplicity, consider the case of single-input-single-output (SISO) models. In the controllable form, the m local models involve, for .

i = 1, 2, . . . , m, Ãi = 2 6 6 6 6 6 4
0 1 0 • • • 0 0 0 1 • • • 0 . . . . . . . . . • • • . . . a (1) i a (2) i a (3) i • • • a (n) i 3
The assumption that the local LTI models in this canonical form are coherent implies that there exists a common invertible matrix T 2 R n⇥n such that, for all i = 1, 2, . . . , m, Ãi = T A ⇤ i T 1 . Let u(M ) denote the sub

1
Transforming the local LTI models to the same LTI canonical form does not lead to a canonical form of the underlying LPV system, as pointed out in [START_REF] Tóth | Discrete time LPV I/O and state space representations, differences of behavior and pitfalls of interpolation[END_REF]. matrix of M excluding its last row, then

u T A ⇤ i T 1 = 2 6 6 6 4 0 1 0 • • • 0 . . . . . . . . . • • • . . . 0 0 0 • • • 1 3 7 7 7 5 , 8i = 1, 2, . . . , m. (26)
By assuming that u T A ⇤ i T 1 are equal to the same matrix for all i = 1, 2, . . . , m, the local LTI models are not structurally independent in the sense of Definition 5.

Coherent LTI models based on the observability matrix

In [START_REF] De Caigny | Interpolated modeling of LPV systems[END_REF] another method is proposed to make local LTI models coherent, in a framework more general than the one of the present paper. In the particular case of SISO observable systems, this method consists in finding different similarity transformations so that the m transformed local LTI models all have the same observability matrix. It implies the assumption that the local LTI systems composing a multi-snapshot of the underlying LPV system all have the same observability matrix, or more explicitly,

C i A s i = C j A s j (27)
for all i, j 2 {1, 2, . . . , m} and s 2 {0, 1, . . . , n 1}. Again the local LTI models are not structurally independent in the sense of Definition 5.

Remark. This observability matrix-based method is incompatible with the previously presented example of canonical form-based method, in the sense that in general the two resulting model sets are not coherent with each other. Moreover, the canonical observable formbased method is incompatible with the canonical controllable form-based method, and similarly the observability matrix-based method is incompatible with the controllability matrix-based method. These incompatibilities between "natural" methods illustrate the fact that there is no generally natural global structural assumption for making estimated local LTI models coherent. 2

4 Data-based transformations for coherent LTI models In order to make local LTI models coherent, as an alternative to global structural assumptions, the information contained in some data sequences can be used. The idea is based on Property 1 that is satisfied by coherent local LTI models. In particular, the state trajectory "continuity" in Definition 2 is essential: when the scheduling variable p(t) switches from one value to another within P, the local LTI model switches accordingly, but the initial state of the new active LTI model is equal to the final state of the previous active LTI model. Given a data set D = {(p(t), u(t), y(t)) : t = 1, 2, . . . , N} with p(t) 2 P, start by segmenting the data sequence into pieces of constant p(t) value. If the local LTI models are estimated from pieces of this data set, the whole state trajectory on each data segment can be obtained as a by-product of the identification procedure [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation Applications[END_REF]. If the local LTI models have been estimated from another data set, then the initial state and the whole state trajectory on each data segment of D can be estimated with the already estimated local LTI models [START_REF] Ljung | System Identification -Theory for the User[END_REF]. Similarity transformations of the local LTI models are then determined so that, when the local LTI model switches, the initial state of the new active LTI model is equal to the final state of the previous active LTI model. A complete data-based algorithm for making local LTI models coherent has been introduced in [START_REF] Verdult | Subspace identification of piecewise linear systems[END_REF] in the framework of piecewise linear hybrid systems. A variant algorithm has been presented in [START_REF] Zhang | LPV System Common State Basis Estimation from Independent Local LTI Models[END_REF].

Conclusion

It has been shown in this paper that structurally independent local LTI models estimated from data collected with fixed scheduling variable of an LPV system do not contain su cient information to make themselves coherent with the perspective of LPV model interpolation. Nevertheless, well designed data sequences can be used to address this problem.

Preprint submitted to Automatica


Masoud Pourasghar
email: mpourasgharlaf@iri.upc.edu

Vicenç Puig
email: vpuig@iri.upc.edu

Carlos Ocampo-Martinez

Qinghua Zhang
email: qinghua.zhang@inria.fr

Reduced-order Interval-observer Design for Dynamic Systems with Time-invariant Uncertainty

Keywords: Uncertain dynamic system, interval observer, set-based observer, reduced-order observer

This paper addresses the design of reduced-order interval-observers for dynamic systems with time-invariant uncertainty. Because of the limitations of using the set-based approach to preserve the time dependency of parameter uncertainty and the wrapping effect to deal with interval-observers, the trajectory-based interval-observer approach is used with an appropriate observer gain. But, there could be some difficulties to satisfy the conditions for selecting a suitable gain to guarantee the positivity of the resulting observer. Then, a reducedorder observer is designed to reduce the computational complexity and to increase the degree of freedom when selecting the observer gain. Finally, a simulation example is employed for illustrating and analyzing the effectiveness of the proposed approach.

INTRODUCTION

Interval-observer approaches make use of the interval hull of the approximated set instead of the exact set. Moreover, in the set-based approach, the propagation of the state set is affected by several problems such as the wrapping effect, temporal variance on uncertain parameters (or uncertain parameter time dependency) and range evaluation of an interval function, especially in the case of using the interval hull of the set at each iteration. Therefore, conservative and unstable results may obtain (for even a stable system) with using the set-based approach in the simulation of the system with parametric time-invariant uncertainties [START_REF] Puig | Observers for interval systems using set and trajectory-based approaches[END_REF]]. On the other hand, the approximated state set can be computed based on a set of point-wise trajectories. This type of approach is called trajectory-based approach [START_REF] Puig | Observers for interval systems using set and trajectory-based approaches[END_REF]]. The advantage of trajectory-based approach in comparison with the set-based approach is to overcome the wrapping effect due to generating the real trajectories based on selecting the particular value of the uncertain parameters. Meanwhile, the uncertain parameter time dependency can be preserved if the set of point-wise trajectories are generated [START_REF] Puig | Observers for interval systems using set and trajectory-based approaches[END_REF]].

Based on the literature, interval dynamic systems are those uncertain systems whose uncertain parameters are bounded by intervals [START_REF] Kolev | Assessing the stability of linear time-invariant continuous interval dynamic systems[END_REF] and [START_REF] Le | Zonotopic set-membership estimation for interval dynamic systems[END_REF]]. Particularly, the wrapping effect appears in that kind of interval systems that are not monotonic [START_REF] Puig | Observers for interval systems using set and trajectory-based approaches[END_REF]]. It means that if there are some negative elements in the state matrix related to the statespace model of the observer (i.e., it is not monotonic), the wrapping effect affects the interval-observer. As it is mentioned before, one of the possible solutions to overcome the wrapping effect is to make use of trajectory-based observers, which will be further reviewed later in this paper. But still, there are some difficulties related to the computational complexity and designing the suitable observer gain in the trajectory-based observer.

In order to reduce the computational complexity of the trajectory-based approach, the order of the observer can be reduced and unknown states of the system can be separately estimated by using a reduced-order observer, since the problem of wrapping effect appears only when those elements of the state matrix that are related to the unmeasurable states are negative values. Besides, by using the reduced-order observer not only the problem of wrapping effect can be solved with less computational complexity, but also, designing the observer gain can be done with more degree of freedom in comparison with the full-order observer.

x(k + 1) = A(θ)x(k) + B(θ)u(k), (1a) y(k) = C(θ)x(k), (1b)
where k ∈ Z ≥0 indicates the discrete time, x ∈ R nx is the state vector, u ∈ R nu and y ∈ R ny denote the input and the output vectors, respectively. The system matrices of appropriate dimensions are A(θ) ∈ R nx×nx , B(θ) ∈ R nx×nu and C(θ) ∈ R ny.×nx . Moreover, θ is the vector of time-invariant uncertain parameters with the bounded values by an compact set Θ as

Θ = θ ∈ R n θ | θ i ≤ θ i ≤ θ i , i = 1, ... , n .
(2) Furthermore, the system matrices A(θ), B(θ) and C(θ) are not exactly known and the uncertainty is included into them based on (2) as

A(θ) ≤ A(θ) ≤ A(θ), (3a)
B(θ) ≤ B(θ) ≤ B(θ), (3b)
C(θ) ≤ C(θ) ≤ C(θ), (3c)
where • and • denote the lower and upper bound of each matrix, respectively. Note that the inequalities in (3) should be understood as element-wise inequalities.

Therefore, the uncertainties in matrices A(θ), B(θ) and C(θ) can be decomposed as A(θ) = A e + ∆A(θ), B(θ) = B e + ∆B(θ) and C(θ) = C e + ∆C(θ), where A e , B e and C e are the nominal part of the system matrices A(θ), B(θ) and C(θ), respectively. Furthermore, ∆A(θ), ∆B(θ) and ∆C(θ) denote the effect of uncertainties in A(θ), B(θ) and C(θ), respectively. Remark 1. The time invariance of uncertain parameters is usually not taken into account in the literature, being this the main goal of this paper.

Assumption 1. B(θ) and C(θ) are assumed monotonic with respect to θ. It means, there exist two functions F(y(k), u(k), θ) and F(y(k), u(k), θ) where the monotonicity property is used to build them.

If (1) is observable, the Luenberger observer is written as

x(k + 1) = A(θ)x(k) + B(θ)u(k) + L(y(k) -ŷ(k)), (4a) ŷ(k) = C(θ)x(k), (4b)
where ŷ ∈ R ny and x ∈ R nx are the estimated system output and state, respectively. Furthermore, the state observer (4) can be written as x(k + 1) = (A(θ) -

LC(θ))x(k) + [B(θ) L] u(k) y(k) . By denoting A obs = A(θ) - LC(θ), B obs = [B(θ) L] and u obs = u(k) y(k)
, the state observer can be rewritten as x(k + 1) = A obs (θ)x(k) + B obs (θ)u obs (k).

(5) Additionally, the observer gain L for θ ∈ Θ should be defined such that A obs was a Schur matrix. Based on the observer expression in (5), the effect of the uncertainty is introduced in the observer through A obs and B obs matrices. There are two main approaches to take into account the uncertainty in the output/state estimation as it will be recalled in the following sections.

Set-based vs Trajectory-based Approaches

In set-based approaches, the set of states at time instant k + 1 is approximated by using propagation algorithms from the set of states at time k [START_REF] Puig | Observers for interval systems using set and trajectory-based approaches[END_REF], [START_REF] Combastel | Zonotopes and Kalman observers: Gain optimality under distinct uncertainty paradigms and robust convergence[END_REF]]. However, in set-based approaches, the parameter time invariance is not taken into account during the propagation of the uncertainties. Uncertain parameters are unknown but bounded in their uncertainty intervals and can vary arbitrarily at each time instant within the bounded interval in the set-based approach. Moreover, some additional problems appear in the case of using setbased approaches such as wrapping effect and range evaluation of an interval function [START_REF] Puig | Worst-case simulation of discrete linear time-invariant interval dynamic systems[END_REF]] as already discussed. One way to solve the time-dependency problem of the set-based approach is to derive the relation between states and parameters, which can bring the system from the initial state to the current state. Considering this idea, the wrapping effect problem is solved while the time dependency of the uncertain parameter can be preserved. In the trajectory-based approach, the value of parameter uncertainty is unknown but bounded in the interval and the invariance of parameter uncertainty can be guaranteed at each time instant. The interval of the states can be estimated at each iteration by using the particular state trajectories corresponding to particular values of uncertainty θ. Based on [START_REF] Puig | Worst-case simulation of discrete linear time-invariant interval dynamic systems[END_REF], the relation between states and uncertain parameters is derived as

x(k) = (A obs (θ)) k x(0)+ k-1 j=0 (A obs (θ)) k-1-j B obs (θ) u obs (j).
(6) In ( 6), the uncertainty always is propagated from the initial state to avoid the wrapping effect. Furthermore, by assuming that A obs (θ) at time instant k + 1 is equal to the A obs (θ) at time instant k and B obs (θ) at time instant k + 1 So far, however, several studies have reported state esti-mation based on trajectory-based interval observer, but there is still insufficient study for reducing the compu-tational complexity and the design of the suitable ob-server gain in order to satisfy the positivity conditions. Therefore, the main contribution of this paper is to de-sign the reduced-order interval observer to estimate the unmeasurable states. Furthermore, the observer gain can be designed with more degree of freedom to obtain the monotonic observer in spite of non-monotonicity of the system state matrix. Simultaneously, the time invariance of the uncertain parameter can be preserved.

Regarding the structure of the paper, Section 2 deals with preliminaries and the problem formulation of interval state observation. Moreover, the set-based and trajectory-based observers are introduced, respectively. The design of the interval observer is proposed in Section 3. Set-based, trajectory-based and reduced-order observers are implemented and compared by using an academic example for both monotonic and non-monotonic systems and the effect of the observer gain will be discussed in Section 4.

Finally, the main conclusions are drawn in Section 5.

PROBLEM FORMULATION OF INTERVAL STATE OBSERVATION

Problem Statement

The uncertain discrete linear time-invariant model is rep-resented by the following state-space form:

x(k) = max θ∈Θ ( Âobs + ∆A obs (θ)) k x(0) + k-1 j=0 ( Âobs + ∆A obs (θ)) k-1-j ( Bobs + ∆B obs (θ)) u obs (j) , (7a) x(k) = min θ∈Θ ( Âobs + ∆A obs (θ)) k x(0) + k-1 j=0 ( Âobs + ∆A obs (θ)) k-1-j ( Bobs + ∆B obs (θ)) u obs (j) , (7b)
where Âobs , Bobs and ∆A obs (θ), ∆B obs (θ) are nominal and uncertainties of the observer matrices, respectively.

DESIGNING THE INTERVAL OBSERVER

One possible way to avoid solving the optimization problems in ( 7) is to design the observer gain L in such a way to force A obs into an element-wise positive matrix despite of the negative elements in A(θ). Thus, designing the observer gain in this manner imposes some additional constraints. In case the observer is not possible to be designed, a reduced-order observer can be used to estimate only those states that are affected by the negative elements of A(θ). Moreover, the observer gain L can be designed with more degree of freedom for the reduced-order observer in comparison with full-order observer. In other words, the number of conditions and constrains of selecting the observer gain can be reduced by using the reduced-order observer.

Designing Observer Gain

In the case that one element of A(θ) is negative, A obs will satisfy the monotonicity property if

a ij -(LC) ij ≥ 0, (8)
where the subindices ij denote that there exist one state such xi whose variation with respect to another state xj is negative (e.g., a ij < 0). Furthermore, (LC) ij shows the i-th row and j-th column of matrix LC ∈ R nx×nx that is obtained by multiplication of matrix L ∈ R nx×ny into the matrix C(θ) ∈ R ny×nx . Therefore, the condition in (8) will be satisfied if

a ij ≥ (LC) ij . (9)
Consequently, by forcing the observer gain to eliminate the negative elements of A obs , a monotonic interval observer can be obtained. Furthermore, the right side of (9) can be written as (LC) ij = L i C j , where L i shows the i-th row elements of the matrix L and C j indicates the j-th column elements of the matrix C(θ). Obviously, the elements C j can not be chosen freely, because it is determined by the observer model. Therefore, in order to force the condition in (9), only the elements L i can be designed freely.

Moreover, the interval-observer convergence should be considered when designing the observer gain. Therefore, the observer gain matrix L can be divided into two matrices L -and L + , where

• L -shows the elements that are chosen to force the condition (9), • L + shows the elements that are chosen to guarantee the observer convergence.

Thus, the observer gain can be split into

L = L + + L -. (10)
Similarly, the matrix A(θ) can be written as

A(θ) = A + (θ) + A -(θ), (11
) where, A + (θ) and A -(θ) are determined by the positive and negative elements of matrix A(θ), respectively. Moreover, by considering ( 10) and ( 11), A obs can be written as

A obs (θ) = A + (θ) -L + C(θ) + A -(θ) -L -C(θ) . (12)
Furthermore, by denoting the index mn for those elements of matrices A(θ) and LC(θ) that do not have any effect on the monotonicity property of the A obs , each part of (12) will be determined as

if a mn > 0 =⇒    A + (θ) mn = a mn , A -(θ) mn = 0, L + C(θ) mn = L mα C αn , (13)
where L mα are the m-th row elements of matrix L and C αn are the n-th column of matrix C(θ) and

if a ij < 0 =⇒      A + (θ) ij = 0, A -(θ) ij = a ij , L + C(θ) ij = 0. (14)
Furthermore, for the positive elements, A -(θ) mn = L -C(θ) mn . Therefore, the elements of matrix L -C are either positive or null by assuming the condition in (9).

Taking into account that the condition (9) should be satisfied to achieve the positivity of A obs , it can be written as A obs = A -L + C(θ), where A = A(θ) -L -C. Therefore, L -should satisfy the positivity conditions and L + should be designed in order to guarantee the convergence of the observer. Thus, two state observers can be designed to estimate the maximum and minimum state trajectories. Therefore, A obs , for each observer, is

A obs = A -L + C, (15a) A obs = A -L + C, (15b
) where • and • denote that the matrices are related to the upper and lower trajectory observers, respectively.

In this paper, the upper L + and lower L + observer gains are designed using Linear Matrix Inequalities (LMIs). Therefore, L + and L + should be designed in such way that the eigenvalues of the A -L + C and A -L + C are placed in a stable LMI region. Theorem 1 can be used to designed such an observer gains. Theorem 1. Consider the observer (5) satisfying the condition in (9) and the following LMI inequalities:

is equal to B obs (θ) at time instant k, the time invariance is guaranteed all iterations.

The state estimation x ˆ(k) in ( 6) can be bounded by the

interval hull 2X ˆ(k) = [x ˆ(k), x ˆ(k)],
that can be obtained by solving the following optimization problems:

-rK

qK + K T A -W T C qK + A T K -C T W -rK < 0, (16a) -rK qK + K T A -W T C qK + A T K -C T W -rK < 0, (16b)
where q and r denote the center and the radius of a LMI region, respectively, K is the unknown symmetric matrix, W = L + K is related to the upper observer and W = L + K is related to the lower observer. Obtaining L + and L + as

L + = (W K -1 ) T , (17a)
L + = (W K -1 ) T , (17b) the A obs is Schur and element-wise positive matrix and hence the convergence of the upper and lower observers is guaranteed.

Proof. The proof follows from the application of results presented in [START_REF] Chilali | Robust pole placement in LMI regions[END_REF] to the interval observer case.

Reduced-order Observer

All the states can be estimated together by using a fullorder observer. In this regard, the observer gain can be designed to eliminate the effect of negative elements of A on A obs . However, there exist some cases that designing the observer gain to satisfy the condition ( 9) is not possible or the LMI inequalities in ( 16) are to be solved.

Consequently, the trajectory-based observer will produce wrong state estimation when the unmeasurable states are affected by the negative elements of A. As an alternative, a reduced-order observer is used to reduce the number of conditions and constrains in the LMIs in ( 16). Hence, it is not necessary to design the observer to estimate all the states that are directly measured by sensors. If this is the case, Theorem 2 can be used to design a reduced-order observer for estimating the unmeasured states. Theorem 2. By separating the original states x into measurable x a and unmeasurable x b states such that x(k) = [xa(k) . . . x b (k)] T , (1) can be written as

x a (k + 1) • • • x b (k + 1) = A aa (θ) A ab (θ) A ba (θ) A bb (θ) x a (k) • • • x b (k) + B a (θ) • • • B b (θ) u(k), (18a)
y(k) = I . . . 0 x a (k) • • • x b (k) , ( 18b
)
where x a ∈ R nx a and x b ∈ R nx b . Moreover, the dimensions of the matrices are

A aa ∈ R nx a ×nx a , A ab ∈ R nx a ×nx b , A ba ∈ R nx b ×nx a , A bb ∈ R nx b ×nx b , B a ∈ R nx a ×nu , B b ∈ R nx b ×nu .
Considering the vector of uncertain parameters θ whose values are bounded by Θ (i.e., θ ∈ Θ) or by an interval [θ, θ] (i.e., θ ∈ [θ, θ]), the following reduced-order interval observer can be used for estimating x b :

xb (k + 1) = ζ(k + 1) + Lx a (k + 1), (19)
where Therefore, x b can be estimated by ( 19).

ζ(k + 1) = A bb (θ) -LA ba (θ) ζ(k) + A bb (θ) -LA ab (θ) + A ba (θ) -LAaa(θ) y(k) + B b (θ) -LBa(θ) u(k).
Proof. The proof follows from the extension of standard reduced observer results available as e.g. in [START_REF] Ostertag | Mono-and Multivariable Control and Estimation Linear, Quadratic and LMI Methods[END_REF] to the reduced interval observer case.

Hence, the obtained interval observer in ( 19) is called reduced-order interval observer, which will be able to estimate the unmeasurable states. In other words, if A contains the negative elements, this method can be used to separate the measurable and unmeasurable states. Thus, the only concern will be the unmeasurable states. Therefore, designing the observer gain in order to compensate the effect of negative elements on the unmeasurable states can be done for this reduced-order observer with more degree of freedom as in Section 3.1. Remark 2. If there are still some cases where selecting the observer gain in order to eliminate the effect of negative elements is not possible, a change of coordinates as proposed by [START_REF] Raïssi | Interval state estimation for a class of nonlinear systems[END_REF] may be helpful to achieve the positivity conditions and the successful interval-observer design.

ILLUSTRATIVE EXAMPLE

In order to illustrate the approach proposed in previous sections, an academic example is considered based on the dynamical model (1) with

A(θ) = 0 + θ 11 0.1 + θ 12 0.3 + θ 13 0 + θ 21 0.8 + θ 22 0.2 + θ 23 0.01 + θ 31 0 + θ 32 0.8 + θ 33 , B(θ) = 0 + θ 11 0 + θ 21 1 + θ 31 , C(θ) = 0 + θ 11 0 + θ 21 1 + θ 31 ,
where the time-invariant uncertain parameters are bounded by the following intervals: • θ 11 , θ 21 , θ 22 , θ 31 , θ 32 ∈ [-0.0015, 0.0015],

• θ 12 , θ 13 , θ 23 , θ 33 ∈ [-0.0150, 0.0150].

Two different scenarios are considered in this section:

(1) In the case of monotonic system that means all the elements of A are positive, (2) In the case of a non-monotonic system that at least one of the elements in A is negative.

Therefore, the set-based, trajectory-based and reducedorder observers are designed to estimate the states for each scenario. In addition, the observer gain L is designed by the well-known LMI pole placement method for each case. Furthermore, the proposed observer-gain design method in Section 3.1 is used in the second scenario. In Fig. 1, the simulation results that are obtained by using a setbased approach based on zonotopes in both positive and negative scenarios are presented. In this type of intervalobserver approach, the uncertainty is modeled by using the zonotopic set. Moreover, Fig. 1a shows the obtained estimation of the state x 1 in the case that all the elements of A are positive. On the other hand, Fig. 1b shows the result in the case that at least one of the elements of A that is related to the unmeasurable states (in this case x 1 and x 2 ) is negative. The wrapping effect can be observed in the non-monotonic case since the state estimations of upper and lower bounds are diverging in the second scenario.

Hence, one possible solution to avoid the wrapping effect and preserving the time invariance of the observer is to used the trajectory-based observer. In this case, the observer is constructed based on the punctual trajectories for estimating upper and lower bounds instead of using a zonotopic set. In Fig. 2, the obtained results from the simulation of the trajectory-based observer for both presented scenarios with the same observer gain as the set-based method is shown. On the one hand, the results of a monotonic system can be seen in Fig. 2a. On the other hand, Fig. 2b shows the results when A contains negative elements. The inner solution problem can be seen in Fig. 2b implying the estimated value by the observer is tighter than the correct one. The problem can be solved by designing the observer gain to eliminate the effect of the negative element of A on A obs . In Fig. 3, the set-based state estimation of interval-observer approach provide wider intervals than those produced by the trajectory-based approach. In other words, the set-based approach even without wrapping effect is more conservative than the trajectory-based approach because of considering the timevarying parametric uncertainty.

On the other hand, the performance of set-based intervalobserver approach and trajectory-based interval-observer approach can be compared by looking at Fig. 2b and Fig. 1b that allows seeing that the wrapping effect is avoided by using trajectory-based observer. But, this approach provides the inner solution instead of the exact solution. The clue to overcome the undesired problem in trajectorybased observer is to design the observer gain L to enforce A obs into an element-wise positive matrix even in the situation that the system is non-monotonic. In Fig. 4, the state estimation is obtained by designing the observer gain according to Section 3 that forces the resulting observer to be monotonic. But, sometimes the positivity condition could be hard to be satisfied or guarantee the convergence of the observer at the same time is not possible. With considering this point in mind, the idea of this paper is to overcome this limitation by designing a reducedorder observer instead of a full-order observer that only estimates the unmeasurable states. Fig. 5 shows the results obtained with the designed reduced-order observer in both scenarios. Moreover, it can be seen in Fig. 5a, in the case that all the elements of the state matrix are positive, the observer is working well. But, the inner solution problem appears in Fig. 5b that indicates the state estimation of the non-monotonic system. That means, the upper and lower approximated bounds is tighter than the exact one. Therefore, the observer gain can be designed to overcome this problem as in Fig. 6, where the obtained result from the simulation that the reduced-order observer gain is designed to solve the inner solution problem is shown. Therefore, the reduced-order observer can be used instead of the full-order observer with more degree of freedom regarding the observer gain design in comparison with the full-order observer. Therefore, the wrapping effect is avoided, the inner solution problem is solved and the observer gain can be designed with more degree of freedom, simultaneously.

CONCLUSIONS

This paper presents the state estimation of the discretetime linear state-space model with the time-invariant parameter uncertainty using different interval-observer approaches. First, it is shown that the time invariance of parameter uncertainty can not be preserved when the observer is designed by using the set-based method and the wrapping effect affects those systems that are not monotonic. Second, the trajectory-based observer is used to overcome the problems. But, the computational complexity is increased when solving the global optimization problem by using the classical trajectory-based observers in the case of the non-monotonic system. Therefore, a method for designing the observer gain is proposed such that a monotonic interval observer is obtained without solving the optimization problems. Moreover, in the case that the number of conditions and constrains for designing the gain does not let to compute the suitable observer gain, a reduced-order observer is proposed to increase the degrees of freedom for selecting the gain. Additionally, it is shown the set-based approach is more conservative than the trajectory-based approach in the case that the system does not affect by the wrapping effect.

As a future research, the idea of forcing the reducedorder observer matrix to satisfy the positivity condition is proposed by using a matrix transformation in order to have much more degrees of freedom to select the observer gain.

Fig. 1 .

 1 Fig. 1. Set-based interval observer.

  A contains a negative element.

Fig. 2 .

 2 Fig. 2. Trajectory-based interval observer.

Fig. 4 .

 4 Fig. 3. Set-based vs. Trajectory-based interval observer.

  A contains a negative element.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Reduced-order interval observer.

This work has been partially funded by the Spanish Government through the project DEOCS (ref. DPI2016-76493-C3-3-R), by MINECO and FEDER through the project HARCRICS (ref. DPI2014-58104-R) and by Agència de Gestió d'Ajuts Universitaris i de Recerca.


M Strugatsky

K Seleznyova

V Zubov

J Kliava
email: janis.kliava@u-bordeaux.fr

New insight in the nature of surface magnetic anisotropy in iron borate

Keywords:

a b s t r a c t

The theory of surface magnetism of iron borate, FeBO 3 , has been extended by taking into consideration a crystal field contribution to the surface magnetic anisotropy energy. For this purpose, a model of distortion of the six-fold oxygen environment of iron ions in the near-surface layer of iron borate has been put forward. The spin Hamiltonian parameters for isolated Fe 3 + ions in the distorted environment of the near-surface layer have been calculated using the Newman's superposition model. The crystal field contribution to the surface magnetic anisotropy energy has been calculated in the framework of the perturbation theory. The model developed allows concluding that the distortions of the iron environment produce a significant crystal field contribution to the surface magnetic anisotropy constant. The results of experimental studies of the surface magnetic anisotropy in iron borate can be described assuming the existence of relative contractions in the near-surface layer of the order of 1 %.

Introduction

Iron borate, FeBO 3 , has rhombohedral calcite type structure with point group symmetry D 3 d and space group D 6 3 d [1] . From the stand point of magnetic properties FeBO 3 is a two sublattice easy plane anti ferromagnet possessing weak ferromagnetism and the Néel temperature T N ≈ 348 K [1] . Within the accuracy of experiments, the sublattice mag netization vectors lie in the basal plane, (0 0 0 1) in hexagonal coordi nate system, and are nearly antiparallel [1] . However, because of a slight departure from anti parallelism, apart from a strong antiferromagnetic moment 𝑳 = 𝕸 1 -𝕸 2 , a weak ferromagnetic moment 𝑴 = 𝕸 𝟏 + 𝕸 𝟐 occurs, 𝕸 𝟏 and 𝕸 𝟐 being the magnetizations of two sublattices with a norm 𝔐 .

Iron borate can be synthesized by two routes: (i) from the solution in the melt and (ii) from the gas phase [2,3] . The gas phase technique allows obtaining bulk single crystals of iron borate with large non basal faces of optical quality. The existence of natural non basal faces has al lowed finding out and describing a specific magnetic state of a thin near surface layer surface magnetism [4,5] caused by lowering of symme try in the environment of near surface magnetic ions in comparison with those in the crystal volume. Néel was the first to specify the existence of surface magnetic anisotropy in ferromagnets due to this effect [6] . How ever, the manifestation of this anisotropy in conventional ferromagnets is usually obscured by the demagnetizing field and large volume mag netocrystalline anisotropy. In contrast, in the weakly ferromagnetic iron borate crystal the surface magnetic anisotropy can be observed because the demagnetizing field, proportional to the total volume magnetiza tion, is small and the anisotropy in the basal plane is weak [1] . There fore, magnetic characteristics of a thin (0.01 0.1 μm) near surface layer of iron borate drastically differ from those of the volume. Such effects were studied both experimentally, using the magneto optic Kerr effect, and theoretically [4,5] .

The surface magnetic anisotropy energy can be calculated as the dif ference of the magnetic energies of ions in the near surface layer and in the crystal volume. This difference is due to two causes: (i) the oc currence of the crystal surface per se (neglecting a modification of the oxygen environment of near surface iron ions) and (ii) structural distor tions arising in the near surface layer.

In FeBO 3 the exchange energy in a good approximation is isotropic [1] ; therefore, the density of the surface magnetic anisotropy energy, 𝜎 is expected to include dipole dipole, 𝜎 dip and crystal field, 𝜎 cf contribu tions. Obviously, the cause (i) intervenes only in the calculation of the dipole dipole contribution whereas the cause (ii) is expected to modify both the dipole dipole contribution as a result of alteration of mutual disposition of iron ions in the near surface layer and the crystal field contribution in as much as the oxygen environment of iron ions in the near surface layer undergoes additional distortions in comparison with the crystal volume.

Previously, a theoretical description of 𝜎 has been put forward tak ing into account only the cause (i), i.e. neglecting structural distortions in the near surface layer [4] . In this approximation, 𝜎 will include only the dipole dipole contribution, and its maximum value will occur for the ( 1 0 1 4 ) face [4] . It was assumed that the magnetic moments of iron ions in the near surface layer lie in the basal plane; thus, the polar an gle ϑ of the reduced antiferromagnetic moment 𝒍 = 1 2 𝑳 ∕ 𝔐 equals 𝜋 2 [4] . Indeed, a considerable deviation of l from this plane would result in the emergence of a transition magnetic layer, similar to a domain wall, be tween the crystal volume and the near surface layer. In this layer the magnetization will be non uniformly distributed, so that the equilib rium orientation of l will gradually turn from that in the volume to that in the near surface layer [4] . As far as the transition layer would possess a very significant energy density; of the order of magnitude of √ 𝑎 ef f 𝐴 ≈ 1 . 8 × 10 -3 J m -2 , ( 𝑎 ef f ≈ 4 . 85 × 10 5 J m -3 being an effective uni axial anisotropy constant in the volume and 𝐴 ≈ 0 . 7 × 10 -11 J m -1 being the constant in the expression describing exchange interaction in the transition layer), the emergence of this layer would be energetically un favorable [4] . Thus, in the near surface layer l will lie in the basal plane or make a small angle with this plane, and the azimuthal angle of l in the dipole dipole approximation will be determined by surface magnetic anisotropy [4] :

𝜎 dip = 𝑎 𝑆 dip sin 2 𝜑 ( 1
)
where 𝜑 is the azimuthal angle of l with respect to the two fold C 2 axis and 𝑎 𝑆 dip is the dipole dipole contribution to the surface magnetic anisotropy constant. At room temperature 𝑎 𝑆 dip = 1 . 4 × 10 -5 J m -2 [5] .

In equilibrium l is parallel to one of the C 2 axes; therefore, the reduced ferromagnetic vector 𝒎 = 1 2 𝑴 ∕ 𝔐 , m ⊥ l , is perpendicular to this axis, so that C 2 is the hard magnetization axis in the basal plane. This result has been confirmed by experimental observations for ( 1 0 1 4 ) face of the crystal [4] .

Applying a magnetizing field H in the basal plane, where the mag netic anisotropy is small, the direction of l in the transition layer will gradually change from that in the volume, determined by H , to that in the near surface layer, determined by both H and the surface magnetic anisotropy. The saturation field along the hard magnetization axis in the near surface layer, called the critical field 𝐻 c , is considered as the measure of the surface magnetic anisotropy [4] :

𝐻 c = 4 𝑎 𝑆 dip 2 AM . ( 2
)
where M is the spontaneous magnetization of the crystal. For FeBO 3 at room temperature 𝑀 = 11 . 79 G [4] . In fact, an application of such a field would totally erase the transition layer. Zubov et al. have found that on the ( 1 0 1 4 ) face of iron borate 𝐻 c ≈ 1 kOe at 300 K; meanwhile, the 𝐻 c value calculated with Eq. ( 2) is much lower, 𝐻 c = 0 . 2 kOe [5] . From the preceding analysis we can ex pect that this discrepancy could be removed if we allow for structural distortions in the near surface layer, ( "surface reconstruction "). Earlier, an attempt in this direction has been made; however, only changes in the positions of iron ions have been considered [5] . (In such an approx imation, only the dipole dipole contribution to the surface magnetic anisotropy is accounted for.) In this case, the experimental 𝐻 c value could be obtained only for relative extensions as large as 7 -12 % [5] .

Meanwhile, it is evident that in the near surface layer, the oxygen environment of iron ions undergoes additional distortions in comparison with the crystal volume, giving rise to the crystal field contribution to the surface magnetic anisotropy energy. The aim of the present work is to consider these distortions and to provide a complete description of the surface magnetic anisotropy of iron borate, including both the dipole dipole and crystal field contributions.

Surface reconstruction

In the crystal volume, each iron ion is surrounded by six oxygens forming a nearly perfect octahedron. In turn, borons are located at the centers of equilateral oxygen triangles, so that Fe 3 + cations turn out to be sixfold coordinated by flat BO 3 3-groups playing the role of anions. EPR studies have revealed the existence of two non equivalent iron sites with local axes rotated through an angle ± 𝛼 about the C 3 axis, see Fig. 1 [7] .

The smallest rhombohedron with faces of the type ( 1 0 1 4 ) , contain ing iron ions in all vertices and face centers (face centered cell), is shown in Fig. 2 ; the faces parallel to the crystal surface are colored. These faces are perpendicular to yz planes. We assume that the distortions of inter atomic distances in the near surface layer occur in yz planes.

In the near surface layer, oxygen octahedrons of two non equivalent sites of Fe 3 + are distorted. Consider the locations of irons and sur rounding oxygens for these two sites with respect to the crystal surface. Fig. 3 (a) shows the rhombohedron of Fig. 2 turned so that the yz plane is parallel to the plane of the figure, while the crystal surface and the basal plane xy are orthogonal to yz plane. (For simplicity, we show only two non equivalent sites, Fe 1 and Fe 2 and their oxygen environments in the absence of distortions.)

We assume that the distortions in the near surface layer occur only in the oxygen octahedrons nearest to the crystal surface denoted by F in Fig. 3 and passing through O 1 6 and O 2 6 oxygens (the subscripts num bering the iron sites and the superscripts numbering the oxygens linked to a given iron site). We can divide the near surface layer into four par allel sub layers separated by planes passing through oxygen and iron ions, see Fig. 3 . We assume that the B plane remains immobile, the dis tortions occur along the AD edge of the rhombohedron, and the dis placements of ions lying at different distances from B are proportional to these distances. Fig. 3 (b) illustrates this model; for clarity the relative displacements have been exaggerated.

Crystal field surface anisotropy

Previously, we have carried out EPR studies of Fe 3 + in GaBO 3 single crystals and determined the spin Hamiltonian parameters using crystal lographic data and the Newman's superposition model [7,8] .

In order to calculate the crystal field contribution to the surface magnetic anisotropy energy, we have considered the generalized spin Hamiltonian expressed by means of tesseral spherical tensor operators 𝒯 𝑙 𝐵 𝑙 𝑆 𝑙𝑚 ( 𝒏 , 𝑺 ) where n is the unit vector in the direction of the magne tizing field H and l B and l S are powers of H and of the spin operators, respectively [9] .

The required form of the spin Hamiltonian can be adapted from Eq. ( 4) in the paper by Tennant et al. [10] , as follows (cf. the review paper by Kliava and Berger [11] ):

ℌ = ∑ 𝑙=2 , 4 𝑙 ∑ 𝑚 =-𝑙 ℬ 0 𝑙 𝑙𝑚 𝒯 𝑙𝑚 ( 𝑺 ) + 𝑔 e 𝛽𝐻 ∑ 𝑙 𝑆 =1 , 3 , 5 ∑ 𝑙= ⎧ ⎪ ⎨ ⎪ ⎩ 𝑙 𝑆 -1 𝑙 𝑆 + 1 𝑙 ∑ 𝑚 =-𝑙 ℬ 1 𝑙 𝑆 𝑙𝑚 𝒯 1 𝑙 𝑆 𝑙𝑚 ( 𝒏 , 𝑺 ) . ( 3
)
Here the first and second terms on the right hand side are, respec tively, zero field ( l B = 0) and linear Zeeman ( l B = 1) spin Hamiltonians; l = | l S -l B |,…, l S + l B must be even to preserve time inversion invariance. In particular, for l B = 0, l = l S , and for l B = 1, l = l S ± 1. The ℬ 0 𝑙 𝑙𝑚 parame ters in Eq. ( 3) are proportional to the corresponding Stevens parameters 𝐵 𝑚 𝑙 , and the procedure of calculating ℬ

1 𝑙 𝑆 𝑙𝑚 is described in detail by Mc Gavin et al. [9] .

As far as ℬ 𝑙 𝐵 𝑙 𝑆 𝑙𝑚 are components of rank l irreducible tensors [9,10] , they can be consistently expressed within the superposition model, as follows [7,11] :

ℬ 𝑙 𝐵 𝑙 𝑆 𝑙𝑚 = 𝑛 ∑ 𝑗=1 𝑏 𝑙 𝐵 𝑙 𝑆 𝑙𝑚 ( 𝑟 𝑗 ) 𝐶 𝑚 𝑙 ( 𝜗 𝑗 , 𝜑 𝑗 ) (4)
where j enumerates the nearest neighbors of the paramagnetic ion with spherical coordinates r j , ϑ j , 𝜙 j , and 𝑏 𝑙 𝐵 𝑙 𝑆 𝑙𝑚 ( 𝑟 𝑗 ) and 𝐶 𝑚 𝑙 ( 𝜗 𝑗 , 𝜑 𝑗 ) are, respec tively, radial functions and tesseral coordination factors; the latter are given in Appendix of the paper by Seleznyova et al. [7] . For simplicity, the radial functions have been chosen in the following form:

𝑏 𝑙 𝐵 𝑙 𝑆 𝑙𝑚 ( 𝑟 𝑗 ) = b 𝑙 𝑆 ( 𝑟 0 ∕ 𝑟 𝑗 ) 𝑡 𝑙 𝑠 ( 5
)
where b 𝑙 𝑆 are intrinsic parameters, 𝑡 𝑙 𝑆 are the corresponding power ex ponents and r 0 = 2.101 Å is the reference distance corresponding to the average Fe O distance in MgO [12] . For near surface iron ions, the su perposition model parameters are expected to remain the same as for those in the volume, at least as far as distortions in their environment remain weak. We have used the following values of b 𝑙 𝑆 and 𝑡 𝑙 𝑆 [7] :

b 2 = 0 . 4 cm -1 , 𝑡 2 = 8 and b 4 = 3 . 1 ⋅ 10 -5 cm -1 , 𝑡 4 = 5 . ( 6
)
The b 3 and b 5 parameters turn out to be very small [7] , therefore, the corresponding terms in the spin Hamiltonian have been neglected.

On the other hand, the spin Hamiltonian for iron ions in the near surface layer should include terms describing the concomitant lowering of symmetry. Besides, as far as Fe 3 + ions in the near surface layer of FeBO 3 are subject to a strong exchange field, we can use the mean field approximation and substitute an effective exchange field H k for the mag netizing field H ( k = 1, 2 numbering the non equivalent iron sites). For the type of distortions described above, the spin Hamiltonian takes the form:

ℌ 𝑘 = 𝑔 𝛽𝑯 𝑘 ⋅ 𝑺 𝑘 ± 𝐵 -2 2 O -2 2 𝑘 + 𝐵 -1 2 O -1 2 𝑘 + 𝐵 0 2 O 0 2 𝑘 ± 𝐵 1 2 O 1 2 𝑘 + 𝐵 2 2 O 2 2 𝑘 + 𝐵 0 4 O 0 4 𝑘 ± 𝐵 3 4 O 3 4 𝑘 + 𝐵 -3 4 O -3 4 𝑘 ( 7
)
where g is close to the free electron g factor g e = 2.0023, 𝛽 is the Bohr magneton, S k is the electron spin for the k th Fe 4 𝑘 are, respectively, the second and fourth order fine structure parameters. In Eq. ( 7) as well as in the subsequent equations the ± signs before certain terms correspond to k = 1, 2, respec tively.

Assuming the "Zeeman " term ℌ 0 = g 𝛽H k • S k in Eq. ( 7) to be much larger than the remaining terms, we have calculated the energies of the spin levels E mk using the perturbation theory. With this aim in view, we have followed the approach of "correct " zero order wave functions described by Landau and Lifshitz [START_REF] Landau | Quantum Mechanics: Non-Relativistic Theory[END_REF] . In the first order of perturbations, we get:

𝐸 𝑚𝑘 = 𝑔𝛽𝑚 𝐻 𝑘 + [ 3 𝑚 2 -𝑆 ( 𝑆 + 1 ) ] ℎ 1 2 𝑘 + [ 35 𝑚 4 -30 𝑚 2 𝑆 ( 𝑆 + 1 ) + 25 𝑚 2 -6 𝑆 ( 𝑆 + 1 ) + 3 𝑆 2 ( 𝑆 + 1 ) 2 ] ℎ 1 4 𝑘 (8)
where 𝑆 = 𝑆 𝑘 = 5 2 , 𝑚 = -5 2 , -3 2 , … , 5 2 and

ℎ 1 2 𝑘 = ± 1 2 𝐵 -2 2 sin 2 𝜗 𝑘 sin 2 𝜑 𝑘 + 1 4 𝐵 -1 2 sin 2 𝜗 𝑘 sin 𝜑 𝑘 + 1 2 𝐵 0 2 ( 3 cos 2 𝜗 𝑘 -1 ) ± 1 4 𝐵 1 2 sin 2 𝜗 𝑘 cos 𝜑 𝑘 + 1 2 𝐵 2 2 sin 2 𝜗 𝑘 cos 2 𝜑 𝑘 ; ℎ 1 4 𝑘 = 1 8 𝐵 -3 4 cos 𝜗 𝑘 sin 3 𝜗 𝑘 sin 3 𝜑 𝑘 + 1 8 𝐵 0 4 ( 35 cos 4 𝜗 𝑘 -30 cos 2 𝜗 𝑘 + 3 ) ± 1 8 𝐵 3 4 cos 𝜗 𝑘 sin 3 𝜗 𝑘 cos 3 𝜑 𝑘 . ( 9
)
In the strong exchange case, the spins of two non equivalent iron ions are antiparallel: 𝜗 1 = 𝜗, 𝜑 1 = 𝜑, 𝜗 2 = 𝜋 -𝜗, and 𝜑 2 = 𝜑 + 𝜋.

At 𝑇 = 0 K , the only occupied spin level is the lowest one (with 𝑚 = - 5 2 ), so that the anisotropic part of the right hand side of Eq. ( 8) for this level provides the crystal field contribution to the density of the magnetic anisotropy energy of the near surface layer:

𝜎 cf 𝑆 = 1 2 𝑁 [ 10 𝐵 -1 2 𝑆 sin 𝜗 cos 𝜗 sin 𝜑 + ( 30 𝐵 0 2 𝑆 -450 𝐵 0 4 𝑆 + 10 𝐵 2 2 𝑆 ) cos 2 𝜗 + 525 𝐵 0 4 𝑆 cos 4 𝜗 +20 𝐵 2 2 𝑆 sin 2 𝜗 cos 2 𝜑 + 15 𝐵 -3 4 𝑆 cos 𝜗 sin 3 𝜗 sin 3 𝜑 ] ( 10
)
where N is the number of Fe 3 + ions per unit surface ( 𝑁 = 6 . 0036 × 10 18 m -2 for ( 1 0 1 4 ) face) and the S subscript refers to the parameters for the In the absence of structural distortions, i.e. at 𝜀 = 0, the crystal field contribution to the density of the magnetic anisotropy energy of the near surface layer would be the same as for an analogous layer in the crystal volume: where the V subscript refers to the parameter in the volume, and N is the same as in Eq. ( 10) . Thus, 𝜎 cf , defined as 𝜎 c f S -𝜎 c f V , can be expressed as:

𝜎 cf 𝑉 = 1 2 𝑁 [( 30
𝜎 cf = 𝜎 cf 1 sin 𝜗 cos 𝜗 sin 𝜑 + 𝜎 cf 2 cos 2 𝜗 + 𝜎 cf 3 cos 4 𝜗 + 𝜎 cf 4 sin 2 𝜗 cos 2 𝜑 + 𝜎 cf 5 cos 𝜗 sin 3 𝜗 sin 3 𝜑 ( 12
)
where 𝜎 cf 1 , 𝜎 cf 2 , 𝜎 cf 3 , 𝜎 cf 4 and 𝜎 cf 5 are defined as the difference between the parameters of matching symmetry terms in Eqs. (10) and (11) .

Comparison with experimental data

The parameters featuring in Eq. ( 12) depend on relative distortions and have been calculated for | 𝜀 | ⩽3%. The crystal field contribution to the surface magnetic anisotropy energy is a few orders less than the en ergy required for emergence of the transition layer between the crystal volume and the near surface layer. Thus, the emergence of this layer is still energetically unfavorable, see Introduction, so that 𝜗 = 𝜋 2 [4] . So that taking into account Eqs. (10) (12) the crystal field contribution to the surface magnetic anisotropy energy can be expressed as:

𝜎 cf = -10 𝑁𝐵 2 2 𝑆 sin 2 𝜑 ( 13
)
Comparing Eqs. ( 13) and ( 1) , we get the crystal field contribution to the surface magnetic anisotropy constant:

𝑎 𝑆 cf = -10 𝑁𝐵 2 2 𝑆 , ( 14
)
and the total surface magnetic anisotropy constant can be expressed as:

𝑎 𝑆 = 𝑎 𝑆 dip + 𝑎 𝑆 cf . ( 15
)
The 𝑎 𝑆 dip ( 𝜀 ) dependence for the layer containing Fe 3 + ions (D in Fig. 3 ) has been previously calculated at 300 K [5] . As far as the dipole dipole interaction energy is proportional to the square of magnetization, 𝑎 𝑆 dip is expected to depend on the temperature as 𝔐 ( T ) 2 , 𝔐 ( T ) being the sublattice magnetization at temperature T . Taking into account that 𝔐 (0 K)∕ 𝔐 (300 K) ≈ 1 . 47 [START_REF] Kadomtseva | [END_REF] , we have calculated 𝑎 𝑆 dip ( 𝜀 ) at 0 K. Fig. 4 shows the dependences on 𝜀 of 𝑎 𝑆 dip [5] , 𝑎 𝑆 cf and 𝑎 𝑆 at 0 K. Obviously, the distortions much more affect 𝑎 𝑆 cf than 𝑎 𝑆 dip . According to experimental data [4] , in the absence of H , in equilibrium l is parallel to C 2 ; thus, 𝑎 𝑆 is positive, which corresponds to the case of contractions. Indeed, for the near surface layer of ionic crystals contractions rather than extensions are expected [16] .

Next, from the 𝑎 𝑆 ( 𝜀 ) dependence shown in Fig. 4 , we calculate 𝐻 c ( 𝜀 ) , see Fig. 5 , using the equation similar to Eq. (2) :

𝐻 c = 4 𝑎 𝑆 2 AM ( 16
)
The experimental value of 𝐻 c at 300 K, (1 kOe, vide supra ) and the temperature dependence of 𝐻 c , see Fig. 10 in the paper by Zubov et al. [4] , suggest that at 0 K 𝐻 c ≈ 1 . 47 kOe. As one can see from Fig. 5 , such 𝐻 c value can be attained at a contraction of ca. 1 %.

Note that the crystal surface in itself constitutes a structure defect; nevertheless, other types of structure defects can occur in the near surface layer, for instance, vacancies of magnetic ions. Such defects are also expected to contribute to 𝑎 𝑆 cf . As model objects for studying surface magnetism in crystals with such defects one can consider diamagneti cally diluted iron gallium borates, Fe x Ga 1 -x BO , in which case 𝑎 𝑆 cf will depend on x . In Fig. 6 we show the calculated 𝑎 𝑆 cf ( 𝑥 ) dependence. Ex perimental studies of surface magnetism in these crystals are in progress in our laboratories.

Conclusions

The theory of the surface magnetism of iron borate has been ex tended to include, besides the dipole dipole contribution, the crystal field contribution to the surface magnetic anisotropy energy.

A model of structural distortions in the near surface layer for ( 1 0 1 4 ) face of iron borate has been developed assuming that ions located in this layer are displaced proportionally to their distance from a reference plane assumed to remain immobile. In order to account for the lowering of symmetry of the ions in the near surface layer, the generalized spin Hamiltonian expressed through tesseral spherical tensor operators has been applied, and the spin Hamiltonian parameters have been calculated using the Newman's superposition model. The crystal field contribution to the density of the surface magnetic anisotropy energy at 0 K has been calculated in perturbation theory. The distortions of iron environment produce a significant crystal field contribution to the surface magnetic anisotropy constant a S ; indeed, the experimental results on the studies of surface magnetism by the magneto optic Kerr effect can be satisfactorily described assuming relative contractions in the near surface layer of ca. 1 %.

Fig. 1 .

 1 Fig. 1. Two non-equivalent sites of Fe 3 + in iron borate. The Cartesian coordinate axes are directed as follows: 𝑥 ∥ C 2 , y lies in the symmetry plane m and 𝑧 ∥ С 3 , the three-fold axis [1] . The z -axis is perpendicular to the plane of the figure and points toward the reader. Full and empty circles represent ions located above and below this plane, respectively.

Fig. 2 .

 2 Fig. 2. A rhombohedron with edge length a r = 5.862 Å and apex angle 104.2°used to calculate the density of surface magnetic anisotropy energy. The Cartesian coordinate axes are directed as in Fig. 1 .

Fig. 3 .

 3 Fig. 3. Two non-equivalent Fe 3 + sites and their oxygen environments: without distortions (a) and in the case of contraction (b) (full red and empty blue circles for Fe 1 and Fe 2 , respectively). The Cartesian coordinate axes are directed as in Fig. 1 . The x -axis is perpendicular to the plane of the figure and points towards the reader. The nearest-neighboring "top " and "bottom " oxygen triangles for both iron sites have the same z -coordinate and different x -coordinates. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .

 4 Fig. 4. Values of 𝑎 𝑆 cf , 𝑎 𝑆 dip [5] and a S vs. 𝜀 at 0 K.

Fig. 5 .

 5 Fig. 5. Calculated H c ( 𝜀 ) dependence at 0 K.

Fig. 6 .

 6 Fig. 6. Calculated dependence of 𝑎 𝑆 cf on x for 𝜀 = 1% in Fe x Ga 1 -x BO 3 crystals at 0 K.

Acknowledgments

This work was partially supported by the Russian Foundation for Basic Research and the Ministry of Education, Science and Youth of the Republic of Crimea, in the framework of scientific project grant no. 15 42 01008 "р юг а " and by the V.I. Vernadsky Crimean Federal Univer sity Development Program for 2015 2024 including the participation of one of the authors (K. Seleznyova) in the European School on Magnetism 2017 (Corsica, France).