domain: spi math info phys

text:
Jean-Marc Biannic
email: jean-marc.biannic@onera.fr

Clément Roos
email: clement.roos@onera.fr

Samir Bennani
email: samir.bennani@esa.int

Fabrice Boquet
email: fabrice.boquet@esa.int

Valentin Preda
email: valentin.preda@esa.int

Bénédicte Girouart
email: benedicte.girouart@esa.int

Advanced probabilistic µ-analysis techniques for AOCS validation

Keywords:

Monte-Carlo simulations play a key role in the current Attitude and Orbit Control Systems (AOCS) Verification and Validation (V&V) process, but it is generally time-consuming and it may fail in detecting worst-case configurations, especially in the presence of rare events. In such a case, µ-analysis offers a nice alternative, although it cannot measure the probability of occurrence of the identified worst-cases, which can invalidate a control system on the basis of unlikely events. Probabilistic µ-analysis was introduced in this context 20 years ago to bridge the gap between the two techniques, but until recently no practical tools were available. This paper summarizes recent advances on this topic with a particular emphasis on practical applications to space systems. More precisely, the proposed technique is applied to evaluate AOCS controllers in the context of a challenging high accuracy satellite pointing control problem. The way the proposed tools can be integrated into the traditional AOCS V&V process and used to tighten the V&V analysis gap is also highlighted.

Probabilistic µ-analysis, Monte-Carlo simulations, Space systems

Introduction

Novel lightweight materials and deployable structures allow to perform an increasingly wide variety of on-orbit services, which results in stronger interactions between the spacecraft flexible structures (solar arrays, robotic arms, antennas, mirrors) and fuel sloshing in the reservoirs. Efficient AOCS must therefore be designed to ensure a high pointing accuracy. This requires a robust control architecture, but also adequate V&V methods to assess the mission risk and check if performance is guaranteed regardless of the uncertainties and external disturbances affecting the system.

Monte-Carlo simulations [START_REF] Helton | Survey of samplingbased methods for uncertainty and sensitivity analysis[END_REF][START_REF] Landau | A guide to Monte Carlo simulations in statistical physics[END_REF] are the preferred validation means in the space industry. They are able to quantify the probability of sufficiently frequent phenomena, but they are generally time-consuming, provide only soft bounds [START_REF] Zhu | Soft vs. hard bounds in probabilistic robustness analysis[END_REF] and may fail in detecting rare but nevertheless critical events. Less expensive deterministic and simulation-free alternatives exist and have reached a good level of maturity, as is the case for µ-analysis [START_REF] Zhou | Robust and optimal control[END_REF][START_REF] Ferreres | A practical approach to robustness analysis with aeronautical applications[END_REF]6]. But unlike Monte Carlo simulations, if worst-case scenarios are no longer missed, their probability of occurrence is also not measured, which can invalidate an AOCS on the basis of very rare and therefore extremely unlikely events [7,[START_REF] Marcos | Stochastic µ-analysis for launcher thrust vector control systems[END_REF][START_REF] Thai | Probabilistic µ-analysis for stability and H ∞ performance verification[END_REF].

Research to fill the gap between these two approaches is still at a very early stage and only few practical tools are available, although this issue was identified 20 years ago by [START_REF] Khatri | Guaranteed bounds for probabilistic µ[END_REF]. This is all the more surprising since the validation process currently accounts for up to 80% of the AOCS total development time, and is becoming longer as the space missions become increasingly complex. In this context, this paper builds on the work of [START_REF] Thai | Probabilistic µ-analysis for stability and H ∞ performance verification[END_REF][START_REF] Falcoz | Probabilistic µ-analysis for system performances assessment[END_REF]7,[START_REF] Zhu | Improved bounds computation for probabilistic µ[END_REF] on probabilistic µ-analysis and its ambitions are twofold. First, develop new cheap and reliable tools to improve the characterization of rare but nonetheless possible events, so as to tighten the aforementioned V&V analysis gap. Second, apply these tools to a challenging high accuracy satellite pointing control problem, to show how they can be integrated into the traditional AOCS V&V cycle to improve the current industrial standard and fasten the validation process.

The paper is organized as follows. The latest advances in probabilistic µanalysis are first presented in Section 2. The resulting computational tool is then described in Section 3. It is finally applied in Section 5 to a challenging and realistic AOCS benchmark introduced in Section 4.

Latest advances in probabilistic µ

Problem statement

Let us consider a continuous-time uncertain linear time-invariant system (usually including control laws):

ẋ = A(δ)x + B(δ)u y = C(δ)x + D(δ)u (1)
where the parametric uncertainties δ = (δ 1 , . . . , δ N ) ∈ R N are independent random variables with probability density functions f = (f 1 , . . . , f N ). It is assumed that A(δ), B(δ), C(δ), D(δ) are polynomial or rational functions of the δ i . As a result, system (1) can be transformed into a Linear Fractional Representation (LFR) as in Figure 1 (right): the uncertainties are separated from the nominal (closed-loop) system M (s) and isolated in a block-diagonal operator ∆ = diag(δ 1 I n 1 , . . . , δ N I n N ), where I n i is the n i × n i identity matrix.

Note that this paper focuses on real parametric uncertainties, but as for classical µ-analysis, complex uncertainties and neglected dynamics can be considered as well. The set of matrices with the same block-diagonal structure as ∆ is denoted ∆. It is assumed that the uncertainties are normalized, so that the subset B ∆ of ∆ defined as

B ∆ = {∆ ∈ ∆ : σ(∆) < 1} = {∆ ∈ ∆ : |δ i | < 1, i ∈ [1, N ]} corresponds to the set of physically meaningful uncertainties B δ = [-1 1] N .
With these notations in mind, two main problems can be solved using probabilistic µ-analysis: Problem 2.1 (Probabilistic robust stability). Compute the probability P ∆,f (M (s)) that the interconnection of Figure 1

(left) is unstable when ∆ ∈ B ∆ .
Problem 2.2 (Probabilistic worst-case H ∞ performance). Given a performance level γ > 0, compute the probability P γ ∆,f (M (s)) that T u→y (s, ∆) ∞ > γ on Figure 1 (right) when ∆ ∈ B ∆ , where T u→y (s, ∆) is the transfer from u to y.

Once computed, these probabilities can be confronted to a given tolerance level , so as to validate or reject the considered control system, depending on whether P ∆,f (M (s)) and P γ ∆,f (M (s)) are lower or higher than .

Remark 2.1. The uncertainties being bounded, their probability distributions must be supported on a bounded interval. Uniform and truncated normal distributions are often used in practice.

Probabilistic robust stability

Classical µ-analysis [START_REF] Zhou | Robust and optimal control[END_REF][START_REF] Ferreres | A practical approach to robustness analysis with aeronautical applications[END_REF] aims at computing the robust stability margin k r , which satisfies the following properties:

• the interconnection of Figure 1 (left) is stable for all ∆ ∈ k r B ∆ = {∆ ∈ ∆ : |δ i | < k r , i ∈ [1, N ]},
• for all k > k r , there exists at least one ∆ ∈ kB ∆ such that the interconnection is unstable.

This defines an uncertainty box centered at 0 and of radius k r , which touches the instability domain and where stability is guaranteed. So if k r < 1, there are parts of the uncertainty domain B δ that can be stable or unstable, but for which nothing can be concluded. Let us illustrate this with a simple example extracted from [START_REF] Falcoz | Probabilistic µ-analysis for system performances assessment[END_REF]. The considered system is described by the state-space representation:

     ẋ = 0 1 -a 1 (δ 1 ) -a 2 (δ 2 ) x + 0 1 u y = 1 0 x (2)
where a 1 (δ 1 ), a 2 (δ 2 ) are two uncertain parameters defined as:

a 1 (δ 1 ) = 1 + 2δ 1 a 2 (δ 2 ) = 0.8 + δ 2 (3)
and δ 1 , δ 2 are two normalized real parametric uncertainties. In this particular case, the stability and instability domains can be calculated analytically, using e.g. the Routh-Hurwitz stability criterion. They are represented in light/dark green and dark red respectively in Figure 2. The domain k r B δ where stability is guaranteed by µ-analysis is the light green box, and it is clear that there are both stable and unstable zones outside this area, where no information is available at this stage. Remark 2.2. Computing k r is in general NP-hard, so lower and upper bounds are computed instead. Much work has been done to reduce the gap between these bounds, and many efficient algorithms are now available [START_REF] Roos | A detailed comparative analysis of all practical algorithms to compute lower bounds on the structured singular value[END_REF]. It can thus usually be assumed that the (almost) exact value of k r can be computed with a reasonable computational time.

A branch-and-bound (B&B) algorithm can be used to explore the whole uncertainty domain. The idea is to partition B δ into smaller boxes until each box has guaranteed stability or is sufficiently small to be neglected [START_REF] Newlin | Mixed µ problems and branch and bound techniques[END_REF]. Taking the union of the boxes with guaranteed stability leads to an approximation D s ⊂ B δ of the exact stability domain. In practice, this approximation is usually quite accurate, as can be seen in the previous example by comparing the green areas in Figures 2 and3. This strategy can however result in a significant computational time for large-scale problems. In particular, the boxes where the uncertain system is unstable are never identified as such by the algorithm, since µ-analysis requires to start with a nominally stable system. They are partitioned until they reach a negligible size, which leads to an unnecessary generation of boxes. This is materialized in Figure 3 by the use of the blue color to represent these very small undetermined boxes. But this issue can be addressed easily. µ-analysis indeed detects when the poles of the nominally stable system M (s) reach the imaginary axis as the size of ∆ increases. The same strategy can be applied to a nominally unstable system. Integrated into the previous B&B scheme, this directly yields a domain D s ⊂ B δ of guaranteed instability. The domain D s of guaranteed stability conveniently takes the form of a finite union of disjoint boxes D (k) s :

D s = k D (k) s where D (k) s = [δ (k) 1 , δ (k) 1 ] × • • • × [δ (k) N , δ (k) N ]
making the associated probability p(D s ) easy to compute:

p(D s ) = k p(D (k) s ) (4) = k N i=1 δ (k) i δ (k) i f i (δ i )dδ i (5)
The same applies to D s , which finally leads to both lower and upper bounds on the exact probability P ∆,f (M (s)) of instability, thus solving Problem 2.1:

p(D s ) ≤ P ∆,f (M (s)) ≤ 1 -p(D s ) (6)
The considered control system can then be either validated if 1-p(D s ) < or rejected if p(D s ) > , where is the tolerance level introduced in Section 2.1.

In practice, B&B is executed until the gap between the bounds becomes small enough and one of the two previous conditions occurs. The uncertainty domain is finally partitioned as follows:

B δ = D s ∪ D s ∪ D su (7)
where D su denotes the domain of undetermined stability, with probability p(D su ). The B&B algorithm can indeed only approximate D s and D s , and not compute them exactly. The application to the previous academic example is shown in Figure 4, where it can be seen that very good approximations of the exact stability and instability domains are obtained.

Remark 2.3. When applying B&B, the boxes can be divided along the direction with the highest µ-sensitivity (i.e. corresponding to the uncertainty with the greatest influence on stability, see [START_REF] Braatz | µ-sensitivities as an aid for robust identification[END_REF]), instead of the longest edge as usually done. This significantly reduces the number of boxes and the computational time in most cases [START_REF] Thai | Probabilistic µ-analysis for stability and H ∞ performance verification[END_REF][START_REF] Lesprier | Improved µ upper bound computation using the µ-sensitivities[END_REF]. This is observed here by comparing the number of green boxes in Figures 3 and4, obtained without and with the µ-sensitivities respectively.

Probabilistic worst-case H ∞ performance

As with robust stability, B&B is combined with µ-analysis to compute domains of guaranteed performance D γ and guaranteed non-performance D γ , as well as the associated probability measures p(D γ ) and p(D γ ). This leads to bounds on the exact probability P γ ∆,f (M (s)) of non-performance, thus solving Problem 2.2:

p(D γ ) ≤ P γ ∆,f (M (s)) ≤ p(D s ) -p(D γ ) (8)
Performance is guaranteed on a given box

D, i.e. D ⊂ D γ , if: max ∆∈D T u→y (s, ∆) ∞ ≤ γ (9)
This can be reformulated as a skew µ problem using the main loop theorem [START_REF] Packard | The complex structured singular value[END_REF] and solved using existing µ-based tools [6]. On the other hand, checking if non-performance is guaranteed on D, i.e. if D ⊂ D γ , requires to solve: min

∆∈D T u→y (s, ∆) ∞ > γ (10)
This is a minimax problem, which cannot be directly reformulated as a skewµ problem as above. To address this issue, a sufficient condition for inequality [START_REF] Khatri | Guaranteed bounds for probabilistic µ[END_REF] to hold is introduced in [START_REF] Thai | Probabilistic µ-analysis for stability and H ∞ performance verification[END_REF], in the form of a skew µ calculation involving the inverse transfer T -1 u→y (s, ∆). It is very efficient from a computational point of view, but it may introduce conservatism, although this is usually not the case in practice. The main limitation is that this condition is restricted to a scalar performance channel, i.e. u ∈ R and y ∈ R in [START_REF] Helton | Survey of samplingbased methods for uncertainty and sensitivity analysis[END_REF].

A B&B algorithm can then be implemented as for stability. The investigated domain is limited to the domain of guaranteed stability D s , since performance analysis only makes sense for stable systems. This leads to:

D s = D γ ∪ D γ ∪ D γu ( 11
)
where D γu is the domain of undetermined performance. The following partition of B δ is finally obtained by combining ( 7) and ( 11):

B δ = D γ ∪ D γ ∪ D γu ∪ D s ∪ D su (12)

A new computational tool

All the results and algorithms presented in Section 2 have been implemented in the Matlab function mupb, the main routine of the new STOWAT 1 (STOchastic Worst Case Analysis Tools) package to be integrated in a forthcoming version of the SMAC (Systems Modeling Analysis and Control) Toolbox2 . This function solves Problems 2.1 and 2.2 by computing guaranteed lower and upper bounds on the probabilities P ∆,f (M (s)) and P γ ∆,f (M (s)) with the desired accuracy. It is fully compatible with the Generalized State-Space (GSS) Library of the SMAC Toolbox [START_REF] Biannic | Generalized State Space: a new Matlab class to model uncertain and nonlinear systems as LFR[END_REF], which offers a powerful and user-friendly way to model LFR, including the ability to define probability distributions for parametric uncertainties. This library is itself compatible with the standard uss object of the Robust Control Toolbox, and it provides automated tools to convert uss objects to gss objects, as well as to incorporate probability distributions initially not present in uss objects. Finally, the function mupb is also fully interfaced with the Skew Mu Analysis Robustness Tools (SMART) Library of the SMAC Toolbox [6], which allows the user to benefit from a number of state-of-the-art µ-analysis based algorithms. It can be called quite easily as follows:

[pbnds,partout]=mupb(sys,pb,partin,options); where:

• sys is a gss or uss object describing the LFR of Figure 1, • pb specifies whether Problem 2.1 or 2.2 is to be solved,

• partin can be used to provide an initial description of D s , D s , D su , D γ , D γ and D γu coming from a previous call to mupb, • options contains tuning parameters related to the stopping criterion of the B&B algorithm, the use of µ-sensitivities and the accuracy of the µ computation, • pbnds gives the guaranteed bounds on P ∆,f (M (s)) or

P γ ∆,f (M (s)),
• partout provides the final list of boxes which make up D s , D s , D su , D γ , D γ and D γu (same format as partin).

Benchmark description

As already emphasized in the introduction, pointing performance is increasingly demanding on both scientific and observation space missions of future generations [START_REF] Sanfedino | Dynamic modeling and analysis of micro-vibration jitter of a spacecraft with solar arrays drive mechanism for control purposes[END_REF]. In the presence of uncertainties, performance degradation remains very challenging to quantify reliably as it results from complex interactions between external perturbations and structural flexible modes of the spacecraft at very specific frequencies. Inspired by previous works presented in [START_REF] Sanfedino | Dynamic modeling and analysis of micro-vibration jitter of a spacecraft with solar arrays drive mechanism for control purposes[END_REF][START_REF] Preda | Robust microvibration mitigation and pointing performance analysis for high stability spacecraft[END_REF][START_REF] Boquet | Active solutions to the spacecraft jitter problem[END_REF], the proposed benchmark focuses on the effects on pointing accuracy of the micro-perturbations induced by a Solar Array Drive Mechanism (SADM).

Model description

The system under consideration, illustrated by Figure 5, is essentially composed of a main body, two solar arrays SA 1 and SA 2 , an isolated payload P L, and a wheel W whose mass and inertia are neglected here.

The single axis case is considered in this work. Let us denote X = [θ x SA 1 x SA 2 x P L ] the position vector of the plant, where θ is the pointing error. Using the M -D -K formalism, the model is described as:

Ẋ Ẍ = 0 I -M -1 K -M -1 D X Ẋ + M -1 Z Γ B Γ SA (13)
D = S(V D ) K = S(V K ) (14)
where:

S(V ) = sum(V ) -V -V diag(V ) (15)
and:

V D = [D SA 1 D SA 2 D P L ] = 10 2 × [0.5 0.5 8] V K = [K SA 1 K SA 2 K P L ] = 10 5 × [1 1 200] (16)
Finally, the input matrix:

Z = 1 0 0 0 0 1 1 0 ( 17
)
distributes the control torque Γ B applied to the main body and the input perturbation torques Γ SA = Γ SA 1 = Γ SA 2 which, for simplicity, are assumed to affect similarly the two solar arrays.

Uncertainties mainly affect the first three elements of the generalized inertia matrix as follows:

     J B = (1 + 0.1 δ 0 )J B 0 J SA 1 = (1 + 0.2 δ 1 )J SA 10 J SA 2 = (1 + 0.2 δ 2 )J SA 20 (18)
where the δ i -denoting normalized real parametric uncertainties -respectively introduce 10% variations on the main body and 20% on each solar array. Resulting from these uncertainties, rather significant variations (see Table 1) can be observed on the main flexible modes respectively located near 6.8 rad/s and 12.5 rad/s in the nominal case. Such variations require a robust control design, which is described below. The particular structure of equation ( 13) where the uncertainties, through M -1 , enter the model in a rational way, strongly suggests to use the LFR framework. Either using uss or gss objects [START_REF] Biannic | Generalized State Space: a new Matlab class to model uncertain and nonlinear systems as LFR[END_REF], a minimal LFR F u (G(s), ∆), with ∆ = diag(δ 0 , δ 1 , δ 2 ), is readily obtained and can be integrated in a robust H ∞ control design scheme (see Figure 6). In this scheme, A(s) = (1+0.05s) -1 denotes the actuator model, while the W X (s) are the standard input/output weighting functions of the H ∞ framework. Here, the most specific ones are W SA (s) and W AP E (s), which shape the disturbance inputs interacting with the flexible modes of the solar arrays and the absolute pointing error (APE) output respectively. Based on previous work introducing relevant metrics for pointing accuracy [START_REF] Pittelkau | Definition, metrics and algorithms for displacement, jitter and stability[END_REF][START_REF] Ott | ESA pointing error engineering handbook[END_REF][START_REF] Pittelkau | Optical transfer functions, weighting functions, and metrics for images with two-dimensional line-of-sight motion[END_REF], it appears that the latter can be quantified via the H ∞ norm of a weighted transfer. In this application, the proposed weighting functions, adapted from [START_REF] Boquet | Active solutions to the spacecraft jitter problem[END_REF], are all described in Table 2. The main functions of interest (W SA and W AP E ) have been calibrated in such a way that the pointing performance requirement is met when:

Element

γ ape = T w SA →z AP E (s) ∞ ≤ 1 ( 19
)
Note that W SA (s) is a very selective band-pass filter designed to excite the system around 12.45 rad/s, while W AP E (s) is a high-pass filter.

Name

Numerical expression

Main characteristic

W SA W SA (s) = 0.01245 s 155+0.01245 s+s 2 2 selective band-pass W AP E W AP E (s) = 2 × 10 5 × 1+300s 1+3s high-pass W U I W U I (s) = 10 -4 × 1 1+0.1s low-pass W Y W Y (s) = 10 -7 × 1 1+2s low-pass W U O W U O (s) = 10 × 1+0.002s 1+0.8s
low-pass

Controller design

Control design is performed in two steps. A preliminary PD-like first order controller K P D (s) is initially designed by a modal-based approach that robustly stabilizes the double integrator of the plant for all admissible uncertainties δ i . In a second step, the poor performance level of this initial controller is improved by a multi-model structured H ∞ design approach based on systune and a general strategy exposed in [START_REF] Apkarian | Multi-model, multi-objective tuning of fixed-structure controllers[END_REF]. Note however that in the particular case of this challenging application, the optimization problem is not only nonsmooth but also strongly nonconvex when a low order is required. Then, despite its high efficiency, the optimization algorithm implemented in systune generally fails or requires a large number of random starts. This issue is solved by providing the solver with a non-minimal third-order stabilizing initial condition K IN IT (s) directly derived from K P D (s) as follows:

K IN IT (s) = K P D (s) × s 2 + 2 ξ 0 ω 0 s + ω 2 0 s 2 + 2 ξ 0 ω 0 s + ω 2 0 =1 = K P D (s) ( 20
)
where ω 0 = 12.5 rad/s and ξ 0 = 0.001 are chosen close to the frequency and damping parameters of the band-pass filter W SA (s). Interestingly, by this approach, a very low (third) order controller is rapidly obtained despite the relatively high order of the weighted design interconnection (n = 17). As is visible on the Bode plot of Figure 7, this controller exhibits poorly damped poles and zeros near 12.5 rad/s and 13 rad/s respectively. These are the necessary ingredients to ensure a good rejection of the input perturbations on the solar arrays. However, with such a low-order solution, it was not possible to ensure a priori that the performance requirement ( 19) is rigorously met for all admissible uncertainties.

Preliminary robustness analysis

A preliminary analysis is realized to evaluate the robustness of the pointing performance. This is first achieved by a standard evaluation of the H ∞ norm of T w SA →z AP E (s) for 2000 randomly generated configurations according to a uniform distribution for each uncertainty δ i . The results are presented in Figure 8 which also displays a worst-case plot detected by a skew µ-analysis approach. At this stage, the controller seems to be validated by the standard Monte-Carlo approach since the H ∞ norm never gets larger than 1, while it is invalidated by µ-analysis. A refined but also much more time consuming3 analysis involving from N = 5000 to 200000 samples is then performed. As can be seen in Table 3, a minimum number of 50000 samples is required to obtain a reasonably tight approximation PN (γ ape > 1) ≈ 10 -4 of the probability of failure P F = P(γ ape > 1), which as already proven by skew µ-analysis is not zero. Remark 4.1. As recalled in [START_REF] Tempo | Probabilistic robustness analysis: explicit bounds for the minimum number of samples[END_REF], hard bounds on the probability of failure P F cannot be obtained by this approach. However, given tolerance levels and δ, the Chernoff bound enables to predict a minimum number of samples N to ensure a certain accuracy in a probabilistic setting:

Method

N > 1 2 2 log 2 δ ⇒ P |P F -PN | > < δ (21)
But this bound is very conservative here. Assuming for example that = δ = 10 -4 , one actually get N ≥ 5 × 10 8 ! A less conservative point of view is proposed in [START_REF] Tempo | Probabilistic robustness analysis: explicit bounds for the minimum number of samples[END_REF] from which, after minor adaptations to our context, one obtains:

N > log δ log(1 -) ⇒ P (P(γ ape > γN ) > ) < δ (22)
Let us choose N = 15000 so that γN ≈ 1 (see Table 3) and set δ = . From [START_REF] Pittelkau | Definition, metrics and algorithms for displacement, jitter and stability[END_REF], it is then readily checked that = (1 -) N ⇒ ≈ 5 × 10 -4 . This can be summarized as:

P P F > 5 × 10 -4 < 5 × 10 -4 (23)
Interestingly this result is obtained after "only" 15000 simulations and, as is the case with the Chernoff bound [START_REF] Boquet | Active solutions to the spacecraft jitter problem[END_REF], this number of required simulations is independent of the number of uncertain parameters. Unfortunately, equation [START_REF] Pittelkau | Definition, metrics and algorithms for displacement, jitter and stability[END_REF] does no permit to evaluate P F independently of N which must be fixed here such that γN ≈ 1. As a result, the tolerance cannot be freely defined either. However, it provides a good insight into the required number of simulations to be performed, which is confirmed by the numerical experiments summarized in Table 3. ever, as is clarified next, this will not systematically have a huge impact on the probability level to be demonstrated.

Remark 5.1. In this application involving 3 uncertain parameters, the CPU time required by the proposed µ-based approach is very competitive when compared to Monte-Carlo techniques and provides (beyond guaranteed results) much more information. But the computational effort can significantly increase with the number of uncertain parameters, although the latter is not necessarily a good measure of the computational burden. A model with many uncertain parameters, but only a few of which have a large influence on stability or performance, will indeed be easier to analyze than a model with 5 or 10 uncertain parameters of equivalent importance. This is where the µ-sensitivities introduced in Remark 2.3 play a key role. That being said, although there is no general rule as explained above, numerical difficulties are likely to occur beyond 5 to 10 parameters. A promising approach to push this limit is to exploit parallel computing. This is easily achieved for Monte-Carlo techniques, and turns out to be slightly more complex (but feasible!) for the proposed µ-based approach. The proposed algorithm involves a list of domains, each of which can be analyzed independently of the others. We can therefore consider the following implementation: each available CPU first takes the domain with the highest probability from the list, then performs the analysis, and in case (in)stability or (non-)performance cannot be guaranteed finally divides the domain and returns the resulting subdomains to the list. This important issue clearly deserves further investigations and will be the subject of future works.

Probability levels cheap reassessment

Beyond the guaranteed results (through hard bounds on the probability of failure) provided by probabilistic µ-analysis, the most interesting feature is that the computed sets of domains (illustrated in Figures 9 and10) can be reused at a very low cost to evaluate new probability levels for any distribution of the uncertain parameters in their respective intervals. In this application, this is achieved by considering truncated normal distributions with varying standard deviations σ ∈ [0.3 , 5] for each uncertainty. The results, obtained in a few seconds without recalling the analysis routine mupb, are presented in Figures 11 and12, which show the evolution of the guaranteed upper and lower bounds on P(γ ape > 1) as a function of σ. As expected, the gap between upper and lower bounds is reduced in the refined case (Figure 12). Moreover, a convergence towards the bounds presented in Table 4 (obtained with uniform distributions and visualized by the black solid and dashed lines) is observed when σ → ∞. Conversely, for small values of σ, the probability levels become unsurprisingly so small that their evaluation through Monte-Carlo strategies would require a prohibitive number of simulations. In this respect, the proposed approach clearly bridges a gap between classical methods and µ-analysis. Moreover, the obtained results permit to establish an explicit and constructive relationship between the probability to be demonstrated and the characteristics of the uncertainties. As an example, if the targeted probability of failure is 10 -6 , then the controller is validated as long as σ ≤ 0.44. Finally, a comparison of the upper bounds (see Figure 13) between the default and refined settings reveals a reasonably small gap -especially for low σ -which might not justify the additional CPU time. The refined analysis has essentially permitted to increase the lower bound on the probability of failure. As a result, for the given probability level to be demonstrated the controller is validated by the default analysis as long as σ ≤ 0.42. The difference between the two results based on default or refined settings is then not really significant. Relevant conclusions can thus be obtained here at a very low cost (a few seconds).

CONCLUSIONS AND PERSPECTIVES

Based on recent algorithmic advances in probabilistic µ-analysis and a new Matlab implementation, promising performance assessment results have been obtained on a realistic AOCS benchmark. This opens the way towards the integration of probabilistic µ in a more efficient design-and-validation process. A proposed integration scheme to be further evaluated through future studies is exposed in Figure 14. As emphasized in Section 5.2, a very interesting feature of the proposed approach is to enable very fast iterations between the probability level and the uncertainties characterization. At a medium cost, uncertainty resizing is also possible. Future work will be devoted to more advanced performance evaluation involving several criteria simultaneously such as H ∞ performance on multiple channels but also ro- [7] G. Balas, P. Seiler, A. Packard, Analysis of an UAV flight control system using probabilistic µ, in: Proceedings of the AIAA GNC Conference, 2012.

Figure 1 :

 1 Figure 1: Standard interconnections for robust stability (left) and worst-case performance (right) analysis

Figure 2 :

 2 Figure 2: Exact stability domain (light/dark green), exact instability domain (dark red), and guaranteed stability domain obtained with µ-analysis (light green)

Figure 3 :

 3 Figure 3: Guaranteed stability domain obtained with µ-analysis and B&B (light green), and domain of undetermined stability (blue)

Figure 4 :

 4 Figure 4: Domains of guaranteed stability D s (light green), guaranteed instability D s (light red) and undetermined stability D su (blue) obtained with µ-analysis and B&B

Figure 5 :M

 5 Figure 5: Simplified view of the nominal plant

Figure 6 :

 6 Figure 6: Weighted LFT-based closed-loop model

Figure 7 :

 7 Figure 7: Controllers Bode diagrams

Figure 8 :

 8 Figure 8: Robust performance: Monte-Carlo vs skew µ-analysis

Figure 9 :Figure 10 :

 910 Figure 9: Cross-sectional view (δ 0 = 1) of the performance regions (default settings)

Figure 11 :

 11 Figure 11: Bounds on P(γ ape > 1) w.r.t standard deviation based on the initial set of performance regions (default settings)

  bound (normal) Lower bound (normal) Upper bound (uniform) Lower bound (uniform) Probability threshold =0.4

Figure 12 :

 12 Figure 12: Bounds on P(γ ape > 1) w.r.t standard deviation based on the refined set of performance regions (refined settings)

  upper bound (CPU time = 2 s) Refined upper bound (CPU time = 10 s) Probability threshold

Figure 13 :

 13 Figure 13: Comparison of upper-bounds on P(γ ape > 1) w.r.t standard deviation

Figure 14 :

 14 Figure 14: Integration of probabilistic µ-analysis in a design & validation process

Table 1 :

 1 Flexible modes characteristics

	Damping	Frequency (rad/s)
	Mode #1 [1.6 1.8] × 10 -3	[6.35 7.39]
	Mode #2 [3.0 3.3] × 10 -3	[11.9 13.1]
	Mode #3 [4.7 4.8] × 10 -3	[231 238]

Table 2 :

 2 Weighting functions

Table 3 :

 3 # samples (N ) PN (γ ape > 1) Worst case (γ N ) CPU Refined Monte-Carlo (MC) robustness analysis

	MC	5000	0	0.76	20 s
	MC	15000	4.9 × 10 -4	1.00	60 s
	MC	25000	2.4 × 10 -4	1.34	100 s
	MC	50000	1 × 10 -4	1.76	190 s
	MC	100000	9 × 10 -5	1.89	410 s
	MC	200000	9.5 × 10 -5	2.23	860 s
	skew µ	NA	NA	γ ape = γ∞ = 3.08	6 s

These tools were developed under ESA contract RFP/3-16071/19/NL/CRS/hh

The SMAC toolbox is available at w3.onera.fr/smac

All computations have been performed without parallelization on a standard laptop equipped with a processor Intel i5-8400H, 2.50GHz with 16Gb RAM installed.

This work is the result of the research project ESA/RFP/3-16071/19/NL/CRS/hh

Application of probabilistic µ-analysis

In this example involving a limited number of uncertain parameters and a priori a very low probability of failure of the control system, Monte-Carlo techniques are time-consuming and certainly not the best suited. On the other hand, µ-analysis can quickly prove non-performance but does not provide any associated risk probability. To fill this gap, the new probabilistic µ tool presented in Section 3 is then used to address Problem 2.2 as detailed in Section 2.3.

Initial characterization of the performance regions

Considering first uniform distributions for all parametric uncertainties δ i , the routine mupb is called on the uncertain (but robustly stable) closedloop LFR described in Section 4 to evaluate the performance regions D γ , D γ and D γu , from which guaranteed probability bounds are deduced. A first analysis is performed with default settings (first row in Table 4). Next, the algorithm is resumed with a refined stopping criterion to obtain tighter (but still guaranteed) bounds on the probability of failure P F while keeping a reasonable computational time (second row in Table 4). Note that in both cases, the computed intervals confirm the previous Monte-Carlo analysis (with 50000 samples), which is recalled on the last row of Table 4. Crosssectional visualizations (with δ 0 = 1) of the performance regions in the plane < δ 1 × δ 2 > are then presented in Figures 9 and 10 for both settings. Figure 10 clearly shows that the undetermined region (visualized in blue) is significantly reduced in the second case using some refined settings. How-

Method


Lamy Jean-Baptiste

Abdelmalek Mouazer

Karima Sedki

Rosy Tsopra

Translating the Observational Medical Outcomes Partnership -Common Data Model (OMOP-CDM) electronic health records to an OWL ontology

Keywords: Medical Records, Electronic Health Records, Biological Ontologies, SPARQL

The heterogeneity of electronic health records model is a major problem: it is necessary to gather data from various models for clinical research, but also for clinical decision support. The Observational Medical Outcomes Partnership -Common Data Model (OMOP-CDM) has emerged as a standard model for structuring health records populated from various other sources. This model is proposed as a relational database schema. However, in the field of decision support, formal ontologies are commonly used. In this paper, we propose a translation of OMOP-CDM into an ontology, and we explore the utility of the semantic web for structuring EHR in a clinical decision support perspective, and the use of the SPARQL language for querying health records. The resulting ontology is available online.

Introduction

Electronic health records (EHR) [START_REF] Kataria | Electronic health records: a critical appraisal of strengths and limitations[END_REF] lead to a major progress in the storage, the transmission and the standardization of clinical patient data. However, today, many EHR models and formats exist, each software vendor proposing its own. This heterogeneity is a huge problem for research studies that need to collect data from many EHR, but also for clinical decision support systems that need to be interfaced with many different EHR.

In the last ten years, OMOP-CDM (Observational Medical Outcomes Partnership -Common Data Model) [START_REF] Reich | OMOP Common Data Model Specifications[END_REF] from the OHDSI (Observational Health Data Sciences and Informatics) community emerged as a common and simple EHR model, used to structure clinical data extracted from various other EHR, in order to facilitate clinical research studies. This model is proposed as a relational database schema. A study showed that OMOP-CDM has a higher content coverage than three other similar data models [START_REF] Garza | Evaluating common data models for use with a longitudinal community registry[END_REF]. More recently, OMOP-CDM has been considered for clinical decision support [START_REF] Unberath | EHR-Independent Predictive Decision Support Architecture Based on OMOP[END_REF][START_REF] Gruendner | KETOS: Clinical decision support and machine learning as a service -A training and deployment platform based on Docker, OMOP-CDM, and FHIR Web Services[END_REF].

In parallel, formal ontologies and the semantic web [START_REF] Schulz | Formal ontologies in biomedical knowledge representation[END_REF] have emerged as a standard for the formalization of medical knowledge. Ontologies permit formal reasoning but also facilitate the reuse of the data and the knowledge. In particular, the interest of ontologies for validating EHR models have been shown in the literature [START_REF] Martínez-Costa | Validating EHR clinical models using ontology patterns[END_REF], and ontologies are commonly used in decision support systems.

In this paper, we propose an OWL translation of the OMOP-CDM model, and we explore the utility of the semantic web for structuring EHR in a clinical decision support perspective. Our objective is not to maintain a full compatibility with OMOP-CDM database, but rather to structure an EHR as a formal ontology, grounding on the experience of OMOP-CDM. Consequently, we will focus on the clinical part of OMOP-CDM, and we will not consider the vocabulary part, because ontologies offer native support for structuring hierarchical terminologies. We will also consider the use of the SPARQL language for querying health records, and compare it to the SQL language.

Material and methods

Material

We used OMOP-CDM version 6.0 [START_REF] Reich | OMOP Common Data Model Specifications[END_REF]. In OMOP-CDM (Figure 1), patients and healthy volunteers are represented by the Person table. Each Person may have zero, one or several Visit Occurrence, e.g. visits to a GP or hospital stays. Each Visit may be associated with some diagnoses (Condition Occurrence), tests (Measurement), medical procedures (Procedure Occurrence), drug prescriptions (Drug Exposure), etc. A higher level of abstraction, Eras, is also provided, for facilitating epidemiological studies. An Era groups one or more similar time periods in a single entity; e.g. if a patient was prescribed metformine for 3 months, and then after 3 months, metformine was prescribed again, there are two Drug Exposures (one per prescription), but a single Drug Era. OMOP-CDM provides procedures for computing Eras from the Drug Exposures and Condition Occurrences. Both Eras and lowerlevel entities (Condition Occurrence, Drug Exposure,...) are associated with a concept from a medical terminology. We used the Python programming language for parsing OMOP-CDM specification and generating the OWL ontology, with the Owlready ontology-oriented programming module [START_REF] Lamy | Owlready: Ontology-oriented programming in Python with automatic classification and high level constructs for biomedical ontologies[END_REF][START_REF] Lamy | Ontologies with Python[END_REF][START_REF] Lamy | Ontology-Oriented Programming for Biomedical Informatics[END_REF].

Translating the database model to OWL

We translated the OMOP-CDM database model into an OWL ontology, using a automatic Python script. Each table was translated into a class, each field corresponding to an identifier into an object property, and each non-identifier field into a data property. SQL datatypes were translated into XML Schema datatypes and assigned to the range of data properties. We also added universal class restrictions, and existential class restrictions for fields marked as required in the OMOP-CDM model. The translation may seem rather straightforward, but two difficulties were encountered.

First, contrary to SQL relational databases, OWL ontologies support inheritance. In the OMOP-CDM model (Figure 1), there are some obvious situations where inheritance could be used, e.g. both Person and Provider are humans, and share identical attributes, such as gender or year of birth. Similarly, many OMOP-CDM tables have attributes related to time: the date of an event (e.g. a Condition Occurrence) or the start and end date of a duration (e.g. a Drug Exposure).

We manually added superclasses, such as Base Person, Event or Duration. Then, we wrote an automatic procedure in Python that moves any attribute present in all the subclasses of a superclass to that superclass. For instance, since Base Person has two subclasses, Person and Provider, and since both Person and Provider have a gender attribute, then the gender attribute will be moved to the Base Person class.

Second, in the OMOP-CDM model, the direction of relation is often dictated by the relational database model, e.g. the relation between Person and Drug Era is in the Drug Era Person direction, through the "person_id" field, because each Drug Era is associated with a single Person (*-1 relation) while a Person may be associated with several Drug Era (1-* relation), and the relational model supports only *-1 relations. On the contrary, OWL ontologies support both types of relations. Here, we found more intuitive to reverse the relation, in order to start from the Person and then to obtain his Drug Era via the has_drug_era relation. This places the patient at the center of the model. Moreover, in the database model, fields are defined within a given table, e.g. the "person_id" field in the Drug Era table is distinct from the "person_id" field in the Condition Era table, despite they share the same name. This permits to search for the Drug Era associated with a given patient, using the "per-son_id" field in the Drug Era table. On the contrary, in ontologies, properties are first-order entities, independent from classes. Thus, if we create a "has_person" object property (corresponding to "person_id" but named using ontology conventions), it will be shared by all Era classes. If we search for the Drug Era associated with a given patient, the "has_person" property will thus return all Era, and thus we need to select only those that are Drug Era. In a SPARQL query, this requires two RDF triple patterns: SELECT ?drug_era { ?drug_era a DrugEra . ?drug_era has_person <patient_x> . } If we reverse the direction of the relation, leading to the "has_drug_era" property, we can now select specifically the Drug Era associated with a given patient, simplifying the query as follows:

SELECT ?drug_era { <patient_x> has_drug_era ?drug_era . } Therefore, we manually defined a list of relation to reverse, and then an automatic procedure in Python for reversing those relations.

Importing data

We imported the sample dataset with 1,000 patients from CMS SynPUF, proposed on the OMOP-CDM website. When importing the data, terminological concepts were mapped to UMLS. We used the UMLS import functions from Owlready2 to extract the terminologies present in OMOP from UMLS version 2020AA, and to translate them to ontologies. In these ontologies, concepts are represented by classes, in order to allow inheritance between concepts. When a concept is present in the OMOP data, the corresponding class is instantiated. For instance, if a patient has a Condition Occurrence associated with the "59621000" SNOMED CT code for Essential hypertension, we create an instance of the "59621000" class and we associate it with the Condition Occurrence. The original class can be obtained via the rdf:type relation between the instance and the class.

Querying the ontology

OMOP-CDM data are usually accessed via SQL queries, while SPARQL is commonly used for ontologies. We compared the length and the complexity of the queries in both languages, using various queries inspired by those proposed on the OMOP website (http://cdmqueries.omop.org/), or by the recommendations of the STOPP/START guidelines [START_REF] O'mahony | STOPP/START criteria for potentially inappropriate prescribing in older people: version 2[END_REF] for detecting potentially inappropriate prescribing in the elderly.

Person Drug Era has_drug_era

Results

OMOP-CDM ontology translation

Figures 2 and3 show the general model of the ontology. Compared to Figure 1, notice the presence of inheritance, but also the fact that the direction of many relations was changed. The ontology consistency was verified using the Pellet reasoners. It contains 49 classes, 226 properties and 5,016 RDF triples. After importing the OMOP-CDM sample data, we obtained a total of 10,795,221 RDF triples, including 3,311,359 for terminologies.

The ontology is available online 1 (Apache License 2.0, same license as OMOP-CDM), as well as the Python scripts that generated it from the OMOP-CDM CSV specifications 2 (GNU LGPL license). The script can be customized through global variables, e.g. to export only some parts of OMOP-CDM.

Using SPARQL for querying health records

Figure 4 shows an example of query, in both SQL and SPARQL, inspired by the examples found in the OMOP website. The SQL query is longer and more complex, due to the presence of nested queries. The first nested query (SELECT DISTINCT condition.person_id...) is motivated by the fact that "age" appears three times in the outer query (in the SELECT, the GROUP BY and the ORDER BY clauses). But SQL does not allow creating variables and to assign value to them. The nested query is thus here to avoid duplicating three 1 http://www.lesfleursdunormal.fr/static/_downloads/ omop_cdm_v6.owl 2 https://bitbucket.org/jibalamy/owlready2/src/master/ pymedtermino2/omop_cdm/ times the formula for computing the "age" value. On the contrary, SPARQL allows the creations of variable, using the BIND statement, hence removing the need for the nested query.

The second nested query (SELECT DISTINCT descendant_concept_id...) is used to select all descendants of the desired concept in the terminology. Thanks to property path expressions, SPARQL offers an easier way to select descendants, e.g. using the expression "rdfs:subClassOf*", where "*" means that zero, one or more subClassOf relations must exist between a descendant concept and the original concept. In addition, "/" can be used in property path expressions to chain several relations, e.g. "has_concept/a/rdfs:subClassOf*/rdfs:label "Fracture of bone of hip region"" means that ?condition has for concept an instance of a class that is a descendant of a class associated with the "OMOP Hip Fracture 1" label.

Consequently, with SPARQL, nested queries are no longer required here. This arguably simplifies the query.

Using Owlready, the SPARQL query can be executed on the OMOP sample dataset (1,000 patients, about 10 million RDF triples) in about 0.31 second on a recent laptop computer. In contrast, the SQL query can be executed on the same dataset in about 0.35 second using PostgresQL. This suggests that SPARQL and ontologies are as efficient as SQL and relational databases.

The semantic web also allows linking the data with other data or knowledge very easily. For instance, the sample OMOP dataset use RXNORM for coding drug prescriptions, but the ATC terminology (Anatomical, Therapeutical Chemical classification of drugs) may be more practical when dealing with therapeutical classes (e.g. proton pump inhibitors, PPI, instead of specific active principles such as omeprazole), especially when implementing the rules found in clinical guidelines. With ontologies, it is easy to add relations with ATC in addition to the existing relations with RXNORM: when importing the data, we instantiated the RXNORM classes; we can state that the resulting instances also belong to the corresponding ATC classes by adding new RDF triples. On the contrary, when using database, one would require to add an extra table mapping RXNORM to ATC for adding support for the ATC terminology. This would complicate queries, with additional joints between the concept table and the mapping table.

Figure 5 shows two examples of rules extracted from the STOPP/START clinical guideline [START_REF] O'mahony | STOPP/START criteria for potentially inappropriate prescribing in older people: version 2[END_REF], and their implementation in SPARQL. We used SNOMED CT codes for disorders and ATC codes for drugs. The first rule detects a simple drugdisorder interaction between digoxin and heart failure. The second rule is more complex. It involves aspirin, which has 3 ATC codes; thus, we used a UNION clause for testing the 3 codes. Moreover, the rule should not be triggered when a PPI is prescribed concomitantly. Thus, we used a FILTER NOT EXISTS clause to verify the absence of a PPI, with conditions on the start and end dates of the two Drug Eras to verify the co-occurrence of the two treatments. Both rules were implemented as SELECT queries that returns the patients and the Drug Era that should be stopped.

When implementing STOPP/START rules, we found that SPARQL lacks the IN SQL keyword. This keyword allows testing whether a value is one of a set of given values. It is a shorthand for multiple OR conditions. IN is commonly used when several terms correspond to the desired concept, e.g. for testing the 3 codes for aspirin in the ATC terminology, one may use "aspirin IN ("BO1AC06", "AO1AD05", "NO1BA01")". On the contrary, SPARQL has no such keyword. Thus, in the same situation, we used UNION as seen above, which is much longer and less practical.

Finally, notice that the precision of the queries remains the same with SQL or SPARQL, both queries being semantically equivalent and returning the same results.

Discussion

OMOP-CDM was initially developed for gathering in a single model clinical data from heterogeneous sources, such as EHR from different vendors, in order to facilitate clinical research. However, clinical decision support is another situation where one may need to merge clinical data from different models, e.g. for a given patient, clinical data may be found in the EHR of the GP, but also in the EHR of the hospital and even in the pharmacy. In that situation, the use of OMOP-CDM seems a relevant option. Moreover, the high-level abstractions proposed in OMOP-CDM with Era are useful for clinical studies, but also for decision support. For instance, Drug Eras allow computing drug treatment durations, which are sometimes required for supporting decision (e.g. for rule START E2 in STOPP/START).

The use of ontologies for structuring EHR may not be relevant for clinical research, because of the huge volume of data in-volved, and the relational database format is well established today in that community. On the contrary, for decision support, the volume of data is often lower (after the eventual learning phase for machine learning-based systems), because decision support deals with a single patient at a time, and ontologies are frequently used.

In the literature, a previous tentative exists for translating the OMOP-CDM model to an ontology [START_REF] Pacaci | A Semantic Transformation Methodology for the Secondary Use of Observational Healthcare Data in Postmarketing Safety Studies[END_REF]. However, it does not focus on the latest version of the model (6.0), and the translation was limited to a raw conversion from database to OWL, without adding inheritance and restrictions as we did, nor reversing relations. Finally, since our translation is almost entirely automatic, and performed by a Python script, it will be easy to update the ontology for future versions of OMOP-CDM.

In the near future, we plan to use the proposed ontology for structuring heterogeneous clinical data in a decision support system for medication reviews. We also plan to implement tools for importing into the proposed ontology clinical data in standard formats such as HL7 and FHIR.

Conclusions

In this paper, we proposed an OWL translation of the OMOP-CDM relational database model for electronic health records. We successfully used the resulting ontology for importing the OMOP sample dataset. We also compared the use of the SQL and SPARQL language for querying EHR data. We found that SPARQL often permitted simpler queries, thanks to its ability to deal with recursion and to define variables, and thanks to the ease with which ontologies can be enriched and connected to other resources in the semantic web.

1 Figure 1 -

 11 Figure 1-UML diagram showing the main tables and relations in OMOP-CDM.

Figure 2 -

 2 Figure 2-UML diagram showing the main classes and relations in the ontology translation of OMOP-CDM.

Figure 3 -

 3 Figure 3-Class hierarchy in the ontology translation of OMOP-CDM.

Acknowledgements

This work was funded by the French Research Agency (ANR) through the ABiMed project [grant number ANR-20-CE19-0017-02].

Address for correspondence

Jean-Baptiste Lamy <jean-baptiste.lamy@univ-paris13.


Abdelmalek Mouazer
email: malikmouazer@gmail.com

Karima Sedki

Rosy Tsopra

Jean-Baptiste Lamy

Visual Comparison of Guidelines: Method and Application to Potentially Inappropriate Medication Lists

Keywords: Set visualization, Therapeutic guidelines, Potentially inappropriate medications, Drug classification, ATC classification

Therapeutic guidelines developed by experts are essential tools for improving therapy and drug prescription. Several guidelines often exist that target the same patient, from different organizations and countries. The case of lists for the detection of potentially inappropriate medications (PIMs) is an example which illustrates how these guidelines can be varied and multiple. In order to have an overview to the divergences and similarities between different lists of PIMs, we propose a visual method to compare PIMs lists, based on set visualization, and we apply it to 5 guidelines.

Introduction

Therapeutic guidelines include recommendations for optimizing patient therapy, inspired by evidence-based medicine. However, several guidelines with the same objective can be available at the same time from different institutions and/or countries, and this can lead to potential conflicts between recommendations. The comparison of several guidelines is a complex and tedious task. In the literature, few approaches have been proposed for comparing guidelines [START_REF] Galopin | Using an ontological modeling to evaluate the consistency of clinical practice guidelines: application to the comparison of three guidelines on the management of adult hypertension[END_REF][START_REF] Seroussi | Using the Recommendations Inferred by a Decision Support System to Compare Breast Cancer Clinical Practice Guidelines[END_REF], and they remain limited to two or three guidelines only. Regarding the detection of potentially inappropriate medications (PIMs), several guidelines have been published. The first and one of the best known is the BEERS List [START_REF] Beers | Explicit criteria for determining inappropriate medication use in nursing home residents. UCLA Division of Geriatric Medicine[END_REF], published in its first version in 1991 in the United States, which focuses on American practice. Inspired by this one, several guidelines have been developed in different countries including criteria more suited to local practice. These guidelines are presented in the form of a list of criteria. Each criterion is presented as: either a drug potentially inappropriate in all situations in the elderly, or a drug potentially inappropriate under certain conditions. These conditions may be an underlying disorder (contraindication), a high dose, or an associated drug (association not recommended). An example from Beers criteria is: "History of gastric or duodenal ulcers with aspirin>325 mg/d or COX-2 non selective NSAIDs." This criterion indicates that aspirin at a dose greater than 325 mg per day as well as non-selective COX-2 NSAIDs are inappropriate in elderly patients with a history of gastric or duodenal ulcer. The objective of this work is to propose a method for the quantitative comparison of PIMs guidelines. The comparison will focus on the drugs mentioned in PIMs guidelines. It is aimed at facilitating the selection of the most consensual guidelines for future implementation in a decision support system. We will apply this method to the comparison of 5 PIMs guidelines: 4 European guidelines, LAROCHE [START_REF] Laroche | Potentially inappropriate medications in the elderly: A French consensus panel list[END_REF], PRISCUS [START_REF] Holt | Potentially Inappropriate Medications in the Elderly: The PRISCUS List[END_REF], GHEOPS [START_REF] Tommelein | Older patients' prescriptions screening in the community pharmacy: development of the Ghent Older People's Prescriptions community Pharmacy Screening (GheOP3S) tool[END_REF], STOPP (Subset of stopp & start criteria version 2 [START_REF] O'mahony | STOPP/START criteria for potentially inappropriate prescribing in older people: version 2[END_REF]) and the BEERS List version 2019 [START_REF]By the 2019 American Geriatrics Society Beers Criteria® Update Expert Panel. American Geriatrics Society 2019 Updated AGS Beers Criteria® for Potentially Inappropriate Medication Use in Older Adults: 2019 AGS BEERS CRITERIA® UPDATE EXPERT PANEL[END_REF], which remains an international reference and the first list to be published.

Method

Extraction and standardization of potentially inappropriate drugs

Our guidelines are lists of criteria in textual format containing the drug descriptions and additional information relating to the conditions for applying the criterion. Standardization of the drug descriptions is necessary in order to be able to compare them. We focused on the targeted drugs in a criterion, disregarding the eventual conditions under which the criteria were applied. We used automatic natural language processing (NLP) tools for extracting the ATC codes (Anatomical Therapeutic Chemical Classification System) associated with each criterion. We then refined this result by manually adding the missed codes. Drugs are expressed at various levels of granularity in PIMs guidelines, e.g., the aforementioned criterion targets both a molecule, aspirin, and a drug class, COX-2 non selective NSAIDs. The ATC classification has a 5-level hierarchical structure. We expressed all drugs with one or more codes of the fifth level. Higher-level codes representing drug classes were replaced by the fifth-level codes of their descendants.

During this step, each molecule is associated with a single ATC code representing the molecule alone, without taking into account the association codes that exist in ATC. For molecules that have more than one ATC code even without associations, we have kept only one, ensuring that it is the same for all the guidelines. Some drugs (e.g., ginkgo biloba) have no ATC code. So we assigned them an arbitrary code.

.Visualization

Comparing guidelines in terms of the inappropriate drugs they contain may seem simple. However, this is only true when the number of guidelines is very small. A given drug may belong to one or more guidelines. Thus, if we have n guidelines, we have 2 n -1 sets of drugs to consider, one per each non-empty subsets of the set of all guidelines, called intersections. As the number of intersections increases exponentially with n, the comparison rapidly becomes very complex. Sophisticated visual techniques are thus required for presenting intersections visually. Since our data is in the form of sets, the appropriate method should be a set visualization approach. Here, we chose the RainBio approach, developed by Lamy & al. [START_REF] Lamy | RainBio: Proportional visualization of large sets in biology[END_REF]. Compared to the well-known Venn diagram, Rainbio has the advantage to allow proportional visualization. In this approach, it is a matter of visually highlighting the membership relationships of a set to one or more elements. In our case, we consider PIMs guidelines to be sets and ATC codes as elements. For this purpose, we have created a data set with the syntax imposed by RainBio and which contains the data extracted from the PIMs guidelines.

Results

989 drugs were identified as potentially inappropriate in the elderly by at least one of our guidelines. Each guideline is associated with the drugs it detects as inappropriate. A drug can be detected as inappropriate by one or more guides which represents 2 5 -1=31 possible intersections between one or more guidelines. All 31 non-empty intersections were present, which suggests a great interconnection between the guidelines. We visualized these intersections using RainBio (Figure 1). Now when comparing the 5 guidelines, we can : (a) The BEERS and STOPP guides cover a larger number of drugs compared to the other guidelines. (b) There are 37 drugs for which there is an overall consensus as to their inappropriateness (red box at the bottom), this roughly represents only 4% of the total number of drugs and 5% to 20% of the total number of drugs in each guideline. This suggests some discrepancies between the guides when viewed together. (c) The dark blue boxes at the top show the specific drugs for each guideline. 175 and 213 of drugs present in STOPP and BEERS are specific, compared with 27, 29 and 15 drugs for GHEOPS, PRISCUS and LAROCHE. Furthermore, if we take a closer look at these last 3 guides, we notice that almost all of their content that is not specific to them is present in BEERS and/or STOPP (90% to 94% of their content). This suggests that BEERS and STOPP have a high of drugs that are absent from all other guidelines, and thus might be in excess. (d) 362 drugs are shared between BEERS and STOPP (sum of all the boxes that cross the two columns). This suggests that these two guidelines share some similarities, but which still have their specificity. (e) The STOPP guideline shares 452 drugs with the other guides (all of the boxes that cross the STOPP column), it is the guide that shares the most drugs with one or more other guides. It has been placed in the center of the figure by Rainbio's algorithm, which shows its central position. (f) While STOPP, BEERS and GHEOPS share two large yellow boxes, LAROCHE and PRISCUS are each excluded from one of these boxes, thus missing some drugs considered as inappropriate by all other guidelines. To conclude, STOPP and BEERS have a lot of specific drugs, and LAROCHE and PRISCUS miss a lot of drugs present in all other guidelines, so GHEOPS seems to be the more consensual one. Additionally, if we are interested in drugs, RainBio makes it easy to identify drugs in a given intersection by clicking on the corresponding box. guidelines by columns ordered by similarity (each is placed next to the ones with whom it shares the most drugs). The intersections between the guidelines (i.e. shared drugs) are represented by rectangular boxes. The color of the boxes indicates the degree of the intersection, it is a color gradient going from dark blue (intersection with a single guideline; the drugs placed in these boxes are thus specific to that guideline) to red (intersection between all guidelines, drugs in this box are present in all PIMs guidelines). The height of the box is proportional to the number of drugs in the intersection. The boxes are finally organized vertically by degree of intersection and therefore also by color, so we find the dark blue boxes at the top and the red boxes at the bottom.

Discussion & Conclusions

We proposed a visual method for the comparison of therapeutic guidelines. This method makes it easy to compare multiple PIMs guidelines. It allows for the exploration of their divergences and similarities. The comparison appears to be intuitive and very easy to understand, we have a clear and proportional overview. The results observed with this method are consistent with the development logic of the guidelines. Indeed, GHEOPS is based on 18 guidelines including those taken into account in this work, which explains its consensual nature. STOPP and BEERS also take into account drug interactions, contraindications and in particular those related to the renal function, which may explain the fact that they have a lot of specific drugs. That said, STOPP and BEERS are quite different concerning their specific drugs; this might be explained by the fact that BEERS is developed for American practice when STOPP is developed for European practice. The similarities observed are not to be taken as being absolute similarities, a drug may be deemed inappropriate without conditions in one guide while being associated with a condition in another guide. The comparison we propose should not lead to a judgment of the quality of one guideline compared to another. Our approach remains focused solely on a comparison based on the drugs to which the various guides are concerned. Comparing PIMs guidelines on the basis of the number of drugs they are interested in have some limits. Indeed, a guideline may contain many drugs rarely prescribed, or not available on a given market. There are also guidelines that offer alternatives while others do not. Taking these parameters into account seems relevant for an in-depth analysis and qualitative comparison. However, the quantitative comparison we proposed here is complementary to the usual qualitative comparison. Our perspectives include the improvement of the comparison by considering drug interactions and patient conditions present in some guidelines. We are also thinking of improving our method by adding a filter on medicines available on a European market. We also intend to broaden our comparison by adding other guidelines found in the literature. Finally, our method can be extended to be applied to the comparison of clinical practice guidelines for the management of certain diseases. These guidelines provide recommendations for treating a disease based on the patient profile. We could therefore consider the patient profiles as elements and compare the guidelines on the basis of the patient profiles they support; or rather compare them on the basis of the recommended treatments by taking the recommendations as elements. In the future, our method will be implemented into a decision support system for helping health professionals to develop or revise guidelines and thus to improve care quality.

Figure 1 .

 1 Figure 1. Visualizing Potentially Inappropriate Drugs From 5 Guidelines Using RainBio. RainBio representsguidelines by columns ordered by similarity (each is placed next to the ones with whom it shares the most drugs). The intersections between the guidelines (i.e. shared drugs) are represented by rectangular boxes. The color of the boxes indicates the degree of the intersection, it is a color gradient going from dark blue (intersection with a single guideline; the drugs placed in these boxes are thus specific to that guideline) to red (intersection between all guidelines, drugs in this box are present in all PIMs guidelines). The height of the box is proportional to the number of drugs in the intersection. The boxes are finally organized vertically by degree of intersection and therefore also by color, so we find the dark blue boxes at the top and the red boxes at the bottom.

Acknowledgments

This work was partly funded by the French Research Agency (ANR) through the ABiMed project [grant number ANR-20-CE19-0017-02].