Mathieu Ritou

Farouk Belkadi

Zakaria Yahouni

Catherine Da Cunha

Florent Laroche

Benoit Furet

CIRP Template v4

Keywords: Decision Making, Machining, Knowledge based system

In the context of Industry 4.0, data management is a key point for decision aid approaches. Large amounts of manufacturing digital data are collected on the shop floor. Their analysis can then require a large amount of computing power. The Big Data issue can be solved by aggregation, generating smart and meaningful data. This paper presents a new knowledge-based multi-level aggregation strategy to support decision making. Manufacturing knowledge is used at each level to design the monitoring criteria or aggregation operators. The proposed approach has been implemented as a demonstrator and successfully applied to a real machining database from the aeronautic industry.

Introduction

In modern factories, manufacturing activities generate large volumes of digital data, which are partially stored and remain underexploited. The efficient use of manufacturing data is a key leverage point for the management and the improvement of the company performance, notably through a decision-aid system.

Such systems are nowadays based on the analysis of a large volume of heterogeneous data. Such an analysis is impossible manually, so automatic analyses by Data Mining (DM) are required. It consists in using specific algorithms in order to discover interesting hidden information (patterns) in data. The main DM goals are descriptive (e.g. clustering) or predictive (behaviour in the future). More generally, approaches of Knowledge Discovery in Database (KDD) are used, where DM is a particular step in this process. The KDD steps are the data collection, selection, pre-processing, transformation, data mining and, lastly, the interpretation and visualization of the results [START_REF] Fayyad | From data mining to knowledge discovery in databases[END_REF]. This leads to new and useful manufacturing knowledge. KDD approaches are developed in many domains and particularly in manufacturing [START_REF] Harding | Datamining in manufacturing: A review[END_REF]. As in prognosis techniques, KDD approaches can be classified into three categories: purely data-driven (e.g. statistics or Artificial Intelligence technics), model-based or hybrid [START_REF] Gao | Cloud enabled Prognosis for Manufacturing[END_REF]. More generally, a model-based approach can be considered as a kind of knowledge-based approach. It can be based on mechanical or empirical models. It can also integrate other types of expert knowledge. To satisfy reliability and interpretability of results, which are the two key points in industry, a combined data-driven and knowledge-based approach is advised. It has been proven that knowledge integration enables a better consideration of expert expectations and of the business context [START_REF] Laroche | Accessing enterprise knowledge: a context-based approach[END_REF]. To do so, existing works addressing the issue of manufacturing knowledge capitalization and reuse can be considered [START_REF] Ammar-Khodja | Processing Knowledge to Support Knowledge-based Engineering Systems Specification[END_REF][START_REF] Verhagen | A critical review of Knowledge-Based Engineering: An identification of research challenges[END_REF].

As in every manufacturing domain, there is a growing interest for data mining in the machining industry. Global approaches are especially suitable for flexible production systems. Abundant literature deals with the instrumentation and monitoring of machine-tools [START_REF] Teti | Advanced monitoring of machining operations[END_REF]. Interesting data can also be collected from the CNC (Computer Numerical Control), through communication protocols such as OPC-UA [START_REF] Monostori | Cyber-physical systems in manufacturing[END_REF]. Approaches are proposed for the analysis of machining data, based on technics of Machine Learning (ML) and signal processing. The objectives are generally limited to classical monitoring issues at machine level (e.g. estimation of tool wear or surface roughness) [START_REF] Kim | Smart Machining Process Using Machine Learning: A Review and Perspective on Machining Industry[END_REF]. Only a few recent works have expressed interest in global management aspects, such as an application for adaptive scheduling or process improvement [START_REF] Lenz | Holistic approach to machine tool data analytics[END_REF][START_REF] Mourtzis | A cloud-based cyber-physical system for adaptive shop-floor scheduling and condition-based maintenance[END_REF][START_REF] Godreau | Continuous improvement of HSM process by data mining[END_REF]. Furthermore, most of DM or ML works are applied to a database collected during a cutting test campaign in laboratories. Due to the varied situations encountered in industrial flexible production, the complexity of DM is higher, making knowledge-based approach particularly suitable for it.

The analysis of a large volume of manufacturing data leads to Big Data challenges. Indeed, the KDD approaches generally consist in direct and successive computations of the different steps, requiring huge computing power, especially for the DM step [START_REF] Tao | Data-driven smart manufacturing[END_REF]. It can rapidly become computationally inefficient to analyse such a large volume of data. One solution is to split the computational problem into smaller ones. Distributed and parallel computing can be envisaged, especially in a context of cloud computing [START_REF] Gao | Cloud enabled Prognosis for Manufacturing[END_REF]. Another solution is to reduce the volume of data. Technics of dimension reduction (e.g. principal component analysis) are proposed. Data aggregation can be performed at the transformation step of the KDD process. In this way, a smaller volume of simplified and more meaningful data is obtained. Moreover, rather than a unique aggregation, a multi-level approach (i.e. granular computing) is a good solution to address the Big Data problem [START_REF] Chen | Data-intensive applications, challenges, techniques and technologies: a survey on big data[END_REF].

The paper presents an original knowledge-based multi-level aggregation approach for data mining in machining. Manufacturing knowledge is integrated at different aggregation steps, from real-time signals to smart data and then KPI. It enables a decrease of data volume while increasing their significance. The approach is tested through a decision-aid demonstrator that is applied to an industrial machining database.

Decision-Aid Conceptual Framework for machining

By embedding sensors into a machine-tool and connecting them to intelligent engines, new data and knowledge can be generated and interpreted. This knowledge, combined with other know-how and expert rules, is necessary to understand the current status of the machining process and its future changes. Therefore, one significant challenge is how to exploit this digital content for decision making and operational management perspectives. The proposed knowledge-based decision-aid framework provides indicators to support the operational management in the machining industry. As shown in Figure 1, three software modules are interconnected to support the decision-making process. The first module performs signal processing on sensor measures and collects the in-process heterogeneous data. The next step is a knowledge-based data mining. It uses contextual clustering and unsupervised machine learning algorithms introduced in [START_REF] Godreau | Continuous improvement of HSM process by data mining[END_REF]. Finally, new decision-aid indicators are computed and delivered to various departments according to their needs in terms of process management and improvement. A knowledgeengineering process is used to support the extraction, the formalization and the integration of other manufacturing knowledge.

To be efficient, a decision-aid framework should send the right information at the right moment to decision-makers. To do so, it should be able to manage large quantities of data and ensure flexible circulation of this volume of data throughout the enterprise network. To solve this issue, expert knowledge can be exploited for data aggregation. Only small quantities of meaningful machining-process information are then used for the construction of decision-aid indicators according to the expert intention and the business context. A knowledge-based and multilevel approach is proposed in the next section.

Knowledge-based multi-level aggregation

Contrary to pure data-driven approaches, the manufacturing knowledge is used in the proposed approach at each level to design the monitoring criteria or aggregation operators (Fig 2). It leads to meaningful information and to more reliable analyses of the production process. Different kinds of knowledge are considered and combined according to the expert point of view, e.g. mechanical model of machining, workpiece properties, business rules, as well as empirical knowledge. In this way, the real-time in-process machining signals lead to smart machining data and then to key performance indicators.

Level 0: real-time signals

A monitoring and data collection system, called EmmaTools, has been developed for machining [START_REF] De Castelbajac | Monitoring of distributed defects on HSM spindle bearings[END_REF]. Four accelerometers are integrated into the spindle of the machine-tool (in each radial direction, at the front and rear bearings). The signals of inprocess machining vibration are measured at a sampling rate of 25 kHz, which is suitable due to the high cutting speeds used in industry. The cutting power is also measured.

Level 1: monitoring and collected data

A record of 4 vibration signals at 25 kHz during months or years of industrial productions would result in an unnecessarily and incommodiously large volume of data. In addition, at each instant, in the frequency domain, only a few vibration components are relevant. Consequently, signal processing in realtime is highly suitable. The idea is to detect the occurrence of unsuitable specific events (rather than a simple collection of statistical indicators). Dedicated monitoring criteria were designed based on mechanical model of machining, in order to detect tool breakage, spindle failure or chatter (unstable cuts resulting in a too poor surface quality). Some workpiece properties (e.g. material) are also considered as knowledge to tune signal processing parameters (e.g. bandwidth). Order tracking is performed on the vibration signals in the frequency domain. Tool breakage is detected by increased tool unbalance (at spindle frequency), chatter by asynchronous cutting vibration and spindle failure by faulty bearing-induced vibration. Vibration level VRMS, a piece of empirical knowledge, is also evaluated. The monitoring criteria are computed every 0.1s. Many machining operational context data (such as machine-tool axes position and speed, or tool and program names) are also collected by the device at a sampling rate of 10 Hz (from the CNC of the machine-tool, by field bus). These operational context data are useful for the following aggregation steps. In this way, collected data (monitoring + operational context data) are synchronous.

Level 2: smart data

On the shopfloor, several workpieces are produced daily on a machine-tool, using several cutting tools for different machining operations. In the aeronautic industry, a machining operation with a given tool can last for several hours, with a large volume of corresponding data. It was chosen to aggregate the monitoring data into few smart data, at the time scale corresponding to one usage period of a given cutting tool. One tool usage period can include several machining operations.

The smart data, related to a given tool usage period, result from an accurate data selection (notably through the operational context data) and from aggregation operators. The latter can be a simple statistic operator e.g. the average power or feedrate during cuts, or the maximum spindle temperature. But these basic operators are not relevant for machining incidents or for the machining vibrations. In fact, machining incidents are very rare (e.g. one severe event for a few seconds per week) and need to be clearly emphasized. It is also known that a moderate vibration level is acceptable, contrary to severe vibration that can damage the tool or the spindle [START_REF]Machine tool spindle -Evaluation of machine tool spindle vibrations by measurements on spindle housing -Part 1: Spindles with rolling element bearings and integral drives operating at speeds[END_REF]. These are the reasons why a Criticality Operator (CO) has been proposed. It quantifies the excessive magnitudes, e.g. of vibrations, over a critical threshold Ti. Let Xi={xi(k), k=0…n} the time series representing an i th monitoring criterion. The operator is defined as: with ti and tf the initial and final time of the considered period (based on operational context data) and dt = 0.1s the sampling period. A second proposed operator T evaluates the time during which a criterion exceeds a critical threshold. Note that, in the case of vibrations, the monitoring values obtained from each of the four accelerometers are aggregated by a quadratic mean, before the aforementioned aggregation.

The critical thresholds Ti are obtained by unsupervised Machine Learning from the machining database (since it is generally not possible to collect labelled training data in industry, such as the few exact instants where chatter occurred). They become new manufacturing explicit knowledge that can be reused.

The definition of smart data and the way to select and aggregate the data are determined based on business rules and some tacit knowledge of manufacturing experts (issued from their experience). For instance, the smart data for a given tool usage period can be the duration of chatter (T[Nh>TNh], where Nh is a chatter detection criterion), the critical forced vibration (CO[VRMS>TVrms]), or the average cutting power. This aggregation step is time-consuming, but it can be performed in concurrent operation time. By this means, meaningful and easy-to-use smart data are obtained.

Level 3: decision-aid indicators

The last step is the aggregation of smart data (associated to a given tool-usage period), into Key Performance Indicators (KPI) for decision aid, dedicated to answer the specific managerial needs of different departments of a machining company. A decision-aid indicator is defined as a contextual instantiation of a set of KPI and smart data.

As for smart data, the definition of KPI and the way to select and aggregate smart data are determined based on manufacturing knowledge such as business rules and best practices. For example, the definition can be the sum or mean (weighted or not) of contextually selected smart data. It can also take the form of a comparison with the usual value of a given performance indicator that has been previously assessed by machine learning. Examples of KPI are provided in the use case (next section). Thereby, models defining KPI and smart data for machining constitute new formal knowledge.

The instantiation of KPI models rely on accurate data selection, which is enabled by operational context data. The definition of a decision-aid indicator is obtained according to four instantiation parameters, which describe the context of the decision:

 the objective of the decision for which decision-aid indicators are defined;  the decider who exploits the resulting decision-aid indicator according to a given business point of view;  the scope of analysis in terms of machine, workpiece, program, tool, or other entities;  the instantiation mode that indicates when a given decisionaid indicator is computed. Three modes are considered: periodic, on demand, and on event.

Industrial use case

For validation purpose, the decision-aid system has been firstly tested with 2 well-known industrial use cases. The present section details the data exploitation scenarios. The industrial use cases concern the high speed machining of structural aeronautic parts in aluminum alloy. These high added-value parts present thin walls and floors, subject to vibration issues (the main cause of non-quality, requiring a costly manual polishing on the part). The production system is also very flexible: a machine-tool almost never produces the same part twice consecutively. Consequently, a global approach is suitable for decision aid.

Application scenario: reporting for Manufacturing Dpt

A reporting scenario was conducted as an example of a decision-aid strategy. The reporting strategy is based on the instantiation parameters described in Section 3.4. The role of the reports is to provide useful knowledge to the manufacturing department, which enables a better understanding of the production process. Two examples were chosen for the demonstration in the aim to solve the chatter problems (unstable cuts that result in unacceptable surface roughness on the workpiece). The objective of the first one is to identify the worst cutting tools (for which the cutting conditions should be optimized). The second one highlights the main faulty workpiece programs (requiring updates). The reporting strategy follows a periodic mode, where decision-aid indicators are evaluated every 3 months.

To support the demonstration, an EmmaTools monitoring device was installed on a machine-tool of an aeronautic company. 39 Gb database of monitoring and operational context data were collected in-process, over 426 days of industrial production, where 80 cutting tools and 346 programs machined 534 workpieces. The volume of the database is reasonable (approximately 100 Mb of machining data per machine-tool and per day), thanks to the real-time processing of accelerometer signals that computes monitoring criteria, such as Nh for chatter detection. The second step aggregates the monitoring data into smart data, based on the chatter duration rule: T[Nh>20m/s²] for each tool usage period. It leads to a total volume of smart data of approximately 2.4 Mb per year and per machine-tool. The third aggregation step results in a decision-aid indicator as a combination of the KPIs for tool and program. The KPIs are evaluated through the total duration of chatter per tool and per program respectively, over 3 months of production.

Demonstrator results

To implement this scenario, a software demonstrator was developed for decision aid, based on the in-process machining data, and applied to the database collected in industry. Multi-Agent technology is used for information management between software modules composing the proposed decision-aid framework [START_REF] Kravari | A Survey of Agent Platforms[END_REF]. An agent is usually defined as a virtual autonomous entity that perceives the environment and acts upon an environment. The originality of the proposed approach is to adapt the concept of multi-agent to support interoperability and software interactions management. In fact, the communications are only between agents or between an agent and a software module, while the different software modules have no connection between them.

In the reporting scenario, four agents are used to monitor the execution of software components. The first agent is a Human-Machine Interface (HMI) that enables interaction with users. The Traceability agent aims to manage the connection to various data and knowledge bases. The third category is called a Computing agent and supports the remote execution of various Matlab programs in charge of data aggregation and performance indicator evaluation. The last category is the Reporting agent which is in charge of the definition and sharing of reports based on the four instantiation parameters of the reporting strategy. The heterogeneous agents interact in order to carry out the reporting process. Firstly, the decider is connected to the HMI agent for the configuration of the traceability and reporting parameters. Based on the request, the Traceability agent identifies in the repositories all data and knowledge required for the construction of decision-aid indicators. When the monitoring data are extracted from the database, the Computing agents manage the smart data aggregation and the KPI assessment. The Traceability agent stores the results. Finally, the Reporting agent extracts the computed indicators and generates a PDF report that presents the results in the appropriate format depending on the profile of the decider.

Figure 3 shows the demonstrator interfaces of the decision-aid process for the example of cutting tool scenario (defined in Section 4.1). The accelerometer signals (Fig. 3a) are processed in real time, obtaining monitoring criteria that are collected (e.g. VRMS, Fig. 3b). The smart data are computed in concurrent operation time. The HMI (Fig. 3c) enables the configuration of the reporting parameters, notably by selecting the traceability objective, the scope (program, tool, etc.) and the instantiation mode. The computations of KPIs are then carried out. The real time execution of the different agents is shown in Fig. 3d and3e. It results in the automatic sending of an e-mail to the decider with the appropriate PDF report (Fig. 3f). The report contains a visualization of the decision-aid indicator, here a pie chart of the chatter statistics per tool (Fig. 3g). It indicates that the tool reference 10 026 causes most of the chatter occurrences and needs an optimization of its cutting conditions. Fig. 3h provides more details about the summed chatter durations per tool that are visualized in the pie chart.

Since two well-known application scenarios had been chosen, the results were compared with the experts' expectations and have shown consistency. In this way, the quality of the system was validated.

Conclusion

This paper introduces an original knowledge-based and multilevel approach for data management and mining in machining. The different data aggregation steps are conducted with the integration of manufacturing knowledge, in order to pass from real-time signals to smart data, KPI and lastly to contextual decision-aid indicators. This strategy progressively reduces the volume of data and increases their meaningfulness.

The approach was implemented in a demonstrator of the decision-aid system. The architecture is based on a multi-agent system. An example of reporting scenario for the manufacturing department was successfully implemented. The related KPI were proposed and evaluated on a real machining database from the aeronautic industry. It led to the automatic sending of an e-mail, with the resulting decision-aid indicator, to the decider. The different levels of knowledge-based aggregated data are illustrated in the example. The efficiency of the proposed approach to tackle the big data issue was shown. Furthermore, the positive feedback from the industrial partners of the project provided a first validation.

More complex decision-aid scenarios are under development, concerning the analysis and the diagnosis of critical events and the prediction of failure.

Figure 1 .

 1 Figure 1. Decision-aid framework for machining.

Figure 2 .

 2 Figure 2. Knowledge-based multi-level data aggregation.

Figure 3 .

 3 Figure 3. Demonstrator results for the reporting scenario.

Acknowledgement

The authors acknowledge the financial support of the French government of the ANR SmartEmma project (ANR-16-CE10-0005) and they would also like to thank the contributions made by the industrial partners.


Andrés F López

Lopera St

John Nicolas

Durrande Mines Saint-Étienne

Gaussian Process Modulated Cox Processes under Linear Inequality Constraints

Gaussian process (GP) modulated Cox processes are widely used to model point patterns. Existing approaches require a mapping (link function) between the unconstrained GP and the positive intensity function. This commonly yields solutions that do not have a closed form or that are restricted to specific covariance functions. We introduce a novel finite approximation of GP-modulated Cox processes where positiveness conditions can be imposed directly on the GP, with no restrictions on the covariance function. Our approach can also ensure other types of inequality constraints (e.g. monotonicity, convexity), resulting in more versatile models that can be used for other classes of point processes (e.g. renewal processes). We demonstrate on both synthetic and real-world data that our framework accurately infers the intensity functions. Where monotonicity is a feature of the process, our ability to include this in the inference improves results.

INTRODUCTION

Point processes are used in a variety of real-world problems for modelling temporal or spatiotemporal point patterns in fields such as astronomy, geography, and ecology [START_REF] Baddeley | Spatial Point Patterns: Methodology and Applications with R[END_REF][START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF]. In reliability analysis, they are used as renewal processes to model the lifetime of items or failure (hazard) rates [START_REF] Cha | Point Processes for Reliability Analysis: Shocks and Repairable Systems[END_REF].

Proceedings of the 22 nd International Conference on Artificial Intelligence and Statistics (AISTATS) 2019, Naha, Okinawa, Japan. PMLR: Volume 89. Copyright 2019 by the author(s).

Poisson processes are the foundation for modelling point patterns [START_REF] Kingman | Poisson Processes. Oxford Studies in Probability[END_REF]. Their extension to stochastic intensity functions, known as doubly stochastic Poisson processes or Cox processes [START_REF] Cox | Some statistical methods connected with series of events[END_REF], enables non-parametric inference on the intensity function and allows expressing uncertainties [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF]. Moreover, previous studies have shown that other classes of point processes may also be seen as Cox processes. For example, [START_REF] Yannaros | On Cox processes and Gamma renewal processes[END_REF] proved that Gamma renewal processes are Cox processes under non-increasing conditions. A similar analysis was made later for Weibull processes [START_REF] Yannaros | Weibull renewal processes[END_REF].

Gaussian processes (GPs) form a flexible prior over functions, and are widely used to model the intensity process Λ(•) [START_REF] Møller | Log Gaussian Cox processes[END_REF][START_REF] Adams | Tractable nonparametric Bayesian inference in Poisson processes with Gaussian process intensities[END_REF][START_REF] Teh | Gaussian process modulated renewal processes[END_REF][START_REF] Gunter | Efficient Bayesian nonpara-metric modelling of structured point processes[END_REF][START_REF] Lasko | Efficient inference of Gaussianprocess-modulated renewal processes with application to medical event data[END_REF][START_REF] Lloyd | Variational inference for Gaussian process modulated Poisson processes[END_REF][START_REF] Fernandez | Gaussian processes for survival analysis[END_REF][START_REF] Donner | Efficient Bayesian inference of sigmoidal Gaussian Cox processes[END_REF]. However, to ensure positive intensities, this commonly requires link functions between the intensity process and the GP g (•). Typical examples of mappings are Λ(x) = exp(g(x)) [START_REF] Møller | Log Gaussian Cox processes[END_REF][START_REF] Diggle | Spatial and spatio-temporal log-Gaussian Cox processes: Extending the geostatistical paradigm[END_REF][START_REF] Flaxman | Fast Kronecker inference in Gaussian processes with non-Gaussian likelihoods[END_REF] or Λ(x) = g(x) 2 [START_REF] Lloyd | Variational inference for Gaussian process modulated Poisson processes[END_REF][START_REF] Kozachenko | Simulation of Cox random processes[END_REF]. The exponential transformation has the drawback that there is no closed-form expression for some of the integrals required to compute the likelihood. Although the square inverse link function allows closed-form expressions for certain kernels, it leads to models exhibiting "nodal lines" with zero intensity due to the non-monotonicity of the transformation (see [START_REF] John | Large-scale Cox process inference using variational Fourier features[END_REF], for a discussion). Furthermore, current approaches to Cox process inference cannot be used in applications such as renewal processes that require both positivity and monotonicity constraints.

Here, we introduce a novel approximation of GPmodulated Cox processes that does not rely on a mapping to obtain the intensity. In our approach we impose the constraints (e.g. non-negativeness or monotonicity) directly on Λ(•) by sampling from a truncated Gaussian vector. This has the advantage that the likelihood can be computed in closed form. Moreover, our approach can ensure any type of linear inequality constraint everywhere, which allows modelling of a broader range of point processes. This paper is organised as follows. In Section 2, we briefly describe inhomogeneous Poisson processes and some of their extensions. In Sections 3 and 4, we introduce a finite representation of GP-modulated Cox processes and the corresponding Cox process inference under inequality constraints. In Section 5, we apply our framework to 1D and 2D inference examples under different inequality conditions. We also test its performance in reliability applications with hazard rates exhibiting monotonic behaviours. Finally, in Section 6, we summarise our results and outline potential future work.

POISSON POINT PROCESSES

A Poisson process X is a random countable subset of S ⊆ R d where points occur independently [START_REF] Baddeley | Case Studies in Spatial Point Process Modeling[END_REF]. Let N ∈ N be a random variable (r.v.) denoting the number of points in X. Let X 1 , • • • , X n be a set of n independent and identically distributed (i.i.d.) r.v.'s on S. The likelihood of (N = n, X 1 = x 1 , • • • , X n = x n ) under an inhomogeneous Poisson process with non-negative intensity λ(•) is given by [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF]

f (N,X1,•••,Xn) (n, x 1 , • • • , x n ) = exp(-µ) n! n i=1 λ(x i ), (1)
where

µ = S λ(s) ds (2)
is the intensity measure or overall intensity.

When S is the real line, the distance (inter-arrival time) between consecutive points of a Poisson process follows an exponential distribution. Renewal processes are a generalisation of Poisson processes where interarrival times are i.i.d. but not necessarily exponentially distributed. An example is the Weibull process where inter-arrival times are distributed following λ(x) = αβx β-1 [START_REF] Cha | Point Processes for Reliability Analysis: Shocks and Repairable Systems[END_REF].

Cox processes [START_REF] Cox | Some statistical methods connected with series of events[END_REF] are a natural extension of inhomogeneous Poisson processes where λ(•) is sampled from a non-negative stochastic process Λ(•).

Previous studies have shown that many classes of point processes can be seen as Cox processes under certain conditions [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF][START_REF] Yannaros | On Cox processes and Gamma renewal processes[END_REF][START_REF] Yannaros | Weibull renewal processes[END_REF]. For example, Weibull renewal processes are Cox processes for β ∈ (0, 1] [START_REF] Yannaros | Weibull renewal processes[END_REF]. This motivates the construction of GPs with non-negative and monotonic constraints, so that they can be used as intensities Λ(•) of Cox processes.

APPROXIMATION OF GP MODULATED COX PROCESSES

In this work, we approximate the intensity Λ(•) of the Cox process by a finite-dimensional GP Λ m (•) subject to some inequality constraints (e.g. boundedness, monotonicity, convexity). Since positiveness constraints are imposed directly on Λ m (•), a link function is no longer necessary. This has two main advantages. First, the likelihood (1) can be computed analytically. Second, as our approach ensures any linear inequality constraint, it can be used for modelling a broader range of point processes.

Finite Approximation of 1D GPs

Let Λ(•) be a zero-mean GP on R with arbitrary covariance function k. Consider x ∈ S, with compact space S = [0, 1], and a set of knots

t 1 , • • • , t m ∈ S.
Here we consider equispaced knots t j = (j -1)∆ m with ∆ m = 1/(m -1). We define Λ m (•) as the finitedimensional approximation of Λ(•) consisting of its piecewise-linear interpolation at knots

t 1 , • • • , t m , i.e., Λ m (x) = m j=1 φ j (x)ξ j , (3)
where ξ j := Λ(t j ) for j = 1, • • • , m, and φ 1 , • • • , φ m are hat basis functions given by

φ j (x) := 1 - x-tj ∆m if x-tj ∆m ≤ 1, 0 otherwise. ( 4
)
Similarly to spline-based approaches (e.g., [START_REF] Sleeper | Regression splines in the Cox model with application to covariate effects in liver disease[END_REF], we assume that Λ(•) is piecewise defined by (first-order) polynomials. The striking property of this basis is that satisfying the inequality constraints (e.g. boundedness, monotonicity, convexity) at the knots implies that the constraints are satisfied everywhere in the input space [START_REF] Maatouk | Gaussian process emulators for computer experiments with inequality constraints[END_REF]. Although it is tempting to generalise the above construction to smoother basis functions, it makes this property difficult to enforce.

We aim at computing the distribution of Λ m (•) under the condition that it belongs to a convex set of functions E defined by some inequality constraints (e.g. positivity). This piecewise-linear representation has the benefit that satisfying Λ m (•) ∈ E is equivalent to satisfying only a finite number of inequality constraints. More precisely,

Λ m (•) ∈ E ⇔ ξ ∈ C, (5)
where

ξ = [ξ 1 , • • • , ξ m ] ,
and C is a convex set on R m .

For non-negativeness conditions E + , C is given by and for non-increasing conditions E ↓ , C is given by

C + := {c ∈ R m ; ∀ j = 1, • • • , m : c j ≥ 0}, (6)
C ↓ := {c ∈ R m ; ∀ j = 2, • • • , m : c j-1 ≥ c j }. (7)
Constraints can be composed, e.g. the convex set of non-negativeness and non-increasing conditions is given by The effect of different constraints on samples from the prior Λ m (•) can be seen in Figure 1. Here we set m = 100 and use a squared-exponential (SE) covariance function 1 with covariance parameters σ 2 = 1, = 0.2. The samples were generated via Hamiltonian Monte Carlo (HMC) [START_REF] Pakman | Exact Hamiltonian Monte Carlo for truncated multivariate Gaussians[END_REF].

C ↓ + = C + ∩ C ↓ . Assuming that ξ is zero-mean Gaussian-distributed with covariance matrix Γ = (k(t i , t j )) 1≤i,j≤m ,

Application to 1D GP-Modulated Cox Processes

The key challenge in building GP-modulated Cox processes is the evaluation of the integral in the intensity measure. By considering Λ m (•) as the intensity of the Cox process, the intensity measure (2) becomes

µ m = 1 0 Λ m (x) dx = 1 0 m j=1 φ j (x)ξ j dx = m j=1 c j ξ j ,
where

c 1 = c m = ∆m 2 and c j = ∆ m for 1 < j < m. The likelihood of (N = n, X 1 = x 1 , • • • , X n = x n ) is f (N,X1,•••,Xn)|{ξ1,•••,ξm} (n, x 1 , • • • , x n ) = 1 n! exp - m j=1 c j ξ j n i=1 m j=1 φ j (x i )ξ j . (8) 1 SE covariance function: k(t, t ) = σ 2 exp(-(t-t ) 2 2 2 ).
Since (8) depends on r.v.'s ξ 1 , • • • , ξ m , it can be approximated using samples of ξ. To estimate the covariance parameters θ of the vector ξ, we can use stochastic global optimisation [START_REF] Jones | Efficient global optimization of expensive blackbox functions[END_REF].

Extension to Higher Dimensions

The approximation in (3) can be extended to grids in d dimensions by tensorisation. For ease of notation, we assume the same number of knots m and knot-spacing

∆ m in each dimension, but the generalisation to dif- ferent m 1 , • • • , m d or ∆ m1 , • • • , ∆ m d is straightforward. Consider x = (x 1 , • • • , x d ) ∈ [0, 1] d , and a set of knots per dimension (t 1 1 , • • • , t 1 m ), • • • , (t d 1 , • • • , t d m ). Then Λ m is given by Λ m (x) = m j1,•••,j d =1 i={1,•••,d} φ i ji (x i ) ξ j1,•••,j d , (9)
where ξ j1,•••,j d := Λ(t j1 , • • • , t j d ) and φ i ji are the hat basis functions defined in (4). Inequality constraints can be imposed as in [START_REF] López-Lopera | Finite-dimensional Gaussian approximation with linear inequality constraints[END_REF]. By substituting (9) in (2), we obtain

µ m = 1 0 Λ m (x) dx = m j1,•••,j d =1 i={1,•••,d} c ji ξ j1,•••,j d ,
with c ji defined as in 1D, and the likelihood is

f (N,X1,•••,Xn)|ξ (n, x 1 , • • • , x n ) = 1 n! exp - m j1,•••,j d =1 i={1,•••,d} c ji ξ j1,•••,j d (10) × n i=1 m j1,•••,j d =1 k={1,•••,d} φ ji (x i,k ) ξ j1,•••,j d , where x i,k is the k-th component of the point x i .
Due to the tensor structure of the finite representation, it becomes costly as the dimension d increases. The HMC sampler for truncated multivariate Gaussians from [START_REF] Pakman | Exact Hamiltonian Monte Carlo for truncated multivariate Gaussians[END_REF] follows the same dynamics as a classical HMC sampler, but the particle "bounces" on the boundaries if its trajectory reaches one of the inequality constraints. The computational complexity of each iteration scales linearly with the number of inequality conditions (e.g. m d for positiveness constraints) if the iteration does not require any reflection, but also increases with each bounce. Hence, in the best case, the computational complexity is O(m d ). However, this drawback could be mitigated by using sparse representations of the constraints [START_REF] Pakman | Exact Hamiltonian Monte Carlo for truncated multivariate Gaussians[END_REF], or using other types of designs of the knots (e.g. sparse designs).

COX PROCESS INFERENCE

Having introduced the model, we now establish an inference procedure for Λ(•) using the approximation Λ m (•). For readability, we only assume nonnegativeness constraints, i.e. ξ ≥ 0, but the extension to other types of constraints can be made by constructing a set of linear inequalities of the form l ≤ Aξ ≤ u, where A is a full-rank matrix encoding the linear operations, and l and u are the lower and upper bounds. In that case, results for Aξ|{l ≤ Aξ ≤ u} are similar as for ξ|{0 ≤ ξ < ∞}, and samples of ξ can be recovered from samples of Aξ, by solving a linear system.

Consider the non-negative Gaussian vector ξ and its sample χ. The posterior distribution of ξ conditioned on a point pattern (N = n, X

1 = x 1 , • • • , X n = x n ) is f ξ|{N =n,X1=x1,•••,Xn=xn} (χ) (11) ∝ f (N,X1,•••,Xn)|{ξ=χ} (n, x 1 , • • • , x n ) f ξ (χ),
where the likelihood is defined in (8) and f ξ (χ) is the (truncated) Gaussian density given by

f ξ (χ) = exp -1 2 χ Γ -1 χ ∞ 0 exp -1 2 s Γ -1 s ds , for χ ≥ 0. ( 12
)
Since the posterior distribution ( 11) can be approximated using samples of ξ, it is possible to infer Λ m (•) via Metropolis-Hastings.

Metropolis-Hastings Algorithm with Truncated Gaussian Proposals

The implementation of the Metropolis-Hastings algorithm requires a proposal distribution q for the next step in the Markov chain. In practice, Gaussian proposals are often used, leading to the famous randomwalk Metropolis algorithm [START_REF] Murphy | Machine Learning: A Probabilistic Perspective (Adaptive Computation And Machine Learning[END_REF]. However, since inequality constraints are not necessarily satisfied using (non-truncated) Gaussian proposals, the standard random walk can suffer from small acceptance rates due to constraint violations. We propose as an alternative a constrained version of the randomwalk Metropolis algorithm where inequality conditions are ensured when sampling from the proposal q. As ξ is (non-negative) truncated Gaussian-distributed (with covariance matrix Γ), we suggest the truncated Gaussian proposal q given by q

(χ k+1 |χ k ) (13) = exp -1 2 [χ k+1 -χ k ] Σ -1 [χ k+1 -χ k ] ∞ 0 exp -1 2 [s -χ k ] Σ -1 [s -χ k ] ds
, where χ k+1 , χ k ≥ 0 are samples of ξ and Σ is the covariance matrix. Sampling from q can then be performed via MCMC [START_REF] Pakman | Exact Hamiltonian Monte Carlo for truncated multivariate Gaussians[END_REF].

We use Σ = ηΓ, where η is a scale factor. This has the benefit that we are sampling from a distribution with similar structure to the true one, while η controls the step size of the Metropolis-Hastings procedure and can be manually tuned to obtain a trade-off between mixing speed and acceptance rate of the algorithm.

The acceptance probability is given by

α k = f ξ|{N =n,X1=x1,•••,Xn=xn} (χ k+1 ) f ξ|{N =n,X1=x1,•••,Xn=xn} (χ k ) × β k , ( 14
)
where β k = q(χ k |χ k+1 )/q(χ k+1 |χ k ), and 4) and ( 8). We now focus on the term β k . Since the truncated Gaussian density has the same functional form as the non-truncated one, apart from the differing support and normalising constants, this yields

f ξ|{N =n,X1=x1,•••,Xn=xn} (χ) (15) = exp - 1 2 χ Γ -1 χ -c χ n i=1 φ (x i )χ is the (unnormalised) posterior distribution. φ(•) = [φ 1 (•), • • • , φ m (•)] and c = [c 1 , • • • , c m ] are defined in (
β k = ∞ 0 exp -1 2 [s -χ k ] Σ -1 [s -χ k ] ds ∞ 0 exp -1 2 [s -χ k+1 ] Σ -1 [s -χ k+1 ] ds . ( 16
)
The orthants

∞ 0 exp -1 2 [x -µ] Σ -1 [x -µ]
dx cannot be computed in closed form, but they can be estimated via MC [START_REF] Genz | Numerical computation of multivariate normal probabilities[END_REF][START_REF] Botev | The normal law under linear restrictions: Simulation and estimation via minimax tilting[END_REF]. Algorithm 1 summarises the implementation of the Metropolis-Hastings algorithm for the Cox process inference using the finite approximation of Section 3.

Algorithm 1 Metropolis-Hastings algorithm for Cox process inference with truncated Gaussian proposals k) , ηΓ such that χ ∈ C + .

1: Input: χ (0) ∈ (R m ) + , Γ (covariance matrix of ξ), η (scale factor). 2: for k = 0, 1, 2, • • • do 3: Sample χ ∼ N χ (

4:

Compute α k as in ( 14).

5:

Sample u k ∼ uniform(0, 1).

6:

Set new sample to 7:

χ (k+1) = χ , if α k ≥ u k χ (k) , if α k < u k . 8: Compute λ (k) m (x) = m j=1 φ j (x)χ (k) j
at location x with φ j (•) defined in (4).

Inference with Multiple Observations

For N o independent observations (X ν,1 , • • • , X ν,nν ) with ν = 1, • • • , N o , the acceptance probability follows α k = No ν=1 f ξ|{Nν =nν ,•••,Xν,n ν =xν,n ν } (χ k+1 ) No ν=1 f ξ|{Nν =nν ,•••,Xν,n ν =xν,n ν } (χ k ) β k , (17)
with posterior f ξ|{Nν =nν ,Xν,1=xν,1,•••,Xν,n ν =xν,n ν } and β k given by ( 11) and ( 16). Then, Algorithm 1 can be used with (17).

EMPIRICAL RESULTS

We test the performance of the finite approximation of GP-modulated Cox process on 1D and 2D applications. In the following, we use the squared-exponential covariance for the Gaussian vector ξ so that we can compare to [START_REF] Lloyd | Variational inference for Gaussian process modulated Poisson processes[END_REF]. We estimate the covariance parameters θ = (σ2 , ) by maximising the likelihood (8). For all numerical experiments, we fix m such that we obtain accurate resolutions of the finite representations while minimising the cost of MCMC (see [START_REF] Bay | Generalization of the Kimeldorf-Wahba correspondence for constrained interpolation[END_REF][START_REF] Maatouk | Gaussian process emulators for computer experiments with inequality constraints[END_REF], for discussion about the convergence of the finitedimensional approximation of GPs). 2 For simulating ξ, we use the exact HMC sampler proposed by [START_REF] Pakman | Exact Hamiltonian Monte Carlo for truncated multivariate Gaussians[END_REF]. To approximate the Gaussian orthant probabilities from (16), we use the estimator proposed by Botev (2017) using 200 MC samples. We run Algorithm 1 with a scale factor η between 10 -3 and 10 -4 for a good trade-off between the mixing speed and the acceptance rate for each experiment. 3 The number of discarded burn-in samples until the Markov chains became stationary varied between 10 3 and 10 4 samples. The code was implemented in the R programming language based on the package lineqGPR (López-Lopera, 2018).

Examples with Multiple Observations

Here, we test our approach using the three toy examples proposed by [START_REF] Adams | Tractable nonparametric Bayesian inference in Poisson processes with Gaussian process intensities[END_REF],

λ 1 (x) = 2 exp{-x/15} + exp{-[(x -25)/10] 2 },
λ 2 (x) = 5 sin(x 2 ) + 6, λ 3 (x) = piecewise linear through (0, 2), (25, 3), (50, 1), (75, 2.5) and (100, 3).

The domains for λ 1 , λ 2 and λ 3 are S 1 = [0, 50], S 2 = [0, 5] and S 3 = [0, 100], respectively.

Figure 2 shows the inference results using N o = 1, 10, 100 observations sampled from the ground truth. With increasing number of observations the inferred intensity converges to the ground truth. Here, we fixed m = 100 and η = 10 -3 .

In Table 1, we assess the performance of our approach under non-negativeness constraints (cGP-C + ). We compare our inference results to the ones obtained with a log-Gaussian process (log-GP) modulated Cox process [START_REF] Møller | Log Gaussian Cox processes[END_REF] and Variational Bayes for Point Processes (VBPP) [START_REF] Lloyd | Variational inference for Gaussian process modulated Poisson processes[END_REF] using the Q 2 criterion. This criterion is defined as

Q 2 = 1 -SMSE(λ(•), λ(•))
, where SMSE is the standardised mean squared error [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning[END_REF]. Q 2 is equal to one if the inferred λ(•) is exactly equal to the true λ(•), zero if it is equal to the average intensity λ, and negative if it performs worse than λ. We compute the Q 2 indicator on a regular grid of 1000 locations in S. Then, we compute the mean µ and one standard deviation σ of the Q 2 results across 20 different replicates. Table 1 shows that our approach outperforms its competitors, with consistently higher means of the Q 2 results and lesser dispersion σ.

We assess the computational cost of our approach using the third toy example λ 3 for N o = 100 (which has the largest number of events with on average 22500 events in total). Obtaining one sample using our approach takes around 60 milliseconds, and generating all 10 4 samples takes 10 minutes in total (in contrast to the 18 minutes required by VBPP). 4 The multivariate effective sample size (ESS) [START_REF] Flegal | mcmcse: Monte Carlo standard errors for MCMC[END_REF] was estimated at 322, corresponding to an effective sampling rate of 0.536 s -1 .

Modelling Hazard Rates in Renewal Processes

Poisson processes have been extended to model renewal processes where intensity functions are seen as hazard rates defining the probability that an operating object fails [START_REF] Serfozo | Renewal and Regenerative Processes[END_REF][START_REF] Cha | Point Processes for Reliability Analysis: Shocks and Repairable Systems[END_REF]. However, in many application, e.g. reliability engineering and survival analysis, hazard rates exhibit monotonic behaviours describing the degradation of items or lifetime of organisms. For example, the hazard functions for the failure of many mechanistic devices and the mortality of adult humans tend to exhibit monotonic behaviours. Thus, taking monotonicity constraints into account in renewal processes is crucial for the study of many applications. Moreover, it is known that introducing monotonicity information in GPs can lead to more realistic uncertainties [START_REF] Riihimäki | Gaussian processes with monotonicity information[END_REF][START_REF] Maatouk | Gaussian process emulators for computer experiments with inequality constraints[END_REF].

As discussed in Section 2, some renewal processes can be seen as Cox processes under certain conditions. In order to demonstrate that we can model other types of point patterns, here we use two toy examples where hazard rates are known to be monotonic. Both examples are inspired by two classical renewal process: Weibull process and Gamma process.

For the first class, the Weibull hazard function is

λ W (x) = αβx β-1 for x ≥ 0, (18)
where α and β are the scale and shape parameters, respectively. Depending on β, λ W can be either nonincreasing (0 < β < 1), constant (β = 1), or nondecreasing (β > 1). Moreover, for β ∈ (0, 1], the Weibull renewal process can be seen as a Cox process [START_REF] Yannaros | On Cox processes and Gamma renewal processes[END_REF]. For numerical experiments, we consider the case of non-increasing conditions in the domain S = [0, 100] by fixing α = 1 and β = 0.7 (see Figure 3). We test our framework using N o = 100 observations from λ W , and we consider non-negativeness conditions, with (cGP-C ↓ + ) or without (cGP-C + ) taking into account the non-increasing constraint. We also consider the case where λ W is non-increasing and convex (cGP-C ↓ + ). For the Gamma class, the hazard function is given by where Γ(•) and Γ x (•) are the Gamma function and the incomplete Gamma function, respectively [START_REF] Cha | Point Processes for Reliability Analysis: Shocks and Repairable Systems[END_REF], and α and β are the scale and shape parameters. As for the Weibull process, different behaviours can be obtained using different values of β.

λ G (x) = α x β-1 e -x Γ(β) -Γ x (β) , for x ≥ 0, (19)
Since similar profiles are obtained for β ∈ (0, 1], here we are interested in the case where λ G exhibits nondecreasing constraints (β > 1). We fix S = [0, 5], α = 5 and β = 1.7 obtaining a non-decreasing profile as shown in Figure 3. Here, we consider non-decreasing (cGP-C ↑ + ), and non-decreasing and concave (cGP-C ↑ + ) constraints. Since λ G (x) < α for x ∈ S, we add the constraint λ G ∈ [0, α].

Figure 3 shows the inferred intensities of λ W and λ G under the different conditions previously discussed. In both experiments, we fixed m = 100 and η = 10 -4 . For the Weibull class λ W , the performance of all three models, cGP-C + , cGP-C ↓ + and cGP-C ↓ + , tends to be similar. However, the model without monotonicity constraint exhibits undesired oscillations, whereas the other two approaches provide more realistic decreasing profiles and more accurate inference results for x > 50. We can also observe that the three models cannot learn the singularity at x = 0. Note that the proposed methodology does not make any assumption on the kernel, and it would be possible to consider a covariance function such as k(x, y)/(xy) in order to improve the model behaviour for small and large values of x. For the Gamma hazard function λ G , one may clearly observe the benefits of adding the non-decreasing and concave constraints, obtaining absolute improvements between 0.8% and 3.5% of the Q 2 indicator. Both examples of Figure 3 show that the monotonicity and convexity conditions found in certain point processes can be difficult to learn directly from the data. This suggests that including those constraints in the GP prior is necessary to get accurate models with more realistic uncertainties.

2D Redwoods Data

We now assess the performance of the proposed approach for a 2D spatial problem. We use the dataset provided by [START_REF] Ripley | Modelling spatial patterns[END_REF] which describes the locations of redwood trees. The dataset contains n = 195 events scaled to the unit square (see Figure 4). Here we choose m = 15, obtaining 225 knots in total, to obtain a good trade-off between resolution and computational cost. We use the product of two SE kernels with covariance parameters θ = (σ 2 , 1 , 2 ) as the covariance function of the Gaussian vector ξ, and we choose η = 10 -4 in Algorithm 1. Following the burn-in step, we keep 10 5 samples for the inference of λ(•), yielding a total running time of 7.6 hours (i.e. a sampling rate of approximately 4 s -1 ).

Figure 4 shows the normalised inference results for the redwood dataset for different values of the lengthscale parameters. Since in our approach we directly impose the inequality conditions on the Gaussian vector ξ instead of using a link function, the interpretation of the lengthscale parameters ( 1 , 2 ) are the same as for standard GPs: one can find a trade-off between fidelity and regularity by tuning . One can note, from Figures 4(a properly learn the point patterns but more regularity is exhibited when 1 = 2 = 10 -1 . For the case 1 = 2 = 10 -2 , although the model follows the point patterns, one may observe noisy behaviour in regions without points, e.g. around (x 1 , x 2 ) = (0.30, 0.85), as small values of lead to more oscillatory Gaussian random fields. Finally, we infer λ(•) when the covariance parameters θ are estimated via maximum likelihood using (10). According to the estimated lengthscales ( 1 = 0.055, 2 = 0.084), one can conclude that the estimated intensity λ(•) is smoother along the second dimension x 2 . This is in agreement with the inference results by [START_REF] Adams | Tractable nonparametric Bayesian inference in Poisson processes with Gaussian process intensities[END_REF], where more variations of λ(•) were exhibited across x 1 .

CONCLUSIONS

The proposed model for GP-modulated Cox processes is based on a finite-dimensional approximation of a GP that is constrained to be positive. This approach shows several advantages. First of all, it is based on general linear inequality constraints so it allows us to incorporate more information, such as monotonicity and convexity, in the prior. As seen in the experiments, this appears to be particularly helpful when few data are available. Second, imposing directly the positivity constraint on the GP makes the use of a link function unnecessary. Both the likelihood and the intensity measure can be computed analytically, which is not always the case when using a link function. Finally, the fact that our model is based on a finite-dimensional representation ensures that the computational burden grows linearly with the number of observations.

There are two key elements that make the method work: (a) the finite-dimensional representation of the GP that ensures that the constraints are satisfied everywhere, and (b) the dedicated MCMC proposal distribution based on a truncated normal distribution which allows us to have high acceptance rates compared to a naive multivariate Gaussian proposal.

The main limitation regarding the scaling of the proposed method lies in the dimension of the input space. This is due to the construction by tensorisation of the basis functions used to obtain the finite-dimensional representation. Moreover, our model is also sensitive to three parameters: the dimensionality of the space in which we perform HMC, the number of constraints, and the number of times the HMC particles violate a constraint. However, we believe that these limitations are not inherent to the proposed model and that other types of designs of the knots (e.g. sparse designs) could be used in high dimensions.

Figure 1 :

 1 Figure 1: Samples from the prior Λ m (•) under (a) no constraints, (b) non-negativeness constraints, (c) both non-negativeness and non-increasing constraints. The grey region shows the 95% confidence interval.

  then the distribution of ξ conditioned on ξ ∈ C is a truncated Gaussian distribution. Then, quantifying uncertainty on Λ m relies on sampling ξ ∈ C (see López-Lopera et al., 2018, for further discussion).

Figure 2 :

 2 Figure 2: Inference results with multiple observations (N o = 1, 10, 100) using the toy examples from Adamset al. (2009). Each panel shows the point patterns (black crosses), the true intensity λ (red dashed lines) the intensity inferred by the finite approximation of GP-modulated Cox processes (blue solid lines). The estimated 90% confidence intervals of the finite approximation are shown in grey.

Figure 3 :

 3 Figure 3: Renewal inference examples under different inequality constraints using N o = 100 and m = 100. Inference results are shown for (top row) a Weibull renewal process with α = 1 and β = 0.7, and (bottom row) a Gamma renewal process with α = 5 and β = 1.7. The panel description is the same as in Figure 2.

Figure 4 :

 4 Figure4shows the normalised inference results for the redwood dataset for different values of the lengthscale parameters. Since in our approach we directly impose the inequality conditions on the Gaussian vector ξ instead of using a link function, the interpretation of the lengthscale parameters ( 1 , 2 ) are the same as for standard GPs: one can find a trade-off between fidelity and regularity by tuning . One can note, from Figures4(a) and 4(b), that both profiles tend to

Table 1 :

 1 Q 2 results for the toy examples of Figure2, averaged over 20 ( † 10) replicates. Our results (cGP-

	Toy N o	log-GP	Q 2 (µ ± σ) [%] VBPP	cGP-C +
		1	51.2±30.1 51.9±26.1 65.7±14.3
	λ 1	10	95.1± 3.9	94.6± 3.7	95.4± 2.3
		100 99.5± 0.2 99.5± 0.3 99.5± 0.3
		1	-35.2±43.4	-1.1±28.8	0.7±24.0
	λ 2	10	72.6± 9.1	71.7±10.4 81.9± 7.4
		100	95.4± 0.7	92.1± 3.9	97.8± 0.6
		1	49.2±22.6 49.5±29.9 58.1±21.4
	λ 3	10 100	91.7± 4.4 98.4± 0.4 98.9±0.3 † 98.8± 0.3 93.8± 2.8 94.3± 2.5

C + ) are compared to results for Møller et al. (2001) (log-GP) and Lloyd et al. (2015) (VBPP).

We tested our model for various values of m, observing that, after a certain value, inference results are unchanged.

We observed convergence of Algorithm 1 for a wide range of values of η ∈ [10 -5 , 10 -2 ]. Fine-tuning of η can help the experiment run faster by gradually increasing η until the sampling mixes well.

These experiments were executed on a single core of an Intel R Core TM i7-6700HQ CPU.

Acknowledgements

This work was supported by the Chair in Applied Mathematics OQUAIDO (oquaido.emse.fr, France) and by PROWLER.io (www.prowler.io, UK). We thank O. Roustant (EMSE) and D. Rullière (ISFA) for their advice throughout this work.


B Augier

P Bot

F Hauville

EXPERIMENTAL VALIDATION OF UNSTEADY MODELS FOR WIND/ SAILS/ RIGGING FLUID SRUCTURE INTERACTION

The aim of this paper is to present the work of experimental validation elements of the aero elastic and unsteady model ARAVANTI. Numerical and Experimental results comparison is made on the rigging and sails of a J80 sail boat. Yacht modelling demands to consider unsteady phenomena resulting from the sea state, variations of wind speed and direction, yacht motion or trimming by the crew. A dedicated instrumentation is developed to measure the loads in shrouds and tension points of the sail, the apparent wind, the yacht motion, the sails flying shape and the navigation data. A special effort is made on sensors calibration, physical measurement comprehension and data synchronisation. Comparison with numerical results shows that the loads and flying shapes are well predicted by the model.

INTRODUCTION

A navigating sail boat is a dynamic system exposed to unsteady phenomena as the wind and the sea. Yacht response modeling demands to simulate the complex and changing environment. The present work is driven by the need to get unsteady measurements from an instrumented sail boat in real navigation conditions. Collected data are used to build a database for an aero elastic model validation. To measure loads in the rigging, ship attitudes, navigation data, adjustments and sails flying shape simultaneously a dedicated instrumentation has been developed. Figure 2 shows how full scale data and calculation code ARAVANTI are connected. ARAVANTI is a fluid-structure model using a CST elements membrane model extended in 3dimensions [START_REF] Roux | Strongly coupled VPP and CFD RANSE code for sailing yacht performance prediction' 3rd High Performance Yacht Design Conference[END_REF]. Hypothesis imposed inside this element are constant stresses, constant strains and uniform stiffness of the material. Non-linearities coming from the geometry and compressions are taken into account. The calculation of the flow around the sails is carried out under the hypothesis of an incompressible inviscid fluid, using a particular method developed by Rehbach (Rehbach (1978)) and then Huberson (Huberson (1984). This method is, in essence, unsteady, taking into account the boundary conditions of the displacement velocities of the surface as well as an atmospheric wind gradient. The effects of the interaction are translated into a coupling of the kinematic equation (continuity of the normal component of the velocity at the interface between fluid and structure geometrical domains) and dynamic equations (continuity of the normal component of the external force, pressure forces, on the contact surface of the sail with the fluid). ARAVANTI settings are given by yacht attitude recordings and calculation results are compared to the loads measured on navigation. Previous studies focused on specific sailing parameters like the aerodynamic or hydrodynamic forces with a dedicated sail force dynamometer boat DYNA or Fujin [START_REF] Masuyama | Database of sail shapes versus sail performance and validation of numerical calculations for the upwind condition[END_REF] [4], the flying shape and rig position [START_REF] Pelley | V-SPARS: A combined sail and rig recognition system using imaging techniques' 3 rd High Performance Yacht Design Conference Auckland[END_REF], or the aerodynamic sail force [START_REF] Hansen | Comparison of wind tunnel and full-scale aerodynamic sail force measurements' 2 nd High Performance Yacht Design Conference Auckland[END_REF]. The sailing boat problem is simplified to a small identified number of studied interactions, enough to be a tough work. Wind tunnel experimentation gives very accurate information on the apparent wind angle and apparent wind speed. AWA and trim are controlled. All those parameters are fixed during the set. They cannot represent the unsteady behavior of the wind, the sea state modeled by the heading, the trim or the AWA and the repercussion of a variation on the sail boat equilibrium. Sailors know that adjustments to go fast in a steady state are not the same than in waves and shifty wind ones.

EXPERIMENTAL APPARATUS

SENSORS SITING

Full scale testing is performed on a J80 which principal dimensions are presented in table 1. Sensors have to fit to the normal sailing configuration without disturbing the flow or the load. Figure 1 shows the position of all the sensors on the J80. 7 instrumented turnbuckles take the place of the 6 shrouds turnbuckles, V1, V2 and D1 starboard and portside, and the forestay chain plate. Instrumented shackles are disposed on all the sail load points. Five are mounted on the main sail: outhaul, sheet, halyard, cunningham and boom vang, three on the jib: sheet, halyard and tack. A ninth shackle is placed on the backstay. The 16 load sensors are linked to two dedicated analogical data acquisition and synchronization Spiders8 from HBM located inside the boat. Four analogical cameras are fixed on the top mast, two on the real top recording the main and two just under the forestay hound point recording the jib. Additional cameras are fixed on the roof, recording the crew and the sail foot. The Motion sensor Xsens MTi-G is placed on the point F, the rotation center of the hull for the small angle at 0° heel angle. An ultrasonic 3D anemometer is fixed on the top mast and a loch has been installed on the J80 hull. Wire displacement sensors are fixed between the main car and the boom and the jib car and the clew to measure the sheet length. A ruder angle sensor is fixed on the helm basis. A fluxgate compas and a GPS are used inside the boat. The use of classical force sensors as S sensors has been put apart because their oversize and their weight needed to modify the rigging and restrained to get all the measurement points. That is the reason why turnbuckles and shackles have been instrumented with stress gauges and substituted to the basic fittings in collaboration with HBM, a company specialized in measurement.

Instrumented shackles and turnbuckles are equipped by a load resistive cell giving traction information shown in figure 3 in the final stage of their development. The resistive cell is linked to the data acquisition hardware with a string. Instrumented turnbuckles (Fig. 3) are identical to the one used in navigation, Sparcraft turnbuckle 116mm for Ø5mm wire. The same adjustment for the initial load in the rigging can be done. A load gauge full bridge is stuck on two flat lugs symmetrically machine-cut on the turnbuckle shank linked to the chain plate. Turnbuckles, because of their thin shape, work in pure traction effort.

The maximum designed load is 10000N. Calibration with a HBM measurement standard load sensor presents no repeatability problem and the precision error determination is simple to get. Table 2 gives turnbuckles sensibilities and precision errors.

(c) Instrumented shackles

Important R&D has been done to calibrate instrumented shackles, from a Whishard Ø6mm, which behavior during first calibration was subject to many interrogations. Calibration procedure and shackle shape has been upgraded. The first step was to machine-cut a diabolo, as shown in figure 3, fitted to the pin in order to keep the tension on the symmetry axis. The D shape of the shackle is determined to be a cause of non linearity. The shackle doesn't work in pure traction and the measure is polluted by compression. In order to offset, both branches of the shackle have been equipped by a connected full gauge bridge to average the load in each side and delete the compression effect. The asymmetry of the pin, threaded on a side and linked by a pivot on the other one, was shown to be a source of non repeatability. The shackle was machine-cut to have two bores and the pin was replaced by a bolt with nut, clamp on branches. Calibration of the last upgraded shackles gives really good results and an absolute error inferior to 50N. Results presented in this paper were obtained with the non final version shackle which has an absolute error around 75N.

(d) Shackle calibration

The calibration process was set up to reflect the reality and to make sure that sensors were studied in the same state than in sailing condition. Shackles were pulled by three in serial between two Dynema splices in order to simulate the tension of a sheet. The maximum designed load is 5000N but first navigation recordings permit to recalibrate the measuring range. Shackles supported a basket loaded and unloaded progressively with plumb weight from 0 to 3500N on 7steps. Normal distribution is characterised by a probability density function similar to a bell curve centred on the mean value, by a third order moment, Skewness λ equal to 0 which represents the symmetry and a forth order moment, Kurtosis coefficient β equal to 3. Regarding the condition of calibration and the repartition of the results shown in figure 4, we assume that measurements follow a normal distribution and that errors determination probabilistic tools calculated to a Gaussian variable can be used. With this hypothesis, the precision error is the root of the sum of the squared errors:

E E E E f l

SAIL SHAPE MEASUREMENT

A set of three parallel stripes is applied to each sail located at heights of 20%, 40% and 70% of the luff length. Cameras on the top mast get pictures of full stripes and foot. Images processing from movies are used to measure flying shape parameters defined in figure 5.

(a) Cameras distortion

To correct the lens distortion, cameras are calibrated using a method based on the work of Zhang [START_REF] Zhang | A flexible new technique for camera calibration[END_REF]. The method has been inspired from the camera calibration toolbox for Matlab. Images are corrected in post processing with the distortion coefficient determined for each camera. A camera is fixed on a calibrated support and pictures of the drawn stripes are taken from different cameras positions. Table 3 compares the parameters calculated from the pictures with softwares ASA and ISIS to the values geometrically computed from the cylindrical shape. ASA gives only relative information based on the computation of the chord length from the determination of luff and leech points.

Table 3. Comparison of sail analysis softwares on calibrated drawn visualization stripes. θ is the entry and exit angle which are the same for this cylindrical geometry.

ASA gives really good results, whatever the camera position. ISIS is really camera position sensitive and results suffer when the angle between the optic axis of the camera and the normal to Rf increases. Stripes recognition is a long manual procedure for both software.

A calibration grid is placed on the floor, at the same level of the drawn sail foot, in order to calibrate the camera position determination procedure. The calibration grid is smaller than the real deck grid. The calibrated camera support gives 3D coordinates, compared to extrinsic camera coordinate calculated by a method inspired from calibration toolbox for Matlab and the work of Zhang [START_REF] Zhang | A flexible new technique for camera calibration[END_REF]. This method enables to locate precisely the cameras in space as shown on figure 7.

Figure 7. Comparison between the calibrated and calculated cameras 3D coordinates.

NAVIGATION DATA Table 4. Navigation sensors parameters and accuracy

The instrumented sail boat has the classical navigation equipment found in all cruiser. Some instruments have been updated from the basic fittings to control the acquisition. The GPS position is directly read by the software. The 3D acoustic anemometer WindMaster located on the top mast gives information on vertical flow with the z component. It gives an apparent wind measurement. The sea yacht attitude is recorded by the motion sensor XSENS MTI-G, a GPS aided Attitude and Heading Reference System. This 6DOF measurement unit is used on the auto-pilot on some IMOCA 60' boats. Table 4 sets up navigation sensors parameters and accuracy.

REAL TIME AND SYNCHRONISATION

The aim of the instrumentation is to record unsteady phenomena. Measurements have to be linked to the other in order to be able to "play" synchronized full movie of the navigations. Difficulties are in the important number of sensors with heteroclite physical measurements, in the different signals type and frequency. Recordings are shared on two different computers, one dealing with sail shape, the second dealing with the other sensors. Both are synchronized with a GPS clock. A LabView real-time acquisition and synchronization home made software is based on the principle of a sensor information frame for each system clock top, sized up by the highest acquisition frequency. Oversampling is done for the sensor with low frequency acquisition as the GPS or log.

The software is able to support 10Hz acquisition. Taking apart sheets, all adjustments are announced, voice recorded and dated. Sheet in or case out are short and material. Figure 8 shows wind and boat attitude data for a 15minutes run. Tacks and principal adjustments are noticed.

EXPERIMENTAL PROCEDURE

STEADY STATE

First calculation is made for a steady state. The steady state is 10sec in the 15min run, from time 230s to 240s. The helmsman was keeping the boat in a good port side close haul and the wind and sea condition were steady.

Figure 9 shows the steady wind and boat attitudes recorded data during this period. Figure 10 presents the recorded load on the rigging during this steady state. Loads are shared between leeward and windward shrouds, main and jib loads points. Mean values are calculated over this 10sec duration and the code computes a steady state with these average inputs. In the experimental run shown on figure 12, we study the influence of a material adjustment on the flying shape and rigging loads. During the run, the backstay is shortly eased from 350N to 0. The helmsman and the crew keep the boat straight in a portside close haul in order not to disturb the measurements by other happenings. A duration of 40sec is studied in order to calculate the state of the boat before the easing and the state after. Mean values are calculated for 20 sec before the easing and 20 sec after. The code calculates 2 steady states separated by the backstay adjustment with the average results data.

Figure 13 shows rigging loads of the instrumented J80 during a backstay easing. The adjustment is notified by a vertical black line and occurs at time 179s. Loads gap are clearly visible on the forestay, windward shrouds and jib.

RESULTS

STEADY STATE

ARAVANTI code models the fluid-structure interaction of the steady state by playing the mean data calculated from the 10sec period results, which has been determined to be the best. As presented in figure 2, the model inputs from the experimental data are the trim of the rig and the sail, AWA and AWS, boat attitude and motion. Sails design shape, layout, material, rig mechanical characteristics are fixed model inputs for all J80 fluidstructure modelling. Figure 13 presents the Numerical and Experimental comparison on the loads of the instrumented sail boat based on the mean values calculated from the steady state. Windward shrouds loads is very well evaluated with a relative error <8%. Backstay, main sheet and halyard and jib tack loads have been very well calculated with an error <5%. The code seems to not perfectly spread the effort in the jib because the forestay load is under evaluated and the jib halyard load is over. Nevertheless jib force balance is respected.

Outhaul load comparison gives bad results because the main foot extremity has a strap, which rounds the boom in order to support the vertical effort and make the outhaul effort in the foot axis. This strap and the resultant friction have not been modelled, so the measured value does not represent the calculated load. The code gives a graphical representation of the model in order to study the flying shape calculation and to make some numerical/experimental comparison with the pictures taken from the top mast cameras in navigation.

Figure 14 shows the superposition of the calculated flying shape of the J80 for the steady state and the picture shot during the 10sec run. Calculated sails shape matched very well the recorded one, stripes fitting nicely between each other. The experimental jib was modified and the top batten was changed, which explains the angle in the picture. Rig and boom calculated fit very well to the picture as shows the blue superposition.

BACKSTAY EASE

Figure 15. a) Comparison of the calculated flying shape before and after a 150mm backstay ease. b) Comparison of the recorded flying shape before and after the backstay ease. ARAVANTI code models the Fluid Structure Interaction of the 15sec period before and after the backstay ease. The calculation inputs are average values of the 2 periods. Because of the small difference of sailing condition before and after the ease, the comparison of the loads presents no particular interest. We will focus our interest on the flying shape. Figure 15 presents the comparison of the computed 150mm backstay ease and the comparison of the recorded backstay ease. The stripes, luff and leech have been underlined. Both cases, computed and recorded, start with the same sailing parameters. Both exhibit the expected modification of the flying shape: the twist decreases and the maximum camber increases when the backstay load decreases. Table 5 and 6 present the calculated stripes parameters from ARAVANTI calculation which confirm the tendencies. Here, the main shape changes due to the backstay ease is moderate because the initial load in the backstay was small, for there is only 10 knots breeze. It is remarkable that such a light trim variation is actually measured by the experimental system and well predicted by the model.

Table 5. Stripes parameters with backstay load Table 6. Stripes parameters after backstay ease 5. CONCLUSION A dedicated instrumentation system has been developed on a J80 sail boat to record data in real navigation conditions to be able to validate a fluid structure interaction model applied to the aero elastic problem of yacht sails (ARAVANTI). The system is made to measure the navigation data, the yacht attitudes and motion, the loads in the standing and running rigging and the sails flying shape. Special care was devoted to measurement calibration and uncertainties. Recorded environmental data and yacht attitudes are used as inputs to the model, as well as sails and rig mechanical characteristics. The sails flying shape and loads in the rig computed by the model are then compared to the ones recorded during the experiment. Two runs have been considered: a "steady state" in which the experimental records have been averaged over 10sec in constant conditions; a trim change where the results have been compared with and without tension in the backstay, other parameters remaining the same. In both cases, the experimental and numerical results match very well. The steady state case enables to focus on the loads in the rigging and gives a good comparison. The backstay ease gives a good example of the model capacity to model the modification of the flying shape due to trim change.

The instrumentation system is still under development to improve its capacity to measure strongly unsteady conditions (waves, wind and course changes...). The measurement campaign presented in this paper was made with home made data acquisition software which was limited to a low frequency and not optimal for data synchronisation. In order to improve the system (better synchronisation, higher picture frequency) a new without backstay load with backstay loads acquisition process based on the software RTmaps from INTEMPORA is under development and further experimental campaigns are planned.

Figure 1 .

 1 Figure 1. Plan of a J80 with principal sensors' position.

Figure 2 .

 2 Figure 2. Methodology of the numerical/experimental comparison data.

  standard deviation calculated for both cases.

Figure 3 .

 3 Figure 3. Gauge-fitted shackle and turnbuckles designed for the probe-fitted J80. Load cells are under the black water tight mastic.

  Figure 4. a) Example of the comparison between the measurements and the normal distribution. Coefficients λ and β confirm the Gaussian distribution. b) Linear regression, fidelity and hysteresis gap representation for a full calibration process (schematic enlarged view).

Figure 5 .

 5 Figure 5. Definition of the flying shape parameters 2.2 (b) Cameras calibrationStripes are designed to belong to a plane parallel to a reference plane Rf, at a known altitude. Two calibration grids are fixed on the deck to represent jib and main planes Rf. Images 2D-data generation are transformed in 3D-data in the world coordinate system with this hypothesis. Angles between the optic axis of the camera and the normal to Rf, illustrated in figure6is minimized during the installation for an obvious perspective issue but cannot be deleted. Using the same method, the residual angle and the exact position of the camera retinal

Figure 6 .

 6 Figure 6. Coordinate systems of the calibration grid fixed on the deck 3D transposition and flying shape parameters determination have been calibrated. The equivalent of a sail with 3 parallel stripes is drawn on the wall of a known diameter cylindrical room. Stripes geometry and angle are exactly calculated from the room dimensions.A camera is fixed on a calibrated support and pictures of the drawn stripes are taken from different cameras positions. Table3compares the parameters calculated from the pictures with softwares ASA and ISIS to the values geometrically computed from the cylindrical shape. ASA gives only relative information based on the computation of the chord length from the determination of luff and leech points.

Figure 8 .

 8 Figure 8. Wind and boat attitudes recorded data for a 15min run with the instrumented J80.

  a) In harbour presettings Rigging adjustment is identical to regatta, and Num/Expe comparison is made on shroud loads. Adjustments are based on the gap between the turnbuckle 2 threaded rods and the load is controlled. To be sure that the simulation is made with the same adjustment (rake, hogging) pictures of the boat are made, from a known point, on the port in calm wind condition. Pictures are superposed to the simulation to tune the real position. Mast and shrouds wires mechanic behavior is studied in laboratory. Parameters are loaded in the code.3.1 (b) In navigationLength of main and jib sheets are measured by a wire displacement sensor based on an incremental coder. The main car is kept on the centre line, and the trimmer only adjusts the main sheet, measured from the car to the hound point on the boom. The jib sheave is blocked in a position during all navigation and the length is measure from the car block to the clew.

Figure 9 .

 9 Figure 9. Wind and boat attitudes recorded data for a 10sec steady state run with the instrumented J80.

Figure 10 .

 10 Figure 10. Loads recorded data for a 10sec steady state run with the instrumented J80.

Figure 11 .Figure 12 .

 1112 Figure 11. Loads recorded data of the main sail with the precision error bars. Precision errors of all load sensors are known from the calibration presented in part 2.1. Figure 11 presents a zoom of the main sail load in the sheet and halyard. Errors represent the thickness of the curves. All load results are associated to their error bars.3.3 BACKSTAY EASING

Figure 13 .

 13 Figure 13. Numerical and Experimental loads comparison on a portside close haul steady state calculation

Table 1

 1

		. Principal dimension of a J80
	HULL
	Length over all (m)		8.50
	Length of water line (m)		7.50
	Maximum breadth (m)		2.49
	Draft (m)		1.50
	Disp (ton)		1.45
	SAILS
	I	luff length (m)	8.04	8.41
	P goose neck (m)	8.92	1.04
	J E (m)	2.70	3.23
	2.1	LOAD MEASUREMENT
	2.1 (a) Dedicated force sensors

Table 2 .

 2 Instrumented

	turnbuckles sensibilities and errors
	precision
	2.1 (b) Instrumented turnbuckles

ACKNOWLEDGEMENTS

This work has been done in the VOILENAV project supported by the French Naval Academy, Ecole Navale. The authors are grateful to K-epsilon for continuous collaboration about ARAVANTI model. Technical support from BSG development DeltaVoiles and Intempora is also acknowledged. Thanks also to Didier Munck for his valuable involvement in this project.

8.

AUTHORS BIOGRAPHY