Zachary Al Hoffman
email: zachary.hoffman@str.us

Steve Huntsman
email: steve.huntsman@str.us

Benchmarking an algorithm for expensive high-dimensional objectives on the bbob and bbob-largescale testbeds

Keywords:

We report benchmarks for a recently developed algorithm on the bbob and bbob-largescale benchmarking testbeds in COCO. This algorithm is designed for expensive high-dimensional multimodal objectives (such as arise in hyperparameter optimization or via simulations), and this regime introduces challenges for benchmarking. In particular, while the COCO experimental procedure yields evidence that this algorithm improves on the state of the art, the COCO framework also exhibits shortcomings for the very low evaluation budgets involved here. Consequently, we also report the results of ad hoc experiments that demonstrate a more obvious if less rigorous advantage for the algorithm over its nearest competition.

INTRODUCTION

The paper [START_REF] Huntsman | Parallel black-box optimization of expensive highdimensional multimodal functions via magnitude[END_REF] introduces an optimization algorithm designed to handle intrinsically high-dimensional (𝐷 ⪆ 40, and with all dimensions relevant) optimization problems for very expensive (≥ 1 minute/evaluation) multimodal objective functions. This regime includes problems of fundamental interest (e.g., optimizing the results of simulations) that push against fundamental limits of benchmarking tools: for example, high dimensionality and multimodality coupled with low function evaluation budgets preclude achieving global minima as a realistic goal. Meanwhile, while the algorithm in [START_REF] Huntsman | Parallel black-box optimization of expensive highdimensional multimodal functions via magnitude[END_REF] is carefully designed with such problems in mind, its runtime overhead introduces an additional challenge: benchmarking itself becomes very expensive.

The COmparing Continuous Optimizers (COCO) platform [START_REF] Hansen | COCO: A platform for comparing continuous optimizers in a black-box setting[END_REF] enables automated benchmarking and performance comparisons for optimization algorithms. In this paper, we benchmark the algorithm in [START_REF] Huntsman | Parallel black-box optimization of expensive highdimensional multimodal functions via magnitude[END_REF]-hereinafter referred to as EXPLO2 in the text and EE in the COCO plots-on the bbob [5, 10] and bbob-largescale [START_REF] Elhara | COCO: the large scale blackbox optimization benchmarking (BBOB-largescale) test suite[END_REF] test suites, as well as performing a more ad hoc evaluation to address gaps in COCO for low function evaluation budgets.

ALGORITHM PRESENTATION

If 𝑑 is the distance matrix of a finite subset of R 𝐷 and 𝑡 ≥ 0, then 𝑍 𝑗𝑘 = exp(-𝑡𝑑 𝑗𝑘 ) is a positive definite matrix [START_REF] Steinwart | Support vector machines[END_REF], hence invertible. With this in mind, the corresponding weighting 𝑤 is the solution of 𝑍𝑤 = 1, where 1 denotes a vector of all ones; the magnitude of 𝑍 is ∑ 𝑗 𝑤 𝑗 [START_REF] Leinster | Entropy and Diversity: the Axiomatic Approach[END_REF]. Magnitude is a very general scale-dependent notion of an effective number of points that encodes rich geometrical data [START_REF] Leinster | The magnitude of a metric space: from category theory to geometric measure theory[END_REF]. Meanwhile, the components of a weighting meaningfully encode a notion of effective size per point.

Moreover, 𝑍 is a radial basis function (RBF) interpolation matrix [START_REF] Martin D Buhmann | Radial basis functions: theory and implementations[END_REF] and the equation 𝑍𝑤 = 1 amounts to defining a weighting 𝑤 as the vector whose components are coefficients for interpolating the unit function. So if {𝑥 𝑗 } 𝑛 𝑗=1 are points in R 𝐷 with distance matrix 𝑑 and 𝑦 𝑗 ∶= 𝑓 (𝑥 𝑗 ) for 𝑓 ∶ R 𝐷 → R, then the RBF interpolation is 𝑓 (𝑥) ≈ 𝑦𝑍 -1 𝜁 (𝑥), where 𝜁 𝑘 (𝑥) ∶= exp(-𝑡⋃ 𝑥 -𝑥 𝑘 ⋃) and we treat 𝑦 and 𝜁 respectively as row and column vectors.

EXPLO2 constructs a surrogate on a box in R 𝐷 that convexly combines i) the differential magnitude due to a new point in relation to a dynamically curated subset of prior evaluation points, and ii) the RBF interpolation relative to those same points. These two functions respectively embody notions of exploration and exploitation, and over the course of the algorithm the surrogate shifts from the former to the latter. Successive evaluation points are determined by optimizing the surrogate with an inner large-scale algorithm (default: fmincon).

EXPERIMENTAL PROCEDURE

We ran EXPLO2 with default settings as described in [START_REF] Huntsman | Parallel black-box optimization of expensive highdimensional multimodal functions via magnitude[END_REF] and with 𝑛 ∥ = 32 parallel workers on the entire bbob suite for 25 × 𝐷 evaluations according to [START_REF] Hansen | COCO: The Experimental Procedure[END_REF]. We also ran EXPLO2 in the same way on the multimodal functions 𝑓 15 , . . . , 𝑓 19 in the bbob-largescale suite for 2 × 𝐷 evaluations. 1 The limited budget and function selection were due to runtime overhead, as discussed in the sequel.

As [START_REF] Huntsman | Parallel black-box optimization of expensive highdimensional multimodal functions via magnitude[END_REF] details, a representative set of algorithms to benchmark against for low-dimensional, expensive, multimodal functions is NEWUOA, MCS, GLOBAL, SMAC-BBOB, and *-CMA-ES with * ∈ {lmm, DTS, lq}. Of these, only NEWUOA is well-suited to problems with hundreds of dimensions. Indeed, most of these algorithms have not even been benchmarked for dimension 40 in COCO as of this writing, and even NEWUOA eventually scales poorly with dimension [START_REF] Ros | Benchmarking the NEWUOA on the BBOB-2009 function testbed[END_REF]. Thus it is also important to perform large-scale benchmarking of EXPLO2, and to compare its performance to fmincon, which EXPLO2 is a wrapper for and is itself a large-scale algorithm.

This experimental procedure was useful, but still insufficient to adequately characterize the regime for which EXPLO2 was designed and where it exhibits advantage over competing algorithms, viz. very expensive high-dimensional multimodal objectives. In particular, this procedure did not yield plots in COCO that readily conveyed meaningful information (e.g., an "expensive plot" option does not appear to be available for the bbob-largescale suite). We therefore also performed a more ad hoc comparison of EXPLO2 to VD-CMA-ES [START_REF] Akimoto | Comparison-based natural gradient optimization in high dimension[END_REF], the best performing large-scale algorithm for structured multimodal functions in [START_REF] Varelas | Benchmarking large scale variants of cma-es and l-bfgs-b on the bbob-largescale testbed[END_REF]. As Figures 10111213show, EXPLO2 performs better with respect to evaluation counts, regardless of parallelism. However, it also incurs a very large runtime overhead.

CPU TIMING

In fact, the CPU timing of EXPLO2 is orders of magnitude higher than competing algorithms, as it wraps a large amount of dense linear algebra (albeit with carefully controlled runtime, and also with prospects for eventually achieving sparsity) around an inner optimizer. Therefore, EXPLO2 is only indicated instead of *-CMA-ES for objectives whose evaluation requires on the order of minutes or more. For example, the bbob data presented here took roughly a week to generate; the bbob-largescale data was similarly very time-consuming. For this reason, precise timing measurements were largely pointless, and high function evaluation budgets were impractical.

RESULTS

Results from experiments according to [START_REF] Hansen | COCO: The Experimental Procedure[END_REF] and [START_REF] Hansen | COCO: Performance Assessment[END_REF] on the benchmark functions given in [START_REF] Elhara | COCO: the large scale blackbox optimization benchmarking (BBOB-largescale) test suite[END_REF]5] are presented in Figures 1,23, and 4-8; see also Table 1. 2 The experiments were performed and plots produced with COCO [START_REF] Hansen | COCO: A platform for comparing continuous optimizers in a black-box setting[END_REF], version 2.4. Additionally, our ad hoc experiments with VD-CMA-ES are presented in Figures 10111213.

The expected runtime (ERT), used in Table 1, depends on a given target precision, 𝐼 target = 𝑓 opt + Δ𝑓 , and is computed over all relevant trials as the number of function evaluations executed during each trial while the best function value did not reach 𝑓 t , summed over all trials and divided by the number of trials that actually reached 𝑓 t [START_REF] Hansen | COCO: The Experimental Procedure[END_REF][START_REF] Price | Differential evolution vs. the functions of the second ICEO[END_REF]. Statistical significance is tested with the rank-sum test for a given target Δ𝑓 t using, for each trial, either the number of needed function evaluations to reach 𝑓 t (inverted and multiplied by -1), or, if the target was not reached, the best Δ𝑓value achieved, measured only up to the smallest number of overall function evaluations for any unsuccessful trial under consideration.

Because of the runtime overhead of EXPLO2, our experiments used a fixed budget of 25 function evaluations per dimension (i.e., 25 × DIM function evaluations in total). 3 As [START_REF] Hansen | COCO: The Experimental Procedure[END_REF] points out, algorithms are only comparable up to the smallest budget given to any of them, corresponding in this case to log 10 25 ≈ 1.40 on the horizontal axis for Figures 2-3 and 4-7, and to log 10 2 ≈ 0.3 for Figure 8. This corresponds in our case exactly to the location of crosses (×), which indicate where bootstrapping of experimental data begins to estimate results for larger numbers of function evaluations. At the same time, we used 𝑛 ∥ = 32 parallel workers, 4and log 10 (25⇑32) ≈ -0.107. Thus allowing for parallel resources (see §5.1) suggests comparing EXPLO2 at the value log 10 25 on the horizontal axes with the other algorithms at log 10 (25⇑32), except for *-CMA-ES, which is parallelizable. 5 6 To summarize, as the dimension of structured multimodal problems grows, EXPLO2 outperforms NEWUOA, particularly when taking parallelism into account; meanwhile, the quadratic scaling with dimension of non-large-scale *-CMA-ES algorithms puts them at an increasing disadvantage since EXPLO2 has essentially no runtime dependence on dimension other than through its inner solver, which here is the default large-scale interior-point fmincon algorithm. Finally, on structured multimodal functions, EXPLO2 outperforms VD-CMA-ES, which is otherwise the best-performing large-scale algorithm that we are aware of.

Parallel performance

While EXPLO2 is parallelized (at marginal cost to performance), to the best of our knowledge no competing technique is except for *-CMA-ES. Indeed, as [START_REF] Gao | Distributed Gauss-Newton optimization method for history matching problems with multiple best matches[END_REF] points out, most surrogate-based derivativefree optimizers-including NEWUOA-require sequential function evaluations. Though parallel algorithms exist [START_REF] Raphael T Haftka | Parallel surrogate-assisted global optimization with expensive functions-a survey[END_REF][START_REF] Rehbach | Comparison of parallel surrogate-assisted optimization approaches[END_REF][START_REF] Xia | GOPS: efficient RBF surrogate global optimization algorithm with high dimensions and many parallel processors including application to multimodal water quality PDE model calibration[END_REF], these are relatively few in number outside the context of Bayesian optimization (which is unsuitable for high-dimensional problems), and parallel techniques appropriate for high-dimensional problems have been considered in our design and/or benchmarking.

With this in mind, since our experiments use 𝑛 ∥ = 32, exceeding our per-dimension evaluation budget of 25, it is obvious that EXPLO2 can outperform any of the competing algorithms considered here on the number of rounds of parallel function evaluations, except for *-CMA-ES, which becomes (depending on the variant employed) illsuited for use or less performant than EXPLO2 in high dimensions.

REMARKS

Some of the most important optimization problems in science and industry involve expensive objectives, e.g., the results of simulations. This fact and the experience of the present paper suggests introducing additional functionality in COCO to allow more precise benchmarking for extremely low function evaluation budgets. Experiments on EXPLO2 along similar lines to those presented here indicate that using generic CMA-ES as the inner optimizer incurs a collateral runtime cost orders of magnitude beyond that of fmincon even for 𝐷 = 2, which precludes careful benchmarking of that alternative.

Finally, though the runtime overhead of EXPLO2 scales favorably with dimension, it is still high for any given function evaluation: the surrogate is complicated and even a large-scale inner optimizer takes resources. EXPLO2 is therefore only suited for highdimensional functions that are very expensive to evaluate, and it will generally be advisable to evaluate the suitability of *-CMA-ES as well in any particular application.

15/15 lmm-C 4.9( This ERT ratio and, in braces as dispersion measure, the half difference between 10 and 90%-tile of bootstrapped run lengths appear for each algorithm and run-length based target, the corresponding reference ERT (preceded by the target Δ𝑓 -value in italics) in the first row. #succ is the number of trials that reached the target value of the last column. The median number of conducted function evaluations is additionally given in italics, if the target in the last column was never reached. Entries, succeeded by a star, are statistically significantly better (according to the rank-sum test) when compared to all other algorithms of the table, with 𝑝 = 0.05 or 𝑝 = 10 -𝑘 when the number 𝑘 following the star is larger than 1, with Bonferroni correction by the number of functions [START_REF] Varelas | A comparative study of large-scale variants of CMA-ES[END_REF]. A ↓ indicates the same tested against the best algorithm from BBOB 2009. Best results are printed in bold. [5] S. Finck, N.

.4e5 3/15 EE ∞ ∞ ∞ ∞ ∞ 501 0/15 GLOBA ∞ ∞ ∞ ∞ ∞ 5e4 0/15 MCS h 1.4(1) 1(0.7) 1(0.6) 1(2) 1(2) 3/15 NEWUO ∞ ∞ ∞ ∞ ∞ 2e6 0/15 lmm-C ∞ ∞ ∞ ∞ ∞ 8805 0/15 SMAC-∞ ∞ ∞ ∞ ∞ 2000 0/15 fminc 6.2e-4(2e-5) ⋆4 7.3e-4(7e-7) ⋆4 8.8e-4(7e-7) ⋆4 1.0e-3(9e-5) ⋆4 1.2e-3(9e-5) ⋆4

Figure 1 :

 1 Figure 1: Expected running time (ERT in number of 𝑓 -evaluations as log 10 value) divided by dimension versus dimension. The target function value is chosen such that the best algorithm from BBOB 2009 just failed to achieve an ERT of 10×DIM. Different symbols correspond to different algorithms given in the legend of 𝑓 1 and 𝑓 24 . Light symbols give the maximum number of function evaluations from the longest trial divided by dimension. Black stars indicate a statistically better result compared to all other algorithms with 𝑝 < 0.01 and Bonferroni correction number of dimensions (six). Legend: ○: EE mul25par32lamLin, ♢: GLOBAL pal noiseless, ⋆: MCS huyer noiseless, ▽: NEWUOA ros noiseless, : SMAC-BBOB hutter noiseless, △: fmincon pal noiseless, : lmm-CMA-ES auger noiseless.

Figure 2 :

 2 Figure 2: Bootstrapped empirical cumulative distribution of the number of objective function evaluations divided by dimension (FEvals/DIM) for all functions and subgroups in 10-D. The targets are chosen from 10 (-8..2⌋ such that the best algorithm from BBOB 2009 just not reached them within a given budget of 𝑘 × DIM, with 31 different values of 𝑘 chosen equidistant in logscale within the interval {0.5, . . . , 50}. As reference algorithm, the best algorithm from BBOB 2009 is shown as light thick line with diamond markers. multi-modal fcts weakly structured multi-modal fcts all fcts

Figure 3 :

 3 Figure 3: As in Figure 2, but for 20-D.

Figure 5 :Figure 6 :

 56 Figure 5: As in Figure 4, but for 𝑓 9 , . . . , 𝑓 16 .

Figure 9 :Figure 10 :Figure 11 :

 91011 Figure 9: As in Figure 8, but for dimensions 320 (top row) and 640 (bottom row).

Figure 12 :Figure 13 :

 1213 Figure 12: As in Figure 10, but for the Griewank-Rosenbrock function F8F2 (without rotations or shifts, but otherwise corresponding to 𝑓 19 in the BBOB testbed) on (-5, 5⌋ 𝐷 for 𝐷 = 20. Note that VD-CMA-ES often terminates before the budget is reached.

Table 1 :

 1 Expected runtime (ERT in number of function evaluations) divided by the respective best ERT measured during BBOB-2009 in dimension 20.

	15/15

  As in Figure4, but for dimension 40 and 𝑓 3 , 𝑓 4 , 𝑓 7 , 𝑓 15 , . . . , 𝑓 19 . As in Figure4, but for dimension 80 (top row) and 160 (bottom row) and 𝑓 16 , . . . , 𝑓 19 on the bbob-largescale testbed. The Rastrigin function 𝑓 15 has degenerate results for all algorithms and so is not shown. Note that[START_REF] Hansen | COCO: The Experimental Procedure[END_REF] points out that algorithms are only comparable up to the smallest budget given to any of them, corresponding in our case exactly to the location of crosses (×), which indicate where bootstrapping of experimental data begins to estimate results for larger numbers of function evaluations. These locations are very far to the left, though EXPLO2 (here indicated as EE) is consistently a top performer there.

	Fraction of function,target pairs	0.0 0.2 0.4 0.6 0.8 1.0	0 bbob f3, 40-D log10(# f-evals / dimension) 2 31 targets RLs/dim: 0.5..50 from refalgs/best2009-bbob.tar.gz 15 instances v2.4 3 Rastrigin separable	NEWUOA ro EE mul25p best 2009 Fraction of function,target pairs	0 bbob f4, 40-D log10(# f-evals / dimension) 2 31 targets RLs/dim: 0.5..50 from refalgs/best2009-bbob.tar.gz 15 instances v2.4 4 Skew Rastrigin-Bueche separ 0.0 0.2 0.4 0.6 0.8 1.0	EE mul25p NEWUOA ro best 2009 Fraction of function,target pairs	0.0 0.2 0.4 0.6 0.8 1.0	0 bbob f7, 40-D log10(# f-evals / dimension) 2 31 targets RLs/dim: 0.5..50 from refalgs/best2009-bbob.tar.gz 15 instances v2.4 7 Step-ellipsoid	EE mul25p NEWUOA ro best 2009 Fraction of function,target pairs	0.0 0.2 0.4 0.6 0.8 1.0	0 15 instances log10(# f-evals / dimension) 2 v2.4 from refalgs/best2009-bbob.tar.gz 31 targets RLs/dim: 0.5..50 bbob f15, 40-D 15 Rastrigin	NEWUOA ro EE mul25p best 2009
	Fraction of function,target pairs	0.0 0.2 0.4 0.6 0.8 1.0	0 bbob f16, 40-D log10(# f-evals / dimension) 2 31 targets RLs/dim: 0.5..50 from refalgs/best2009-bbob.tar.gz 15 instances v2.4 16 Weierstrass	EE mul25p NEWUOA ro best 2009 Fraction of function,target pairs	0.0 0.2 0.4 0.6 0.8 1.0	0 bbob f17, 40-D log10(# f-evals / dimension) 2 31 targets RLs/dim: 0.5..50 from refalgs/best2009-bbob.tar.gz 15 instances v2.4 17 Schaffer F7, condition 10	NEWUOA ro EE mul25p best 2009 Fraction of function,target pairs	0 bbob f18, 40-D log10(# f-evals / dimension) 2 31 targets RLs/dim: 0.5..50 from refalgs/best2009-bbob.tar.gz 30, 15 instances v2.4 18 Schaffer F7, condition 1000 0.0 0.2 0.4 0.6 0.8 1.0	NEWUOA ro EE mul25p best 2009 Fraction of function,target pairs	0 30, 15 instances log10(# f-evals / dimension) 2 v2.4 from refalgs/best2009-bbob.tar.gz 31 targets RLs/dim: 0.5..50 bbob f19, 40-D 19 Griewank-Rosenbrock F8F2 0.0 0.2 0.4 0.6 0.8 1.0	NEWUOA ro EE mul25p best 2009
	Fraction of function,target pairs	0.0 0.2 0.4 0.6 0.8 1.0	Figure 7: 0 2 4 6 log10(# f-evals / dimension) m2DLBFGS EE on bbo LBFGS Var LMCMA14 V R2ES Vare R10ES Var V2D-CMA V VkD-CMA V LMCMA17 V sepCMA Va VD-CMA Va CMA Varel bbob-largescale f16, 80-D 51 targets: 100..1e-08 1, 15 instances v2.4 16 Weierstrass 0 0.0 0.2 0.4 0.6 0.8 1.0 Fraction of function,target pairs bbob-largescale f17, 80-D 2 log10(# f-evals / dimension) 4 6 51 targets: 100..1e-08 1, 15 instances v2.4 17 Schaffer F7, condition 10	EE on bbo m2DLBFGS LBFGS Var LMCMA14 V R10ES Var R2ES Vare LMCMA17 V V2D-CMA V VkD-CMA V sepCMA Va CMA Varel VD-CMA Va Fraction of function,target pairs	0 bbob-largescale f18, 80-D 2 log10(# f-evals / dimension) 4 6 51 targets: 100..1e-08 1, 15 instances v2.4 18 Schaffer F7, condition 1000 0.0 0.2 0.4 0.6 0.8 1.0	EE on bbo LBFGS Var m2DLBFGS R10ES Var LMCMA14 V R2ES Vare LMCMA17 V V2D-CMA V VkD-CMA V sepCMA Va VD-CMA Va CMA Varel Fraction of function,target pairs	0 1, 15 instances 2 log10(# f-evals / dimension) 4 6 v2.4 51 targets: 100..1e-08 bbob-largescale f19, 80-D 19 Griewank-Rosenbrock F8F2 0.0 0.2 0.4 0.6 0.8 1.0	EE on bbo m2DLBFGS LBFGS Var LMCMA14 V R10ES Var R2ES Vare sepCMA Va VkD-CMA V V2D-CMA V LMCMA17 V CMA Varel VD-CMA Va
	Fraction of function,target pairs	0.0 0.2 0.4 0.6 0.8 1.0	0 bbob-largescale f16, 160-D 2 log10(# f-evals / dimension) 4 6 51 targets: 100..1e-08 1, 15 instances v2.4 16 Weierstrass	m2DLBFGS LBFGS Var EE on bbo LMCMA14 V R10ES Var R2ES Vare VkD-CMA V V2D-CMA V LMCMA17 V sepCMA Va CMA Varel VD-CMA Va Fraction of function,target pairs	0.0 0.2 0.4 0.6 0.8 1.0	0 bbob-largescale f17, 160-D 2 log10(# f-evals / dimension) 4 6 51 targets: 100..1e-08 1, 15 instances v2.4 17 Schaffer F7, condition 10	EE on bbo LBFGS Var m2DLBFGS LMCMA14 V R10ES Var R2ES Vare LMCMA17 V V2D-CMA V VkD-CMA V sepCMA Va VD-CMA Va CMA Varel Fraction of function,target pairs	0 bbob-largescale f18, 160-D 2 log10(# f-evals / dimension) 4 6 51 targets: 100..1e-08 1, 15 instances v2.4 18 Schaffer F7, condition 1000 0.0 0.2 0.4 0.6 0.8 1.0	EE on bbo LBFGS Var m2DLBFGS R2ES Vare LMCMA14 V R10ES Var LMCMA17 V V2D-CMA V VkD-CMA V sepCMA Va CMA Varel VD-CMA Va Fraction of function,target pairs	0 1, 15 instances 2 log10(# f-evals / dimension) 4 6 v2.4 51 targets: 100..1e-08 bbob-largescale f19, 160-D 19 Griewank-Rosenbrock F8F2 0.0 0.2 0.4 0.6 0.8 1.0	EE on bbo LBFGS Var m2DLBFGS LMCMA14 V R2ES Vare R10ES Var V2D-CMA V LMCMA17 V sepCMA Va VkD-CMA V VD-CMA Va CMA Varel
	Figure 8: #FEs/D	0.5	1.2	3		10	50	#succ	#FEs/D	0.5	1.2	3	10	50	#succ
	f21		6.3e+1:36	4.0e+1:77	4.0e+1:77		1.6e+1:456	4.0e+0:1094	15/15		f23	6.3e+0:29	4.0e+0:118	2.5e+0:306	2.5e+0:306	1.0e+0:1614	15/15
	EE	3.5(2)		2.9(3)	2.9(2)		1.3(1.0)	6.6(12)	1/15	EE	1.5(2)		11(14)	∞	∞	∞ 501	0/15
	GLOBA 9.2(0.7)		5.2(0.8)	5.2(0.8)		1.0(0.2)	1.00(0.9)	15/15	GLOBA 1.8(2)		2.9(1)	1.7(0.2)	1.7(0.2)	1.00(0.7)	15/15
	MCS h 0.99(0.5)		1.0(0.7)	1.0(0.7)		31(44)	36(39)	14/15	MCS h 1.5(0.9)		18(26)	29(31)	29(31)	124(103)	5/15
	NEWUO 1.8(0.2)		0.99(0.2)	0.99(0.2)		1.4(3)	2.9(2)	15/15	NEWUO 7.2(19)		2.8(6)	2.1(3)	2.1(3)	3.5(4)	15/15
	lmm-C 3.5(1)		2.5(1)	2.5(1)		1.3(0.1)	4.3(9)	11/15	lmm-C 1.9(3)		8.2(5)	408(548)	408(296)	∞ 8823	0/15
	SMAC-7.5(4)		4.2(7)	4.2(8)		2.7(4)	5.2(3)	4/15	SMAC-1.5(2)		5.0(4)	46(68)	46(75)	∞ 2000	0/15
	fminc 0.81(0.3)	0.64(0.2)	0.64(0.3)		0.96(0.9)	0.92(1)	15/15	fminc 4.9(5)		4.8(7)	3.2(2)	3.2(2)	4.4(4)	15/15
	f22		6.3e+1:45	4.0e+1:68	4.0e+1:68		1.6e+1:231	6.3e+0:1219	15/15		f24	2.5e+2:208	1.6e+2:918	1.0e+2:6628	6.3e+1:9885	4.0e+1:31629	15/15
	EE	2.4(0.7)		2.8(3)	2.8(3)		2.6(3)	2.9(3)	2/15	EE	0.40(0.2)	∞	∞	∞	∞ 501	0/15
	GLOBA 7.9(0.6)		5.8(0.9)	5.8(0.9)		2.0(0.9)	1.0(1)	15/15	GLOBA 31(95)		∞	∞	∞	∞ 3e4	0/15
	MCS h 1(0.3)		0.99(0.4)	0.99(0.3)		33(76)	20(27)	13/15	MCS h 1.8(1)		2.7(4)	12(33)	∞	∞ 8e4	0/15
	NEWUO 1.5(0.2)		1.2(0.3)	1.2(0.4)		1.00(1)	1.4(3)	15/15	NEWUO 1.1(2)		1.2(2)	4.3(5)	247(334)	∞ 2e5	0/15
	lmm-C 3.1(2)		5.8(2)	5.8(11)		5.9(13)	3.8(8)	11/15	lmm-C 0.74(0.2)	1.1(0.5)	1.4(0.8)	1.2(0.9)	1.2(1)	3/15
	SMAC-10(13)		7.2(22)	7.2(9)		4.2(7)	2.0(2)	7/15	SMAC-0.65(1)		10(5)	∞	∞	∞ 2000	0/15
	fminc 0.99(0.4)	2.5(0.2)	2.5(7)		2.9(4)	3.5(9)	15/15	fminc 0.51(0.3)	12(1)	886(1345)	∞	∞ 4e5	0/15
	Data produced with COCO v2.4
	produced with COCO v2.4

Table 2 :

 2 Continuation of Table 1.

	Fraction of function,target pairs	0.0 0.2 0.4 0.6 0.8 1.0	0 bbob-largescale f16, 320-D 2 log10(# f-evals / dimension) 4 6 51 targets: 100..1e-08 1, 15 instances v2.4 16 Weierstrass	LBFGS Var m2DLBFGS EE on bbo LMCMA14 V R10ES Var R2ES Vare sepCMA Va VkD-CMA V V2D-CMA V LMCMA17 V VD-CMA Va CMA Varel Fraction of function,target pairs	0.0 0.2 0.4 0.6 0.8 1.0	0 bbob-largescale f17, 320-D 2 log10(# f-evals / dimension) 4 6 51 targets: 100..1e-08 1, 15 instances v2.4 17 Schaffer F7, condition 10	EE on bbo LBFGS Var m2DLBFGS LMCMA14 V R10ES Var R2ES Vare LMCMA17 V V2D-CMA V sepCMA Va VkD-CMA V VD-CMA Va CMA Varel Fraction of function,target pairs	0 bbob-largescale f18, 320-D 2 log10(# f-evals / dimension) 4 6 51 targets: 100..1e-08 1, 15 instances v2.4 18 Schaffer F7, condition 1000 0.0 0.2 0.4 0.6 0.8 1.0	EE on bbo m2DLBFGS LBFGS Var LMCMA14 V R2ES Vare R10ES Var LMCMA17 V sepCMA Va V2D-CMA V VkD-CMA V VD-CMA Va CMA Varel Fraction of function,target pairs	0 1, 15 instances 2 log10(# f-evals / dimension) 4 6 v2.4 51 targets: 100..1e-08 bbob-largescale f19, 320-D 19 Griewank-Rosenbrock F8F2 0.0 0.2 0.4 0.6 0.8 1.0	LBFGS Var m2DLBFGS EE on bbo R2ES Vare R10ES Var LMCMA14 V LMCMA17 V sepCMA Va VkD-CMA V V2D-CMA V VD-CMA Va CMA Varel
	Fraction of function,target pairs	0.0 0.2 0.4 0.6 0.8 1.0	0 bbob-largescale f16, 640-D 2 log10(# f-evals / dimension) 4 6 51 targets: 100..1e-08 1, 15 instances v2.4 16 Weierstrass	m2DLBFGS LBFGS Var EE on bbo R2ES Vare R10ES Var LMCMA14 V sepCMA Va V2D-CMA V VkD-CMA V LMCMA17 V VD-CMA Va Fraction of function,target pairs	0.0 0.2 0.4 0.6 0.8 1.0	0 bbob-largescale f17, 640-D 2 log10(# f-evals / dimension) 4 6 51 targets: 100..1e-08 1, 15 instances v2.4 17 Schaffer F7, condition 10	EE on bbo LBFGS Var R2ES Vare R10ES Var m2DLBFGS LMCMA14 V LMCMA17 V V2D-CMA V VkD-CMA V sepCMA Va VD-CMA Va Fraction of function,target pairs	0 bbob-largescale f18, 640-D 2 log10(# f-evals / dimension) 4 6 51 targets: 100..1e-08 1, 15 instances v2.4 18 Schaffer F7, condition 1000 0.0 0.2 0.4 0.6 0.8 1.0	EE on bbo LBFGS Var m2DLBFGS R10ES Var R2ES Vare LMCMA14 V V2D-CMA V LMCMA17 V VkD-CMA V sepCMA Va VD-CMA Va Fraction of function,target pairs	0 1, 15 instances 2 log10(# f-evals / dimension) 4 6 v2.4 51 targets: 100..1e-08 bbob-largescale f19, 640-D 19 Griewank-Rosenbrock F8F2 0.0 0.2 0.4 0.6 0.8 1.0	EE on bbo LBFGS Var m2DLBFGS R10ES Var R2ES Vare LMCMA14 V LMCMA17 V V2D-CMA V VkD-CMA V sepCMA Va VD-CMA Va

  Hansen, R. Ros, and A. Auger. 2009. Real-Parameter Black-Box Optimization Benchmarking 2009: Presentation of the Noiseless Functions. Technical Report 2009/20. Research Center PPE. http://coco.lri.fr/downloads/download15. 03/bbobdocfunctions.pdf Updated February 2010.

In an appendix of[START_REF] Huntsman | Parallel black-box optimization of expensive highdimensional multimodal functions via magnitude[END_REF], we also produced fixed-budget plots of cumulative best values using the same data (and substituting, e.g. DTS-CMA-ES for lmm-CMA-ES) via the web interface of IOHanalyzer[START_REF] Doerr | IOHprofiler: A benchmarking and profiling tool for iterative optimization heuristics[END_REF] at https://iohanalyzer.liacs.nl/, but these plots produced neither additional nor conflicting insights.

While[START_REF] Tušar | Anytime Benchmarking of Budget-Dependent Algorithms with the COCO Platform[END_REF] points out how an anytime benchmark of a budget-dependent algorithm such as EXPLO2 can be performed with linear overhead, the cost-to-benefit ratio in our case was still prohibitive. In any event, this would not have made larger budgets any easier (or much more relevant) to obtain.

NB. For benchmarking, it was necessary to simulate parallelism, i.e., we replaced parfor loops in MATLAB code with ordinary for loops.

While *-CMA-ES is easily parallelizable[START_REF] Hansen | Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)[END_REF][START_REF] Khan | A parallel implementation of the covariance matrix adaptation evolution strategy[END_REF], the quadratic scaling of non-largescale variants with dimension creates a serious disadvantage for high-dimensional problems. In high dimensions it is appropriate to compare EXPLO2 to large-scale variants of CMA-ES[START_REF] Varelas | Benchmarking large scale variants of cma-es and l-bfgs-b on the bbob-largescale testbed[END_REF][START_REF] Varelas | A comparative study of large-scale variants of CMA-ES[END_REF], and our ad hoc experiments in this regard yielded good results.

[START_REF] Gao | Distributed Gauss-Newton optimization method for history matching problems with multiple best matches[END_REF] Also, in an experiment shown in an appendix of[START_REF] Huntsman | Parallel black-box optimization of expensive highdimensional multimodal functions via magnitude[END_REF], we found that EXPLO2 outperformed differential evolution and had performance almost indistinguishable from CMA-ES on the mixed-integer version of 𝑓 15[START_REF] Tušar | Mixed-integer benchmark problems for single-and bi-objective optimization[END_REF] in dimension 𝐷 = 20 (the only function/dimension pair we tried among the mixed-integer suite[START_REF] Tušar | Mixed-integer benchmark problems for single-and bi-objective optimization[END_REF]).

funding from the Defense Advanced Research Projects Agency (DARPA). The views, opinions and/or findings expressed are those of the author and should not be interpreted as representing the official views or policies of the Department of Defense or the U.S. Government.


Yliès Falcone

Gwen Salaün

Ahang Zuo

Probabilistic Model Checking of BPMN Processes at Runtime

Business Process Model and Notation (BPMN) is a standard business process modelling language that allows users to describe a set of structured tasks, which results in a service or product. Before running a BPMN process, the user often has no clear idea of the probability of executing some task or specific combination of tasks. This is, however, of prime importance for adjusting resources associated with tasks and thus optimising costs. In this paper, we define an approach to perform probabilistic model checking of BPMN models at runtime. To do so, we first transform the BPMN model into a Labelled Transition System (LTS). Then, by analysing the execution traces obtained when running multiple instances of the process, we can compute the probability of executing each transition in the LTS model, and thus generate a Probabilistic Transition System (PTS). Finally, we perform probabilistic model checking for verifying that the PTS model satisfies a given probabilistic property. This verification loop is applied periodically to update the results according to the execution of the process instances. All these steps are implemented in a tool chain, which was applied successfully to several realistic BPMN processes.

Introduction

A business process describes a set of structured tasks that follow a specific order and thus results in a product or service. The business process model and notation (BPMN), proposed by OMG, is the de facto standard for developing business processes [START_REF]ISO/IEC: International standard 19510, information technology -business process model andnotation[END_REF]. BPMN relies on a graphical workflow-based notation that describes the structured tasks in a business process and the relationships between these tasks.

The BPMN standard was quickly adopted by industry and academia, even though several flaws were identified. One of them regards the lack of formal semantics. Several approaches proposed to use Petri nets or automata-based languages for filling this gap. Related to formal semantics, the lack of formal analysis techniques appeared as another weakness. The final goal is to provide (ideally automated) verification techniques and tools for ensuring that processes respect some functional and non-functional properties of interest (e.g. the absence of deadlocks, the execution of the process within a reasonable amount of time, the occupancy of resources, etc.). All these checks are particularly useful for optimising processes and thus reducing the costs associated with their execution.

In this paper, we tackle the problem of computing the probability of executing certain tasks or combination of tasks when running the processes. The possibility of executing one task or another comes from the use of different kinds of gateways in the BPMN process (e.g. exclusive gateways). These probabilities are difficult to determine, especially when multiple instances of the process are executed at the same time. In that case, since resources are necessary for executing some specific tasks, knowing these probabilities is of prime importance for better adjusting the corresponding resources and thus converging to an optimal allocation of resources. It is worth noting that before executing the process multiple times, the developer has often no clear idea regarding the probability of executing some task or a specific sequence of tasks. Therefore, there is a need for automated techniques that can compute (and update) at runtime these probabilities, thus allowing the verification of probabilistic properties (e.g. what is the probability to execute task T? Is the probability to execute task T1 followed by T2 higher than 40%?).

In this work, we define an approach to perform probabilistic model checking of BPMN processes at runtime. To do so, we assume that a process is described using an executable version of BPMN. The process can be executed multiple times, each execution of the process is called an instance. Different instances may perform different tasks in the process. Our approach first monitors these executions to extract from the corresponding logs the probability of executing each individual task. These probabilities are used to build a semantic model of the BPMN process where these probabilities appear explicitly. This model is called a Probabilistic Transition System (PTS). Then, given a probabilistic property expressed in a dedicated temporal logic and this PTS, a probabilistic model checker is called for verifying whether the property is true/false or for computing the expected probability of that property. Note that this approach is not applied once and for all, because more instances of the process can keep executing including variations in terms of frequency of the executed tasks. Based on these variations, the probability of each transition of the LTS evolves over time. Therefore, the PTS is updated periodically, and the model checker is called again. The result of our approach is thus not a single value, but a dynamic curve indicating the evolution of the property evaluation over time.

To summarise, the main contributions of this work are as follows:

• Monitoring techniques for extracting at runtime relevant information about the execution of multiple instances of a process. • Periodic computation of a Probabilistic Transition System by analysing execution logs resulting of the monitoring of the process. • Integrated toolbox for probabilistic model checking of BPMN processes at runtime. • Validation of our approach on a large set of realistic BPMN processes.

The remainder of this paper is organised as follows. In Section 2, we describe the concepts and definitions used in the subsequent sections. In Section 3, we present the approach in detail. Section 4 focuses on the tool support and the experiments performed for validation purposes. Section 5 describes related work. Finally, in Section 6, we present our conclusions and future work.

Models

In this section, we introduce the preliminary concepts.

BPMN. Business process model and notation (BPMN) is a workflow-based notation for describing business processes [START_REF]ISO/IEC: International standard 19510, information technology -business process model andnotation[END_REF]. Originally, it was a modelling notation, but recent frameworks also allow the execution of such processes using a process automation engine or by translation to an executable language. The syntax of a BPMN process is given by a graph-based structure where vertices (or nodes) correspond to events, tasks and gateways, and edges (or flows) connect these nodes.

Figure 1 describes a fragment of the BPMN notation showing the main elements. Events include the initial/start event and the end event, which are used to initialise and terminate processes. We assume there is only one start event, and at least one end event. A task is an atomic activity containing only one incoming flow and one outgoing flow. Gateways are used to describe the control flow of the process. There are two patterns for each type of gateway: the split pattern and the merge pattern. The split pattern consists of a single incoming flow and multiple outgoing flows. The merge pattern consists of multiple incoming flows and a single outgoing flow. Several types of gateways are available, such as exclusive, parallel, and inclusive gateways. An exclusive gateway corresponds to a choice among several flows. A parallel gateway executes all possible flows at the same time. An inclusive gateway executes one or several flows. The choice of flows to execute in exclusive and inclusive gateways depends on the evaluation of data-based conditions. In this paper, we consider multiple executions of a single process. Each execution is called an instance and is characterised by an identifier and the list of consecutive tasks executed by this process. We assume that a BPMN process cannot run infinitely and that each instance terminates at some point. Therefore, the list of tasks associated to an instance is always finite.

LTS. We use Labelled Transition Systems as a semantic model of BPMN processes, as described in [START_REF] Mateescu | Quantifying the Parallelism in BPMN Processes using Model Checking[END_REF][START_REF] Poizat | Checking business process evolution[END_REF][START_REF] Krishna | Checking business process evolution[END_REF].

Definition 1 (LTS). A labelled transition system (LTS) is a tuple Q, Σ, q init , ∆ where: Q is a set of states; Σ is a finite set of labels/actions;

q init ∈ Q is the initial state; ∆ ⊆ Q × Σ × Q is the transition relation.
A transition (q, a, q ) ∈ ∆, written q a -→ q , means that the system can move from state q to state q by performing action a.

PTS. We also need a more expressive model than LTS because we want to associate transitions with probabilities. We therefore rely on Probabilistic Transition Systems [START_REF] Larsen | Bisimulation through probabilistic testing[END_REF], which is a probabilistic extension of the LTS model.

Definition 2 (PTS)

. A probabilistic transition system (PTS) is a tuple S, A, s init , δ, P such that S, A, s init , δ is a labelled transition system as per Definition 1 and P : δ → [0, 1] is the probability labelling function. [START_REF]Activiti: Open source business automation[END_REF] is the probability for the system to move from state s to state s , performing action a. For each state s, the sum of the probabilities associated to its outgoing transitions is equal to 1, that is ∀s ∈ S : s ∈S P (s, a, s ) = 1.

P (s a → s ) ∈ [0,
MCL. Model Checking Language (MCL) [START_REF] Mateescu | A model checking language for concurrent value-passing systems[END_REF] is an action-based branching-time temporal logic suitable for expressing properties of concurrent systems. MCL is an extension of alternation-free µ-calculus [START_REF] Emerson | Efficient model checking in fragments of the propositional mu-calculus[END_REF] with regular expressions, databased constructs, and fairness operators. We rely on MCL for describing probabilistic properties, using the following construct [START_REF] Mateescu | On-the-Fly Model Checking for Extended Action-Based Probabilistic Operators[END_REF]: prob R is op [ ? ] E end prob, where R is a regular formula on transition sequences, op is a comparison operator among "<","≤",">","≥", "=", "<>", and E is a real number corresponding to a probability. MCL is interpreted over a PTS model.

Probabilistic Model Checking of BPMN

This section first gives an overview of the different steps of our approach. Then, we present with more details the solution for monitoring BPMN processes and the computation of a probabilistic model from the execution traces observed during the monitoring step.

Overview

Recall that before executing a process, it is unclear how often a certain task or combination of tasks are executed. This is of prime importance for adjusting the resources necessary for executing the tasks involved in a process. The goal of our approach is to analyse the multiple instances of a process at runtime to precisely measure the probabilities of executing the tasks involved in a process, and thus to evaluate automatically probabilistic properties on that process.

Our approach takes as input a BPMN process and a probabilistic property, and returns as output the verdict returned by the model checker. The verdict indicates whether the property holds on the system. Such a verdict is obtained by passing the process and the property to a model checker. This verdict is periodically updated, since the process keeps on executing, and our approach runs as long as there are new instances of the process completing. Figure 2 overviews the approach. First, we monitor and analyse the multiple instances resulting of the execution of the BPMN process. These instances are used to compute the probability of execution for each task. Then, these probabilities are added to the LTS semantic model obtained from the BPMN process, resulting in a PTS. Finally, we call a model checker to verify that the PTS satisfies the given probabilistic property. Since the process keeps running, the probability of each task and thus the PTS are periodically updated. The period is a parameter of the approach. Every time the PTS is updated, the model checker is called again. Let us now give a little more details on the three main parts of the approach. Monitoring. The monitoring part focuses on the data streams generated by the execution of the BPMN process. A BPMN process may be executed multiple times, each of its executions produces an instance. Each time a new instance completes (meaning that the process has terminated), the information about that instance execution is stored into a database. We have implemented a technique for extracting from this database the events related to a set of instances that have completed, and we convert these events into execution traces (one trace per process instance). This extraction is applied periodically, where the period can be a duration (e.g. every hour) or a fixed number of completed instances (e.g. when 100 instances have completed).

Model Generation. The first step of this part is to generate an LTS from the BPMN process. This LTS stands as a semantic model, and exhibits all possible execution paths for the given process. There are several ways to transform BPMN to LTS. Here, we rely on an existing work, which proposes to first transform BPMN into the LNT (LOTOS New Technology) [START_REF] Garavel | From LOTOS to LNT[END_REF] process algebra. Since LNT operational semantics maps to LTS, the generation of that LTS is thus straightforward. Due to lack of space, the reader interested in more details regarding the transformation from BPMN to LTS can refer to [START_REF] Mateescu | Quantifying the Parallelism in BPMN Processes using Model Checking[END_REF][START_REF] Poizat | Checking business process evolution[END_REF][START_REF] Krishna | Checking business process evolution[END_REF]. Note that this transformation from BPMN to LTS is only computed once. In a second step, by analysing the execution traces built during the monitoring stage, we compute the probabilities of executing each task involved in the process, and add these probabilities to the LTS, which thus becomes a PTS. This PTS is updated periodically, every time a new set of execution traces is provided by the monitoring techniques.

Verification. This step of the approach takes as input a probabilistic model (PTS) and a probabilistic property, and computes as output a Boolean or numerical verdict depending on the property. This check is performed by using an existing model checker (the latest version of the CADP model checker [START_REF] Garavel | CADP 2011: a toolbox for the construction and analysis of distributed processes[END_REF] in this work). Since the PTS is updated periodically, the model checker is thus called whenever this update takes place. Therefore, the final result does not consist of a single value, but all successive values are gathered on a curve, which is dynamically updated every time the model checker is called with a new PTS.

BPMN Process Monitoring

In this section, we introduce monitoring techniques for BPMN processes at runtime. These techniques are useful because a process is usually not executed only once. Instead, a process can be executed multiple times. Each execution of the process is called an instance. An instance of the process can be in one of the following states: initial means that the instance is ready to start (one token in the start event), running means that the instance is currently executing and is not yet completed, completed means that all tokens have reached end events. Tokens are used to define the behaviour of a process. When a process instance is triggered, a token is added to the start node. The tokens move through nodes and flows of the process. When a token meets a split gateway (e.g. parallel gateway), it may be divided into multiple ones, depending on the type of split gateway. On the contrary, when multiple tokens meet a merge gateway (e.g. inclusive gateway), they are merged into a single token depending on the type of merge gateway. An identifier is used to characterise a specific instance, and this identifier is associated to all nodes (e.g. tasks) executed by this instance.

Monitoring techniques (see Fig. 3 for an overview) aim at analysing the information stored in a database, and extracting for each instance the corresponding execution trace. An execution trace corresponds to a list of tasks executed by this specific instance. The order of execution of these tasks is established by using timestamps at which each task is executed. These timestamps are computed by the process execution engine (Activiti [START_REF]Activiti: Open source business automation[END_REF] in this work), which relies on a global clock. The execution trace corresponding to a specific instance can be computed only when the instance is in its completed state. the time-based strategy means that the trace extraction is performed every fixed period of time; the instance-based strategy is based on the number of instances, and the trace extraction is triggered whenever the total number of new completed instances reaches a certain value.

It is worth noting that a hybrid strategy combining these two strategies is also an option, e.g. we extract traces whenever 100 instances have completed or every hour if after one hour less than 100 instances have completed. In addition, the choice of these different strategies may have a different impact on the actual results.

There are two similar algorithms for extracting execution traces depending on the strategy. We illustrate below with the algorithm relying on the time-based strategy. The first goal of this algorithm is to extract the relevant completed instances of this process from the database. These instances are then traversed in order to generate the corresponding execution traces.

Let us now go through the algorithm to give more details. Algorithm 1 describes the execution of the time-based extraction of execution traces. The inputs of the algorithm are the process identifier (pid ), a timestamp (ts), and a time duration (td ). This timestamp is the start time of the period to identify the new instances that have completed. The output is a set of execution traces (T ). T := T ∪ I .computeSortedTasks() return T Algorithm 1 first connects to the database and retrieves all the instances corresponding to the process identifier by using function getInstances(). These instances are stored in variable I. Each instance consists of the identifier of the instance and a set of tasks (lines 1 to 2). These instances are traversed to keep only those that have completed during the last period of time (presence of an end event and completion time lower than timestamp + duration, line 4). The resulting instances are all eligible instances. For each completed instance, function computeSortedTasks() sorts the tasks using their completion times, and returns an execution trace consisting of the instance identifier and an ordered list of tasks (line 5). The algorithm finally returns the set of execution traces T .

Algorithm 1 Algorithm for extracting execution traces

The time complexity of the algorithm is O(n × m × log m), where n is the number of completed instances over a period, and m is the maximum number of tasks executed by an instance (O(m × log m) is the complexity of the timsort algorithm used for sorting tasks).

Transforming LTS into PTS

Given a BPMN process, we can generate its LTS semantic model using existing techniques such as [START_REF] Mateescu | Quantifying the Parallelism in BPMN Processes using Model Checking[END_REF][START_REF] Poizat | Checking business process evolution[END_REF][START_REF] Krishna | Checking business process evolution[END_REF]. The LTS model exhibits all possible execution paths of the input BPMN process. This generated model is non-deterministic, and it has only one final state 1 . In this section, we show how by analysing execution traces (one trace per instance) extracted during the monitoring of the process, we can extend this LTS with probabilities of execution for each transition included in this LTS model. These probabilities are added as annotations to the transitions of the LTS, which thus becomes a PTS.

Before explaining how we generate a PTS given an LTS and a set of execution traces, it is worth noting that, similarly to trace extraction, the PTS should be updated periodically as well due to the execution of multiple instances. Therefore, this part of the approach also relies on one of the two aforementioned strategies (time or instance-based strategy) for defining the period.

Algorithm 2 takes as input the LTS model corresponding to the BPMN process and a set of execution traces, and returns as output a PTS model. The main idea of the algorithm is to count the number of times each transition is executed using the information from the execution traces. This is achieved by associating a counter to each transition and by traversing the execution traces one after the other. Essentially, each time a task appears in an execution trace, we increment the counter of the corresponding transition. After traversing all execution traces, we compute the probability of executing each transition outgoing from a state by using the associated counter value. We augment the LTS model with these probabilities to obtain the PTS model. for each task ∈ T tasks do 7:

qsucc := {q ∈ Q | ∃q ∈ Qcurrent, (q, task, q ) ∈ ∆} 8:

if task = T tasks [T tasks .length() -1] then 9:

Qnext := qsucc, Qcurrent := Qnext 10: else 11:

qnext := {q ∈ Q | ∃q ∈ qsucc and q = q , (q, task, q ) ∈ ∆} 12:

Qnext := qsucc \ qnext 13:

for each (q, task, q ) ∈ ∆, q ∈ Qcurrent, q ∈ Qnext do 14:

Fpaths.append((q, task, q )) 15:

for each task ∈ T tasks .reverseOrder () do 16:

Qpre := {q ∈ Q | ∃q ∈ Qnext, (q, task , q ) ∈ ∆} 17:

for each (q, task, q ) ∈ ∆, q ∈ Qpre, q ∈ Qnext do 18:

Bpaths.append((q, task, q )) 19:

Qnext := Qpre 20:

Path := Fpaths ∩ Bpaths 21:

for each (s, a, s ) ∈ Path do cnt((s, a, s )) := cnt((s, a, s )) + 1 22: P := {(s, a, s ) → cnt((s, a, s ))/ q∈S,a ∈A,(s,a ,q)∈δ cnt((s, a , q)) | (s, a, s ) ∈ δ} return S, A, sinit, δ, P

Let us now present with more details how this algorithm for generating the PTS model works. The PTS model is first initialised, and a counter (initialised to 0) is added to each transition of the LTS model (lines 1 to 3). The algorithm starts by traversing the set of execution traces T . For each execution trace, the algorithm proceeds in three steps: (a) traversing the tasks of the execution trace, (b) finding the corresponding valid path into the LTS model, (c) increasing the value of the counters. As a final step, all execution traces are traversed for computing the probability of each transition. Q current is the set of current states in the LTS during the traversal, Q next is the set of successor states of a current state, and Q pre is the set of predecessor states of a current state. We now present these steps with more details: (a) Traversing the tasks of the execution trace (lines 5 to 14). Since the LTS may exhibit non-deterministic behaviours, this step (and the following one) computes the valid path in the LTS corresponding to an execution trace. This step relies on a forward traversal of the LTS (from initial state to final state). Each execution trace T consists of an identifier and a sequence of tasks T tasks . For each trace, these tasks are handled one after the other, and by using transitions ∆, the successor states for each current state are obtained until all tasks of the current execution trace have been traversed. We use Fpaths (Forward-paths) to record the sequence of transitions in the LTS corresponding to the execution paths of the current execution trace. (b) Finding the corresponding valid path into the LTS model (lines 15 to 20).

This step relies on a backward traversal of the LTS (from final state to initial state). Therefore, we start by reversing the sequence of tasks for the current execution trace. By using this reversed list and the final state which is stored in the last Q next of the previous step, we then traverse backwards to the initial state. We use Bpaths (Backward-paths) to record all the transitions from the final state to the initial state (lines 15 to 19). Next, we take the intersection of each element in Figure 4 illustrates the execution of the algorithm. Figure 4(a) depicts the input of the algorithm: an LTS and a set of execution traces, where the number in the first column (e.g. 1003) is the identifier of the execution trace. In this example, we assume that State 2 is the end/final state of the LTS. 1Task C -→ 4). The final state reached by this path is 4, which is not the final state of the LTS. Therefore, this path is invalid. For the other one, its final state is the final state of the LTS, and hence, this is a valid path. Figure 4(d) shows each relevant transition coming from the valid path (dashed lines) whose counter is incremented by 1. Figure 4(e) then describes the computation of the probability for each transition of the LTS. Finally, the PTS corresponding to the LTS extended with probabilities is returned.

Tool Support

In this section, we present the tool chain automating the different steps of our approach. We then illustrate the application of these tools to a case study, and end with additional experiments to evaluate performance of the tools on a set of realistic examples.

Tool

Figure 5 overviews the tool chain. First, we use the Activiti framework [START_REF]Activiti: Open source business automation[END_REF] for developing and executing BPMN processes. Activiti is a lightweight Java-centric open source tool. When running a BPMN process once or several times, all data related to these executions are stored into a MySQL database.

Beyond a BPMN process, the second input required by our approach is a probabilistic property. In this work, the property is specified using the MCL [START_REF] Mateescu | A model checking language for concurrent value-passing systems[END_REF] temporal logic, which is one of the input languages of the CADP toolbox [START_REF] Garavel | CADP 2011: a toolbox for the construction and analysis of distributed processes[END_REF]. CADP is a toolbox for the design and verification of asynchronous concurrent systems. Note that the approach can take several properties as input, not just a single one. We also use the Script Verification Language [START_REF] Garavel | SVL: A scripting language for compositional verification[END_REF] (SVL), which is convenient for automating all verification tasks, particularly when there are several properties given as input.

The VBPMN tool [START_REF] Krishna | Checking business process evolution[END_REF] is used for transforming BPMN into LTS. The generation of the PTS from the analysis of the execution traces is automated by a Python program we implemented. The property is then evaluated by calling the CADP probabilistic model checker [START_REF] Mateescu | On-the-Fly Model Checking for Extended Action-Based Probabilistic Operators[END_REF]. As a result, it returns either a Boolean or Numerical value. Since the BPMN process keeps executing (multiple instances), the PTS is updated periodically according to an update strategy defined in Section 3. Whenever the PTS is updated, the model checker is called again. The final result is thus not a single value, but a series of values, which we represent as a curve (x: time or number of instances, y: verification result). This curve is drawn using Matplotlib, which is a plotting library for the Python programming language.

Case Study

Let us illustrate our approach with the shipment process of a hardware retailer, which comes from [START_REF] Mateescu | Quantifying the Parallelism in BPMN Processes using Model Checking[END_REF]. This example, shown in Figure 6, describes a realistic delivery process of goods. More precisely, this process starts when there are goods to be shipped (E1). Two tasks are then processed in parallel (PG1), one is the packaging of the goods (T7) and the other decides whether the goods require normal or special shipment (T1). Depending on that decision, a first option checks whether there is a need for additional insurance (T2), followed by the possibility to buy additional insurance (T4) and/or fill in a post label (T5). A second option is to request quotes from carriers (T3), followed by the assignment of a carrier and preparation of the paperwork (T6). Before completing the whole process, the package is moved to a pick-up area. Probabilistic Property. For illustration purposes, we choose a property checking that the probability of executing task T4 after task T2 is less than 0.4. This is important because the choice of taking extra insurance (T4) comes with a cost, and if this decision is taken too often (e.g. more than four times out of ten), this may become a problem in terms of budget. This property is written in MCL as follows: prob true*. T2. true*. T4 is < ? 0.4 end prob. Since we use the ' ?' symbol, the model checker returns both a Boolean value (indicating whether the property is true or false) and a numerical value (indicating the probability to execute T4 after T2).

Simulation. We implemented a simulator in Java in charge of executing many instances of the BPMN process, varying the order and frequency of task executions in order to simulate a realistic operating environment. Figure 7 shows the Boolean and numerical results for a simulation consisting of 1400 instances, executed over a period of about 4 minutes, where the property is the one mentioned earlier. The update strategy used here relies on the number of completed instances. Every time there are ten completed instances, we compute again the execution probability of each transition of the LTS, generate a new PTS, and we call the model checker to obtain a new result. Figure 7 shows a variation of the truth and numerical values of the evaluated property over time. This is due to our simulator, which favours the execution of some specific tasks during some periods, resulting in the curve given in the figure.

Additional Experiments for Performance Evaluation

The goal of this section is to measure the execution times of the different steps of our approach in practice. To do so, we rely on a set of realistic BPMN processes found in existing papers and frameworks shown in Table 1. We used an Ubuntu OS laptop running on a 1.7 GHz Intel Core i5 processor with 8 GB of memory to carry out these experiments. Table 1 shows the results of these experiments. Each process is characterised by its size (number of tasks and gateways), the size of the generated PTS (number of states and transitions), and the execution time of each step is decomposed as follows:

(1) Time for transforming the BPMN process into an LTS (executed only once);

(2) Time for extracting a certain number of execution traces (100 in these experiments) from the database; (3) Time for analysing these execution traces and for computing the PTS; (4) Time for verifying the property on that PTS using the CADP model checker.

Let us now focus on the results presented in Table 1 for each step. The first step focuses on the transformation of the BPMN process to an LTS model. Table 1 shows this can be a time-consuming step compared to the other ones. This computation time depends on the structure of the BPMN process and increases with the number of parallel and inclusive gateways (in particular when they are nested). Rows 8 and 11 in the table illustrate this point. However, it is worth noting that this step of the approach is only executed once at the beginning, so this extra-time is not really a problem. The second step focuses on the computation time for connecting to the database and extracting a certain number of execution traces (100 in these experiments) from it. We can see that To conclude, these experiments show that, except for the LTS computation that might be costly, the other steps of the approach are computed in a reasonable time for realistic processes. This shows that conducting probabilistic model checking of BPMN processes at runtime is feasible. Last but not least, the sum of the times observed for steps (2), ( 3) and ( 4) could be used to obtain a lower bound value to the period of time used by the time-based strategy. Indeed, this would not make sense to use as period a value that would be smaller than this lower bound.

Related Work

In this section, we overview some existing research efforts proposing probabilistic models and analysis for BPMN.

The approaches in [START_REF] Ceballos | A probabilistic BPMN normal form to model and advise human activities[END_REF][START_REF] Ceballos | Discovering causal relations in semantically-annotated probabilistic business process diagrams[END_REF] focus on the use of Bayesian networks to infer the relationship between different events. As an example, [START_REF] Ceballos | A probabilistic BPMN normal form to model and advise human activities[END_REF] introduces a BPMN normal form based on Activity Theory that can be used for representing the dynamics of a collective human activity from the perspective of a subject. This workflow is then transformed into a Causal Bayesian Network that can be used for modelling human behaviours and assessing human decisions.

The approach in [START_REF] Durán | Stochastic analysis of BPMN with time in rewriting logic[END_REF] extends BPMN with time and probabilities. Specifically, the authors expect that a probability value is provided for each flow involved in an inclusive or exclusive split gateway. These BPMN processes are then transformed to rewriting logic and analysed using the Maude statistical model checker PVeStA. This work is extended in [START_REF] Durán | Analysis of resource allocation of BPMN processes[END_REF] to explicitly model and analyse resource allocation. This series of works allows one to compute at design time generic properties, such as average execution times, synchronisation times or resource usage, whereas the goal of this paper is to compute probabilistic properties at runtime by dynamically analysing the executions of multiple process instances.

The approach in [START_REF] Herbert | Automated evolutionary restructuring of workflows to minimise errors via stochastic model checking[END_REF] presents a framework for the automated restructuring of workflows that allows minimising the impact of errors on a production workflow. To do so, they rely on a subset of BPMN extended to include the tracking of realvalued quantities associated with the process (such as time, cost, temperature), and the modelling of probabilistic-or non-deterministic branching behaviour, and the introduction of error states. The main goal of this approach is to reduce the risk of production faults and restructure the production workflows for minimising such faults.

In [START_REF] Poizat | Checking business process evolution[END_REF][START_REF] Krishna | VBPMN: automated verification of BPMN processes (tool paper)[END_REF], the authors first propose to give a formal semantics to BPMN via a transformation to Labelled Transition Systems (LTSs). This is achieved via a transformation to process algebra and use of existing compilers for automatically generating the LTS from the process algebraic specification. Once the LTS model is generated, model checking of functional properties is possible as well as comparison of processes using equivalence checking. This work does not provide any probabilistic model for BPMN nor any kind of quantitative analysis.

In [START_REF] Herbert | Quantitative analysis of probabilistic BPMN workflows[END_REF][START_REF] Herbert | Precise quantitative analysis of probabilistic business process model and notation workflows[END_REF], the authors present a framework for modelling and analysis of business workflows. These workflows are described with a subset of BPMN extended with probabilistic nondeterministic branching and general-purpose reward annotations. An algorithm translates such models into Markov decision processes (MDP) written in the syntax of the PRISM model checker. This enables quantitative analysis of business processes for properties such as transient / steady-state probabilities, reward-based properties, and best-and worst-case scenarios. These properties are verified using the PRISM model checker. This work supports design time analysis, but does not focus on the dynamic execution and runtime verification of processes.

Statistical model checking [START_REF] Herault | Approximate probabilistic model checking[END_REF], which uses simulation and statistical methods, facilitates the generation of approximate results to quantitative model checking. Although it has a low memory requirement, the cost is expensive if high accuracy is required. In comparison, probabilistic model checking produces highly accurate results, despite the potential problem of state explosion.

Conclusion

We have presented a new approach that allows BPMN analysts to automatically carry out probabilistic model checking of BPMN processes at runtime. This approach takes as input an executable BPMN process and one (or several) probabilistic property. To evaluate this property, we build a probabilistic model (PTS) by analysing the execution traces extracted from the multiple execution of this process. This analysis allows us to annotate the LTS semantic model corresponding to the BPMN process with probabilities, thus obtaining a PTS model. Finally, we can call the probabilistic model checker with the probabilistic model and the property. Since the process keeps executing, the probabilistic model is updated periodically and the model checker is called periodically as well. Therefore, we do not return a single value as a result but a curve displaying the successive truth or numerical values returned by the model checker. Our approach is fully automated by a tool chain consisting of existing and new tools. The tool chain was applied to several realistic examples for validation purposes.

As far as future work is concerned, we first plan to take advantage of the results computed by our approach to effectively adjust resource allocation depending on the runtime analysis results. This requires having an explicit description of resources associated with tasks and dynamically modifying the resource allocation with respect to the analysis results. A second perspective is to not only analyse properties at runtime, but predict the result of the evaluation of these properties in the near future. This would allow the anticipation of changes in the resource allocation for instance. This prediction can be achieved by relying on the computed probabilistic model or by using machine learning techniques.

Fig. 2 .

 2 Fig. 2. Overview of the approach

Fig. 3 .

 3 Fig. 3. Runtime monitoring of multiple executions of a BPMN process

  Inputs: A process ID (pid ), a timestamp (ts), and a time duration (td ) Output: A set of execution traces (T ) 1: I := ∅, T := ∅ 2: I := getInstances(pid ) 3: for each I ∈ I do 4: if I .hasEndEvent() and ts < I .endts() ≤ ts + td then 5:

Algorithm 2 2 :

 22 Algorithm for transforming LTS into PTS Inputs: LTS = Q, Σ, qinit, ∆ , a finite set of execution traces T = T1, T2, . . . , Tn Output: PTS = S, A, sinit, δ, P 1: S := Q, A := Σ, sinit := qinit, δ := ∆ Path := [ ], Fpaths := [ ], Bpaths := [ ], T tasks := [ ] /* [ ] indicates an empty list */ 3: for each (s, a, s ) ∈ ∆ do cnt((s, a, s )) := 0 4: for each T ∈ T do 5: Qcurrent := {qinit}, Qnext := ∅, Qpre := ∅, T tasks := T .getTasks() 6:

  Fpaths and Bpaths, and store the result in Path. This intersection operation eliminates the invalid paths, or more precisely the invalid transitions, in Fpaths and Bpaths. Thus, the Path variable finally stores all the transitions of the LTS model corresponding to the current execution trace (line 20). (c) Increasing the value of the counters (line 21). The values of the counters for the transitions in Path are increased by 1. (d) Computing the probability of each transition (line 22). The probability of each transition is computed. To do so, the value of each transition counter is divided by the total number of transitions with the same starting state. When we have traversed all the execution traces, the algorithm finally returns the resulting PTS. The time complexity of this algorithm is O(|T | × n × |∆|), where |T | is the number of execution traces, n is the number of tasks in the longest trace, and |∆| is the number of transitions in the LTS/PTS.

  Figure 4(b) depicts an example of traversing an execution trace, where the dashed lines

Fig. 4 .

 4 Fig. 4. Example describing the execution of the algorithm

Fig. 5 .

 5 Fig. 5. Overview of the tool chain

Fig. 6 .

 6 Fig. 6. BPMN shipment process of a hardware retailer

Fig. 7 .

 7 Fig. 7. Simulation results for the shipment process

Task Initial Event End Event Split gateways: inclusive, exclusive, parallel Merge gateways: inclusive, exclusive, parallel Task Flow

Fig. 1. Excerpt of the BPMN notation

Table 1 .

 1 Experimental results for some case studies the computation time of this step is less than 0.5 seconds for all the examples in the table. The third step aims to analyse these execution traces, calculating the probabilities of each transition for building a PTS by annotating the previously computed LTS. The algorithm (and its complexity) for computing that PTS was presented in Section 3. According to our experiments, we can see that this time increases with the size of the LTS (number of transitions). It takes less than a second to compute this step for most examples and it is slightly longer for a few examples (about 4 seconds for example Multi-InclusiveV2 for instance). The final computation time corresponds to the verification of the PTS model by calling the probabilistic model checker. We can see in the table that it takes about 1 or 2 seconds for each example to make this computation.

				Characteristics			PTS	Time (s)
	No.	BPMN Process	Tasks		Gateways		States Transitions (1) (2) (3) (4)
				Exclusive Inclusive Parallel
	1	Shipment [20]	8	2	2	2	18	38	15 0.32 0.61 1.45
	2	Recruitment [7]	10	1	0	6	19	31	25 0.21 0.64 1.32
	3	Shopping [17]	22	8	2	2	59	127	50 0.26 1.12 1.68
	4 AccoutOpeningV2 [17]	15	3	2	2	20	33	31 0.25 0.71 0.84
	5	Publish [17]	4	0	2	0	16	61	22 0.42 0.58 0.77
	6	Car [17]	10	2	0	2	18	31	18 0.44 0.67 0.84
	7	Online-Shop [17]	19	7	2	0	36	74	41 0.33 0.79 1.21
	8	Multi-Inclusive [17]	9	0	6	0	47	194	42 0.23 2.32 1.39
	9	Multi-Exclusive [17]	8	6	0	0	6	9	22 0.01 0.51 0.82
	10	Multi-Parallel [17]	8	0	0	6	15	22	12 0.12 0.63 0.67
	11 Multi-InclusiveV2 [17]	12	0	6	0	141	1201	78 0.25 4.29 1.58
	12	Booking [17]	11	2	4	0	53	242	22 0.19 2.37 1.24

A final state is a state without outgoing transitions. If an LTS exhibits several final states, these states can be merged into a single one, resulting into an LTS strongly bisimilar[START_REF] Milner | Communication and Concurrency[END_REF] to the original one.

Acknowledgements. This work was supported by the Région Auvergne-Rhône-Alpes within the "Pack Ambition Recherche" programme, the H2020-ECSEL-2018-IA call -Grant Agreement number 826276 (CPS4EU), the French ANR project ANR-20-CE39-0009 (SEVERITAS), and LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01).