Kirill Antonov

Maxim Buzdalov

Arina Buzdalova

Carola Doerr

Blending Dynamic Programming with Monte Carlo Simulation for Bounding the Running Time of Evolutionary Algorithms

With the goal to provide absolute lower bounds for the best possible running times that can be achieved by (1 + λ)-type search heuristics on common benchmark problems, we recently suggested a dynamic programming approach that computes optimal expected running times and the regret values inferred when deviating from the optimal parameter choice.

Our previous work is restricted to problems for which transition probabilities between different states can be expressed by relatively simple mathematical expressions. With the goal to cover broader sets of problems, we suggest in this work an extension of the dynamic programming approach to settings in which the transition probabilities cannot necessarily be computed exactly, but in which they can be approximated numerically, up to arbitrary precision, by Monte Carlo sampling.

We apply our hybrid Monte Carlo dynamic programming approach to a concatenated jump function and demonstrate how the obtained bounds can be used to gain a deeper understanding into parameter control schemes.

I. INTRODUCTION

Running time analysis of evolutionary algorithms and other search heuristics contributes to our understanding of blackbox optimization not only by providing insights into the basic working principles that drive algorithms' performance, but also by providing lower bounds for the performance of broad classes of black-box approaches [START_REF]Theory of Evolutionary Computation-Recent Developments in Discrete Optimization[END_REF], [START_REF] Doerr | Complexity theory for black-box optimization heuristics[END_REF]. Typically expressed as black-box complexity bounds, these lower bounds can be seen as a baseline against which we can compare state-of-theart solvers, with the goal to quantify the potential of further algorithm development. By comparing black-box complexity bounds for different classes of algorithms, we obtain insight into the impact of certain algorithm features, such as their degree of parallelism [START_REF] Badkobeh | Unbiased black-box complexity of parallel search[END_REF], [START_REF] Lehre | Parallel black-box complexity with tail bounds[END_REF], their memory [START_REF] Droste | Upper and lower bounds for randomized search heuristics in black-box optimization[END_REF], [START_REF] Doerr | Playing Mastermind with constant-size memory[END_REF], their selection principles [START_REF] Doerr | Introducing elitist black-box models: When does elitist behavior weaken the performance of evolutionary algorithms?[END_REF], properties of their sampling strategies [START_REF] Lehre | Black-box search by unbiased variation[END_REF]- [START_REF] Rowe | Unbiased black box search algorithms[END_REF], etc. More refined lower bounds can be obtained by specifying a family of algorithms, the available choices for configuring it, and by studying the dependency of the running time on the choices of these configurable components. This has been classically studied in the asymptotic sense [START_REF] Witt | Tight bounds on the optimization time of a randomized search heuristic on linear functions[END_REF]- [START_REF] Oliveto | A tight lower bound on the expected runtime of standard steady state genetic algorithms[END_REF], however, some recent works provide tight bounds also for concrete problem dimensions [START_REF] Chicano | Fitness probability distribution of bit-flip mutation[END_REF], [START_REF] Gießen | Optimal mutation rates for the (1 + λ) EA on OneMax through asymptotically tight drift analysis[END_REF].

Most of these works assume a static setting of the parameters. In practice, however, it is well know that dynamic parameter settings can greatly improve the efficiency of evolutionary algorithms [START_REF] Karafotias | Parameter control in evolutionary algorithms: Trends and challenges[END_REF]- [START_REF] Eiben | Parameter control in evolutionary algorithms[END_REF]. Precise lower bounds for algorithms with dynamic parameter settings are rare, see [START_REF] Doerr | Theory of parameter control mechanisms for discrete black-box optimization: Provable performance gains through dynamic parameter choices[END_REF] for a survey of rigorous results for algorithms with dynamic parameter choices and [START_REF] Böttcher | Optimal fixed and adaptive mutation rates for the LeadingOnes problem[END_REF]- [START_REF] Lissovoi | Simple hyper-heuristics control the neighbourhood size of randomised local search optimally for LeadingOnes[END_REF] for a few exceptions that provide precise lower bounds for algorithms with dynamic parameter settings.

Proving precise lower bounds for algorithms with dynamic parameter choices is challenging for several reasons: (i) lower bounds require to make statements about all possible settings (as opposed to upper bounds, in which one concrete parameter setting is analyzed), (ii) the (often very complex and trajectorybased) dependency of the parameters on the state of the algorithm needs to be taken into account, and (iii) all bounds need to be very tight to obtain an overall bound that is precise enough to give meaningful results for concrete problem dimensions (and not just for the asymptotic behavior). In lack of precise methods to derive optimal algorithm behavior, several author teams have resorted to compute actions that maximize the step-wise progress, in hope that this "driftmaximization" is close to the true optimal parameter setting policy. This step-wise progress maximization was either done exactly [START_REF] Bäck | The interaction of mutation rate, selection, and self-adaptation within a genetic algorithm[END_REF] or via Monte Carlo simulations [START_REF] Fialho | Extreme value based adaptive operator selection[END_REF]. It is known, however, that even for the ONEMAX problem with perfect fitness-distance correlation, the drift-maximization strategy is only close to being optimal [START_REF] Doerr | Optimal parameter choices via precise black-box analysis[END_REF], but not strictly optimal [START_REF] Buskulic | Maximizing drift is not optimal for solving OneMax[END_REF], [START_REF] Buzdalov | Optimal mutation rates for the (1 + λ) EA on OneMax[END_REF]. A different approach to compute precise lower bounds for algorithms with dynamic parameter settings was therefore taken in [START_REF] Buskulic | Maximizing drift is not optimal for solving OneMax[END_REF]. For the problem of minimizing the expected running time of (1+1)-type evolutionary algorithms on ONEMAX, the authors first derive an exact method to compute the optimal dynamic parameter choices and then evaluate this formula for concrete problem dimensions using a dynamic programming approach. The work was greatly extended in [START_REF] Buzdalov | Optimal mutation rates for the (1 + λ) EA on OneMax[END_REF], by considering different (1+λ)-type heuristics and by not only computing the optimal parameter choices at each state, but also the cost of deviating from it.

Our contribution: We demonstrate in this work how the approach taken in [START_REF] Buskulic | Maximizing drift is not optimal for solving OneMax[END_REF], [START_REF] Buzdalov | Optimal mutation rates for the (1 + λ) EA on OneMax[END_REF] can be extended to settings in which we do not necessarily have a closed-form expression from which we can derive exact parameter settings. Our overall approach is still based on dynamic programming, but we replace the exact formulae that model the probabilities to transition between different states by Monte Carlo simulations. That is, we take a similar approach as in [START_REF] Fialho | Extreme value based adaptive operator selection[END_REF], but we do not simply maximize the expected drift, but use these approximated probabilities to derive strategies that minimize the expected running time -the quantity that most running time analyses focus on.

We apply our approach to the RUGGEDNESS function suggested in [START_REF] Doerr | Benchmarking discrete optimization heuristics with IOHprofiler[END_REF]. RUGGEDNESS extends the classical jump function [START_REF] Droste | On the analysis of the (1+1) evolutionary algorithm[END_REF] by introducing several small jumps that the algorithms need to make in order to find the optimal solution, see Section IV-A for a formal definition. We have chosen this example because in [START_REF] Buzdalova | Hybridizing the 1/5-th success rule with q-learning for controlling the mutation rate of an evolutionary algorithm[END_REF] we observed that classical parameter control schemes essentially fail on RUGGEDNESS, in the sense that they show worse performance than the (1 + λ) evolutionary algorithm (EA) with static mutation rate 1/n (with n denoting the dimension of the problem). Interestingly, a better performance is obtained when mutation rates are capped at 1/n instead of 1/n 2 ; see Figure 1 for an illustration. That is, despite having less flexibility, the algorithms perform better. High-performing parameter control mechanisms should not suffer from more flexibility, and we therefore see this example as an interesting use-case to identify deficiencies of state-of-the-art parameter control techniques, and to develop appropriate remedies to mitigate these.

We suggested in [START_REF] Buzdalov | Optimal mutation rates for the (1 + λ) EA on OneMax[END_REF] heatmaps to illustrate the behavior of parameter control schemes. While indeed being very useful to bound the regret (i.e., the cost of deviating from the optimal parameter setting) of a single iteration, these heatmaps turned out to be suboptimal when analyzing the whole search trajectory. We therefore suggest in this work to complement them with plots that illustrate the regret per iteration. With these plots, one can now make a clear distinction between quick cheap "hiccups" and long-term harmful deviations.

Implications for Ruggedness: The application of our method to RUGGEDNESS shows that the (1 + λ) EA with static mutation rate 1/n is around 50-75% worse than the (1 + λ) EA variant using optimal dynamic mutation rates. We also illustrate the usefulness of our regret plots by analyzing why the (1+λ) EA(A, b) with generalized 1/5-th success rule, shown efficient for other problems in [START_REF] Doerr | Simple on-the-fly parameter selection mechanisms for two classical discrete black-box optimization benchmark problems[END_REF], fails on this problem.

Impact and limitations of our work: While the precise values of the lower bound and of the individual regrets for RUGGEDNESS may be of interest to a rather specialized community only, we believe that our general approach of hybridizing Monte Carlo simulation with dynamic programming should be of interest to everyone willing to derive precise lower bounds against which sampling-based optimization algorithms can be compared.

At present, our work is restricted to settings in which states are never visited twice. In practice, this limits our work to greedy ("elitist") algorithms. Extensions to algorithms and/or problems in which states may be visited more than once forms one of the most challenging directions for future work, on which we comment in more detail in Section V.

Structure of the paper: The (1 + λ) EA framework and its main variants studied in [START_REF] Buzdalova | Hybridizing the 1/5-th success rule with q-learning for controlling the mutation rate of an evolutionary algorithm[END_REF] are introduced in Section II. In Section III we describe our hybrid dynamic programming approach using Monte Carlo simulation. In Section IV we present the results of this approach applied to derive the dynamic parameter settings which minimize the expected running time of the (1+λ) EA on RUGGEDNESS. We comment on further extensions of our approach in Section V. Section VI concludes the paper.

II. (1 + λ) EVOLUTIONARY ALGORITHMS

We present in Algorithm 1 a fairly general family of (1+λ) EAs, for the context of maximizing a "fitness" function f : {0, 1} n → R defined on bit strings. Algorithms covered by this framework are initialized by selecting their first solution candidate uniformly at random. In each iteration, λ "offspring" are sampled, by modifying the best-so-far solution x. The modification is done by unary unbiased variation operators [START_REF] Lehre | Black-box search by unbiased variation[END_REF], which -according to a characterization provided in [START_REF] Doerr | Optimal parameter choices via precise black-box analysis[END_REF] -are defined by fixing a distribution D over the set of possible mutation strengths k ∈ [0..n] := {0, 1, 2, . . . , n}. The operator first samples a mutation strength k from this distribution and then flips k uniformly selected, pairwise different bits in x. 1 The best-so-far solution "survives" and determines the center of variation in the following iteration. Ties can be broken arbitrarily; in practice and in the remainder of this work we assume uniform selection among the best offspring (if at least one of them is at least as good as the parent).

We do not specify in Algorithm 1 a stopping criterion, because we are -in line with the majority of running time analysis papers -interested in the expected optimization time, i.e., the average number of function evaluation that an algorithm performs until it evaluates for the first time a solution x * with f (x * ) ≥ f (y) for all y ∈ {0, 1} n .

In [START_REF] Buzdalova | Hybridizing the 1/5-th success rule with q-learning for controlling the mutation rate of an evolutionary algorithm[END_REF], the following (1 + λ) EA variants were studied:

(1) Shift (1+λ) EA: The standard (1+λ) EA with mutation rate p is Algorithm 1 with D(p t i ) = B(n, p) for all t ∈ N, i ∈ [1..λ]. It was argued in [START_REF] Pinto | Towards a more practice-aware runtime analysis of evolutionary algorithms[END_REF], [START_REF] Pinto | A simple proof for the usefulness of crossover in black-box optimization[END_REF] that for (1 + λ) schemes this operator should not be used in practice, since it assigns probability mass to flipping 0 bits, which cannot advance the search. 1 Examples: The two most commonly used unary unbiased variation operators are 1-bit-flips and standard bit mutation (SBM). 1-bit-flips, as used within Randomized Local Search (RLS), correspond to the 1-point distribution P[k = 1] = 1, whereas SBM (used by many evolutionary algorithms) is exactly the operator characterized by the binomial distribution B(n, p).

Two strategies were suggested to mitigate this unfavorable behavior: "shifting" all probability mass from k = 0 to k = 1 and "resampling" k until k > 0. We use the shift strategy in all algorithms considered in this work, i.e., we use the mutation operator characterized by the distribution B 0→1 (n, p), which has

P[ = 1 | ∼ B 0→1 (n, p)] = P[ ≤ 1 | ∼ B(n, p)] and P[ = k | ∼ B 0→1 (n, p)] = P[ = k | ∼ B(n, p)] for k ≥ 2.
Where p is not specified, we tacitly assume p = 1/n.

(

) (1 + λ) EA(A, b): 2
This algorithm, analyzed in [START_REF] Doerr | Simple on-the-fly parameter selection mechanisms for two classical discrete black-box optimization benchmark problems[END_REF], uses a multiplicative parameter update scheme inspired by the 1/5-th success rule [START_REF] Rechenberg | Evolutionsstrategie[END_REF]. Our version uses shift mutation (as introduced above) instead of resampling mutation, but the update of p remains the same as in [START_REF] Doerr | Simple on-the-fly parameter selection mechanisms for two classical discrete black-box optimization benchmark problems[END_REF] with A = 2 and b = 1/2. That is, p is initialized as 1/n. At the end of each iteration p is updated to 2p if max{f (y 1 ), . . . , f (y λ )} ≥ f (x) and to p/2 otherwise.

(3) (1+λ) EA r/2,2r : This algorithm was suggested in [START_REF] Doerr | The (1 + λ) evolutionary algorithm with self-adjusting mutation rate[END_REF]. It creates y 1 , . . . , y λ 2 by the shift mutation with mutation rate p/2 and the other offspring, y λ 2 +1 , . . . , y λ by the shift mutation with mutation rate 2p. Let f h = max{f (y 1 ), . . . , f (y λ 2 )} and

f d = max{f (y λ 2 +1 ), . . . , f (y λ )}. If f h > f d , p is updated to p/2 with probability 3/4 and to 2p otherwise. If f d > f h , p is updated to 2p with probability 3/4 and to p/2 otherwise. If f d = f h , p is updated to either p/2 or 2p equiprobably.
(4) (1 + λ) HQEA: This algorithm was suggested in [START_REF] Buzdalova | Hybridizing the 1/5-th success rule with q-learning for controlling the mutation rate of an evolutionary algorithm[END_REF]. It hybridizes reinforcement learning with the multiplicative update used within the (1 + λ) EA(A, b). We cannot present details, for reasons of space, and refer the interested reader to [START_REF] Buzdalova | Hybridizing the 1/5-th success rule with q-learning for controlling the mutation rate of an evolutionary algorithm[END_REF] for details. For this work, it is interesting to note that the algorithm also uses shift mutation, with a mutation rate p that is updated after each iteration. All λ offspring are sampled from the same distribution.

Capping of the mutation rate p: An advantage of the shift mutation operator is that it converges to the 1-bit-flip operator when p → 0. It can therefore be used to interpolate between global and local search. However, we typically want to maintain some probability of escaping local optima. In practice, p is therefore often capped to remain within an interval [p min , p max ]. In our work, when considering parameter control methods, we fix p max = 1/2. As discussed in the introduction, p min can have a decisive influence on the efficiency of the algorithms, as can be seen in the example on the right side of Figure 1. We study two values, p min = 1/n and p min = 1/n2 .

Scope of our work: We will focus in the following on computing a lower bound for the family of algorithms that follow Algorithm 1, but which use identical distributions to sample the λ offspring (i.e., we require that D(p t i ) = D(p t j ) for all i, j ∈ [λ]). We also restrict the algorithms to those in which D(p t ) (using the previous convention, we will from now on drop the index i) may depend on n and on f (x), but not on x, 2 nor on t, nor on any other information.

III. ALGORITHM DESCRIPTION

In this section we describe our algorithm that computes good choices of distribution parameters p for each possible fitness value of the parent.

A. High-Level Description

Our algorithm is based on dynamic programming. We iterate over the possible fitness values, starting from the second best and continuing towards the smallest value. For each fitness value we aim at finding the best possible parameter value that minimizes the remaining expected running time, as well as the remaining time itself, assuming that the (1 + λ) EA will subsequently choose the optimal parameter values for all better fitness values. As a side effect we also compute the remaining expected running times for a number of parameter values, which will be later useful to evaluate regrets associated with particular parameter control schemes.

While doing it, we assume, similarly to [START_REF] Buskulic | Maximizing drift is not optimal for solving OneMax[END_REF], [START_REF] Buzdalov | Optimal mutation rates for the (1 + λ) EA on OneMax[END_REF], that for all higher fitness values the best possible expected running times are already computed. However, since we aim at dealing with various fitness functions, we use the Monte Carlo approach to approximate transition probabilities instead.

B. Requirements for the Fitness Function

Before diving into details, we discuss the limitation of this algorithm first. Our main limitation follows from our assumption that, apart from the problem size, the optimal choice of the distribution parameter p depends on the fitness value exclusively. Hence, if two individuals have the same fitness but different structure, it may result in different transition probabilities to higher fitness values. Failing to account for that may result in both overestimation and underestimation of running times. For this reason, we formulate our requirement to the fitness function as follows.

Requirement 1. For any two individuals x 1 and x 2 , such that f (x 1 ) = f (x 2 ):

• either both x 1 and x 2 shall be optima;

• or for each fitness value f > f (x 1 ), assuming Y = {y | f (y) = f } is a set of bit strings with the fitness f , there shall exist a bijective mapping π : Y → Y such that such that for each y ∈ Y the following transition probabilities shall be equal:

P[x 1 → y] = P[x 2 → π(y)].
When applied inductively, this requirement informally means that the fitness value unambiguously determines the further stochastic behavior of the evolutionary algorithm.

Note that some popular benchmark fitness functions, such as ONEMAX, as well as the function RUGGEDNESS used in this study, satisfy this requirement. Some other functions satisfy it only partially: for instance, with certain definitions of the JUMP function, the individuals that form the small fitness valley may have the same fitness but be structurally different. However, in this particular case the probability that the parent of the (1+λ) EA becomes such an individual during the optimization run is overwhelmingly small, so we can still Algorithm 2: High-level outline of the hybrid algo.

1 f min , f max ← minimum and maximum fitness values; 2 Initialize optimal times: T * fmax ← 0;

3 for f ← f max -1, . . . , f min do 4 for p ∈ {p (f ) 1 , p (f ) 2 , . . . , p (f ) m f } do 5 Compute approximate probabilities (p i ) i=0,1,...
of increasing fitness by i with mutation rate p using the Monte Carlo approach;

6 T f,p ← 1 1 -p0 1 + i>f T i • p(i-f) ; 7
Store optimal time: T * f ← min p (T f,p );

8
Store optimal rate: P opt f ← arg min p (T f,p ); 9 return {P opt , T * , T } apply our algorithm and get the results that are imprecise only up to a factor whose difference from one is negligible.

C. Detailed Description

The high-level pseudocode of the proposed algorithm is given in Algorithm 2.

We maintain T * f , our approximation of the optimal expected remaining time starting at fitness f until the optimum is hit. This value for the maximum fitness, f max that corresponds to the optimum, is obviously zero: T * fmax = 0. We compute the remaining values starting from f max -1 and stepping down until the minimum fitness value f min is processed. Since the algorithms from the (1 + λ) EA family are elitist, T * f depends only on T * f with f > f , but not on other entries. To evaluate a particular T * f , we evaluate the expected remaining times T f,p starting at fitness f until the optimum is hit, provided that (i) while the parent's fitness is f , the (1+λ) EA uses the mutation rate p, and (ii) when the parent's fitness is updated, the (1+λ) EA uses the previously computed optimal mutation rates. Since we aim at an easy-to-compute approximation scheme that does not depend on the structure of the fitness function, we employ the following simplifications.

• Instead of testing all possible p, which is computationally infeasible, we use a finite set of mutation rates {p

(f ) 1 , . . . , p (f )
m f } that may depend on f . To obtain a good approximation, this set shall be dense enough around the assumed optimal mutation rate, but in the same time it should have a small enough size so that the computation finishes in affordable times. We detail our choices of these sets later in Section III-D.

• Instead of analytically computing transition probabilities of the (1 + λ) EA from the current fitness f to each larger fitness f (which may be very hard, error-prone and computationally demanding for certain problems) we use the Monte Carlo approach, which amounts to simulation of a part of the run of the (1 + λ) EA.

The Monte Carlo evaluation of T f,p is performed as follows.

1) Choose a large number of modeled iterations, N I .

2) Choose a large number N T of successful iterations (the iterations that resulted in a strict fitness increase) to wait for: if the transitions are easy enough, N T successes already make a good estimation of the majority of likely transitions, so we do not execute all N I iterations and hence save computational resources. 3) Find an individual with fitness f and set it as the parent individual for the next iterations of the (1 + λ) EA. 4) Perform N I iterations of the (1+λ) EA (or fewer if N T is hit earlier), however, at the end of each iteration do not update the parent regardless of the outcome. Instead, for the j-th iteration we record the fitness f j of the individual that would otherwise become the next parent. Note that the correctness of not updating the parent even in the case of neutral move is motivated by Requirement 1: any individual with the fitness f induces the same behavior, so execution of a neutral move may be safely omitted.

Let N A ≤ N I be the actual number of simulated iterations. The probability estimations pi of increasing the fitness by i are computed as pi = 1

N A • |{j | 1 ≤ j ≤ N A ; f j = f + i}|.
Note that, as the (1 + λ) EA is elitist, p0 counts also the occasions where the best offspring was worse than the parent.

Finally, T f,p is computed based on pi and T * f for f > f using the equation in line 6 of Algorithm 2, which is a solution of the recurrent relation

T f,p = 1 + p0 • T f,p + i>0 pi • T * f +i .

D. Choice of Mutation Rate Sets

The choice of the mutation rate sets {p

(f ) 1 , . . . , p (f )
m f } has to be a trade-off between the accuracy of the resulting values and the computation time. However, the freedom to choose different sets for different fitness values partially reduces the effects of this trade-off. In particular, one may conduct preliminary experiments to see which probabilities for which fitness values are most promising, and increase the coverage near these probabilities in more involved experiments. Technically, this opens a possibility of a self-adaptive scheme, where first some predefined grid is used for each fitness value, and then it is refined in the most promising regions to get better results. In this paper, however, we used a more conservative scheme detailed below.

In our experiments, we used a union of two grids: the multiplicative grid of the form p i = p base •α i , which we use for large fitness values and which spans relatively small probabilities, and the additive grid of the form p i = p base + i • p step to cover the whole ranges of fitness values and probabilities. As the particular parameter values need to depend on the problem size, we give these values, as well as the Monte Carlo simulation parameters N I and N T , in the next section.

IV. EXAMPLE APPLICATION OF OUR APPROACH

In this section we outline several kinds of insights that can be derived from the results computed by the proposed algorithm. Some of them, namely lower bounds for parameter control methods, plots of optimal parameter values, and parameter efficiency heatmaps, have been previously proposed in [START_REF] Buskulic | Maximizing drift is not optimal for solving OneMax[END_REF], [START_REF] Buzdalov | Optimal mutation rates for the (1 + λ) EA on OneMax[END_REF], and the regret plots are new to this paper.

A. The Ruggedness Function

We have chosen as use-case the RUGGEDNESS function introduced in [START_REF] Doerr | Benchmarking discrete optimization heuristics with IOHprofiler[END_REF] (function F6 there). It is basically a function of concatenated jumps. That is, assuming that the unique optimum is located in some bit string z, with fitness value n, all points at Hamming distance one from z have fitness n -2, while those at distance two have fitness n -1, those at distance three have fitness n -4, those at distance four have fitness n -3, and so on. Formally, letting OM z (x

) := {i ∈ [1..n] | z i = x i }, RUGGEDNESS assigns r(x) = n if OM z (x) = n, r(x) = OM z (x) + 1 if OM z (x) ≡ n mod 2, and r(x) = OM z (x) -1 otherwise.
As mentioned above, all the mutation rate control methods considered in [START_REF] Buzdalova | Hybridizing the 1/5-th success rule with q-learning for controlling the mutation rate of an evolutionary algorithm[END_REF] perform worse on RUGGEDNESS than the shift (1 + λ) EA with static mutation rate p = 1/n, see Fig. 1. This raises two important questions, which we will answer in the remainder of this section: (i) How far are the algorithms benchmarked in [START_REF] Buzdalova | Hybridizing the 1/5-th success rule with q-learning for controlling the mutation rate of an evolutionary algorithm[END_REF] from the best possible performance? (ii) What is the impact of the individual parameter choices that are made during the optimization process?

Our experiments consistently use problem size n = 100. This choice is motivated by three main factors: it is large enough to see absolute differences in algorithms' behavior, yet small enough to finish the Monte Carlo simulations, and it allows a straightforward comparison with the previous works [START_REF] Buzdalov | Optimal mutation rates for the (1 + λ) EA on OneMax[END_REF], [START_REF] Buzdalova | Hybridizing the 1/5-th success rule with q-learning for controlling the mutation rate of an evolutionary algorithm[END_REF].

To compute the necessary data for RUGGEDNESS, we use the following parameters for our Monte Carlo approach: N I = 10 6 , N T = 5 • 10 4 . We use the same set of probabilities for each fitness, which is constructed using a multiplicative grid with the following parameters: p base = 10 -4 to match p min = 1/n 2 = 10 -4 , and α = 10 1/25 , so that 100 probabilities less than 1 are employed, whose logarithms are evenly distributed.

B. Lower Runtime Bounds for Parameter Control Methods

The near-optimal runtime values T * f for each fitness f obtained from the proposed algorithm may be directly used to construct the lower bounds on the expected running times possible for parameter control methods. For this, we compute for each fitness value f the probability p init f of hitting it with the first created individual. The value T = f T * f p init f , which is the expected runtime of the (1+λ) EA assuming near-optimal parameter choices, is then used as lower bound.

Fig. 1 presents the comparison of these lower bounds, measured in iterations, with the performance of the standard (1 + λ) EA, as well as of the (1 + λ) EA employing a few existing parameter control methods, on RUGGEDNESS and on ONEMAX for comparison, all using the shift mutation. For ONEMAX we use the exact method for computing the lower bounds from [START_REF] Buzdalov | Optimal mutation rates for the (1 + λ) EA on OneMax[END_REF]. For the parameter control methods we additionally employ two different lower bounds on the mutation rate, p min = 1/n and p min = 1/n 2 . Now we compare the insights achieved for RUGGEDNESS with the ones achieved before. In [START_REF] Buzdalova | Hybridizing the 1/5-th success rule with q-learning for controlling the mutation rate of an evolutionary algorithm[END_REF] one was able to see only that some parameter control methods are far away from (1 + λ) EA when using p min = 1/n 2 , and the situation becomes better when p min = 1/n is used. More precisely, with the less generous lower bound on the mutation rate, all the considered methods are able to perform just as (1+λ) EA. But is it possible to come up with a more efficient parameter control method that would outperform (1 + λ) EA? With the just obtained lower bounds on expected runtimes, represented with green lines in Fig. 1, one can conclude that even if one makes the optimal mutation rate choices for each fitness value, only a relatively small constant-factor improvement would be possible. This gives an important bound which was not previously known.

One of the interesting observations on RUGGEDNESS, which can be derived from Fig. 1, is that the (1 + λ) EA(A, b) has a much worse performance than both the (1 + λ) EA r/2,2r and the (1+λ) EA. However, the reason for such a behavior cannot be derived from Fig. 1. We return to this question further in this section by looking more closely at how these methods choose mutation rates and how they relate to the best choices.

2 1 2 3 2 5 2 7 2 9 10 2 10 3 λ Iterations p min = 1/n 2 (1 + λ) 2-rate (A, b) HQEA Lower 2 1 2 3 2 5 2 7 2 9 10 2 10 3 λ p min = 1/n 2 1 2 3 2 5 2 7 2 9

C. Plots of Optimal Parameter Values

The observation of running times alone cannot provide deeper insights about the problem structure, including the optimal mutation rates for different fitness values. Such a knowledge can shed some light on which issues a parameter control method may face when the (1 + λ) EA solves the problem in question. The plots of near-optimal mutation rates, which the proposed algorithm provides, may serve this purpose. What is more, one can infer the influence of the population size λ on the optimal parameters from such plots. Fig. 2 presents the plots of near-optimal mutation rates as a function of the fitness value for ONEMAX, computed as in [START_REF] Buzdalov | Optimal mutation rates for the (1 + λ) EA on OneMax[END_REF], and RUGGEDNESS. One can see that for ONEMAX these plots demonstrate a steady tendency of decreasing while getting closer to the optimum, which matches the general expectations. At the same time, most the plots for RUGGEDNESS, while having a similar general trend, show regular oscillations with a period of 2. Indeed, a small mutation rate is better for the (1 + λ) EA to improve from the fitness value that has a different parity than n (by flipping a single bit), but at least two bits need to be flipped to improve from the fitness value with the same parity as n, which requires a larger mutation rate that tends to 2/n as one gets closer to the optimum.

Note that such optimal parameters may be difficult to be tracked precisely by most parameter control methods that assume a slow change of the optimal mutation rate.

D. Parameter Efficiency Heatmaps

Since Algorithm 2 provides expected runtimes T f,p for each fitness value and mutation rate, assuming that for higher fitness values the optimal mutation rates are chosen, we can use this data to estimate the relative efficiency of mutation rate choices for each fitness value. In Fig. 3 we present this information as a heatmap for n = 100 and λ = 512. Each cell of a heatmap corresponds to a pair of f and p, where the following p are chosen: p = 10 -4+i/25 for 0 ≤ i ≤ 100.

Colors of the cells represent the relative efficiency of the corresponding p among all mutation rates for the corresponding f . They are derived from the value C f,p = exp(α f • (T * f -T f,p )), where α f ≤ 1 is chosen in such a way that at least a half of C f,p are at least 0.5. This is done to display enough information about at least 50% of possible choices, while not artificially emphasizing differences if they are negligible.

Using Fig. 3, one can make the following observations about dynamic mutation rates for the (1 + λ) EA on RUGGEDNESS:

• Overall, near-optimal mutation rates start with a rather high value at the beginning of the optimization and show a trend to gradually decrease towards the optimum. • This trend is non-monotone. By comparison with ONE-MAX, one can see that for even fitness values RUGGED-NESS shows behavior similar to ONEMAX. However, for odd fitness values the best mutation rates are higher.

These observations agree well with the findings from the previous subsection and Fig. 2. However, one can also see that in most of fitness values, except for a few last ones, there is a wide enough range of nearly equally good parameters for each fitness value, which shows that most parameter control methods have, in theory, a fair chance of maintaining good mutation rates until the last moments. sonably good values often enough, which ensures some progress, with the (1 + λ) EA r/2,2r doing slightly better.

• For last few fitness values, both algorithms visit odd fitness values exclusively, but the (1 + λ) EA r/2,2r still chooses good mutation rates, whereas the (1 + λ) EA(A, b) gets stuck with the too small rates.

The latter observation clearly points at the behavior that causes the poor performance of the (1 + λ) EA(A, b) seen at Fig. 1. It appears that the (A, b) adaptation rule forces a switch to local search as progress slows down, which converges the mutation rate to p min = 1/n 2 . With such rate it makes the chance of finding two-bit improvements Θ(n 2 ) smaller, which explains the large performance gap observed in Fig. 1.

E. Regret Plots

Overlaying parameter choices of particular parameter control methods over the parameter efficiency heatmaps may indicate the presence of problems in a method, and how large the deviation is. However, it cannot show how much of the performance the method loses from acting suboptimally.

A possible solution is to show how much the method loses to the optimal choice in each iteration. More precisely, for each iteration we take the parameter value p chosen by the parameter control method in question that corresponds to the current parent's fitness value, and plot the value |T f,p -T * f | for that iteration. Such regret plots are presented in Fig. 4. We change the plot color whenever the fitness value changes, so that one can see which iterations share the same fitness value.

We can observe in the example runs in Fig. 4 that both methods initially quickly reduce their regret. However, the (1 + λ) EA r/2,2r keeps it at a low level, whereas the (1 + λ) EA(A, b) shows a more complicated behavior. First it features two peaks of local regret maxima, which match two occasions of too high mutation rates seen in Fig. 3. It then quickly moves through a good region and falls down to small mutation rates, where it spends most of its time with very large regrets, whose sum dominates the total running time.

Note that, occasionally, a parameter control method may enter a region where T f,p reaches too high values, or where it is even infinite (which happens twice in Fig. 4). If the method is quick enough to return to better mutation rates, this does not harm the performance, as is the case for the (1 + λ) EA r/2,2r .

V. GENERALIZATION OF THE HYBRID APPROACH

As discussed in the introduction, our hybrid approach is at the moment restricted to settings in which states are not visited more than once. This can be resolved, at the price of computational complexity, by construction of Markov chains on all states with equal fitness and solving the resulting system of equations. However, this may also require a revision of used definitions, such as the one for regret: if the mutation rate is fixed to some value for one of the states, globally suboptimal choices for other states may actually make the runtime smaller.

Our approach is, in fact, not limited to (1 + λ) type algorithms: as it simulates an iteration of the algorithm as a parameterized black-box, a wide range of algorithms may be modeled, which is much harder with exact approaches.

The number of states may become too large to have them considered explicitly for certain problems, as well as for population-based algorithms. While we see no obvious solution for this problem yet, state-merging approaches might be a solution for some of the cases.

VI. CONCLUSIONS

We have extended in this work the dynamic programming approach for computing optimal state-dependent, dynamic mutation rates suggested in [START_REF] Buskulic | Maximizing drift is not optimal for solving OneMax[END_REF], [START_REF] Buzdalov | Optimal mutation rates for the (1 + λ) EA on OneMax[END_REF] to settings in which the transition probabilities cannot necessarily be computed by closed-form expressions, but where they need to be approximated by Monte Carlo simulations. We have applied this approach to derive optimal parameter choices for the (1 + λ) EA on RUGGEDNESS. We have also introduced regret plots, which not only show deficiencies in parameter control methods, but indicate their impact on the running time.

We plan to extend our work to more complicated setting detailed in Section V, as well as to use the results presented in this work to improve some of the parameter control methods.

Fig. 1 .

 1 Fig. 1. Lower bounds for the expected running times of the (1 + λ) EA using binomial distribution with varying p for mutation, compared with the performance of various parameter control methods averaged over 100 independent runs. Problem size is n = 100 for all four plots. Legend entries are: (1 + λ): the (1 + λ) EA with fixed p = 1/n, 2-rate: the (1 + λ) EA r/2,2r , (A, b): the (1 + λ) EA(A, b), HQEA: the (1 + λ) HQEA.

Fig. 2 .

 2 Fig. 2. Best mutation rates for different values of λ and fitness values. Problem size is n = 100 for both plots.

Fig. 3 Fig. 3 .( 1 +Fig. 4 .

 3314 Fig. 3. Parameter efficiency heatmaps for (1 + λ) EA, n = 100, λ = 512. Heatmap for RUGGEDNESS also contains traces of mutation rate choices in runs of the (1 + λ) EA r/2,2r (black) and the (1 + λ) EA(A, b) (red), both with p min = 1/n 2 . Color denotes the relative efficiency, the larger the better.

This assumption is irrelevant for our use-case, but the structure of x can have an important impact in general, as is easily seen for the LEADINGONES function, for which f (0 . . . 0) = f (01 . . . 1) = 0, but the optimal mutation rate for (0 . . . 0) is one, whereas the optimal one for (01 . . . 1) is p min .

Acknowledgments. The reported study was funded by RFBR and CNRS, project number 20-51-15009, and by the Paris Ile-de-France region.


Tanguy Godquin

Morgan Barbier

Chrystel Gaber

Jean-Luc Grimault

Jean-Marie Le Bars

A semantic approach for comparing Fog Service Placement Problems

Keywords: IoT Security, Service Placement Problem, Ontology, Monadic Existential logic of the Second Order

As fog computing has developed, new challenges have emerged regarding the placement of services within the network. Such issues are becoming especially important since the integration of IoT infrastructure and the relevant security concerns. There is significant variation in the definitions of service placement problems covered in the literature. Comparing these placements is a difficult task since the services, their placement objectives, their relationships with each other and the infrastructure of the network, strongly contribute to this diversity of definitions. This article proposes a methodology for comparing service placement algorithms and solutions supported by an ontology of service placement problems. This ontology is therefore used in our methodology to solve the problems annotated therein and to compare the results of the various algorithms as well as the quality of their solutions.

I. INTRODUCTION

Securing IoT devices may not always be achievable due to both their limited capabilities and restricted access imposed by the manufacturer. An alternative to securing every single device is to deploy security services at the edge of the network [START_REF] Godquin | Applied graph theory to security: A qualitative placement of security solutions within IoT networks[END_REF], following the fog computing paradigm [START_REF] Bonomi | Fog computing and its role in the internet of things[END_REF], [START_REF] Bonomi | Fog Computing: A Platform for Internet of Things and Analytics[END_REF] which leads to Service Placement Problems (SPP) [START_REF] Salaht | An overview of service placement problem in Fog and Edge Computing[END_REF] or Fog Service Placement Problems (FSPP) [START_REF] Skarlat | Optimized IoT service placement in the fog[END_REF]. Brogi et al. and Salaht et al. suggest [START_REF] Salaht | An overview of service placement problem in Fog and Edge Computing[END_REF], [START_REF] Brogi | How to Place Your Apps in the Fog -State of the Art and Open Challenges[END_REF] classifying the algorithms for solving FSPPs into three categories: algorithms with a mathematical approach (mainly focused towards Integer Linear Programming (ILP) [START_REF] Skarlat | Optimized IoT service placement in the fog[END_REF], [START_REF] Skarlat | Towards QoS-Aware Fog Service Placement[END_REF] and Mixed Integer Linear Programming (MILP) [START_REF] Doriguzzi-Corin | Dynamic and Application-Aware Provisioning of Chained Virtual Security Network Functions[END_REF]), search algorithms and heuristics (often based on greedy algorithms or heuristic behaviours [START_REF] Gupta | iFogSim: A Toolkit for Modeling and Simulation of Resource Management Techniques in Internet of Things, Edge and Fog Computing Environments[END_REF]- [START_REF] Xia | Combining hardware nodes and software components ordering-based heuristics for optimizing the placement of distributed IoT applications in the fog[END_REF]) and exotic algorithms (such as algorithms based on deep learning [START_REF] Tang | Migration Modeling and Learning Algorithms for Containers in Fog Computing[END_REF], graph theory [START_REF] Godquin | Applied graph theory to security: A qualitative placement of security solutions within IoT networks[END_REF] or even game theory [START_REF] Zhang | Computing Resource Allocation in Three-Tier IoT Fog Networks: A Joint Optimization Approach Combining Stackelberg Game and Matching[END_REF]). The problems addressed in works on FSPPs follow a similar structure divided into three models (the infrastructure, application and deployment models [START_REF] Salaht | An overview of service placement problem in Fog and Edge Computing[END_REF]). However these problems are not identical, which complicates the comparison of solutions and the methods used. In addition, the objectives pursued in these works are also different. Some focus on improving QoS [START_REF] Skarlat | Towards QoS-Aware Fog Service Placement[END_REF], [START_REF] Yousefpour | QoS-aware Dynamic Fog Service Provisioning[END_REF], average response time [START_REF] Jemaa | QoS-Aware VNF Placement Optimization in Edge-Central Carrier Cloud Architecture[END_REF], the overall security of the network [START_REF] Godquin | Applied graph theory to security: A qualitative placement of security solutions within IoT networks[END_REF] or even the resource usage [START_REF] Donassolo | Fog Based Framework for IoT Service Provisioning[END_REF]. Such diversity also occurs within the interactions among services and the network infrastructure where the number of devices, their type and characteristics fluctuate widely.

In this article, we formalize FSPPs using Monadic Existential logic of Second-Order (MESO). In order to use and manipulate the proposed model, this formalization is then expressed within a semantic structure using an ontology which allows logical relationships. Finally, we provide a methodology for solving and comparing FSPPs solutions based on this ontology, which enables the comparison of placement solutions according to different resolution methods on the same infrastructure.

The structure of the article is as follows. Section II presents the FSPPs. Section III introduces ontologies. Section IV presents our methodology for solving and comparing service placement solutions.

II. THE FOG SERVICE PLACEMENT PROBLEM (FSPP)

A. Defining a FSPP According to Salaht et al. the definition of a service placement problem in the IoT involves defining the application, infrastructure and deployment models [START_REF] Salaht | An overview of service placement problem in Fog and Edge Computing[END_REF].

The infrastructure model describes the IoT network and covers the resources (physical devices such as sensors, actuators, servers, smartphones, etc.), their characteristics (such as processor frequency, memory, etc.), as well as those of the network (latency, bandwidth, etc.).

The application model, on the other hand, describes the relationships between services in the IoT network. Services can be either independent (monolithic) [START_REF] Salaht | An overview of service placement problem in Fog and Edge Computing[END_REF], [START_REF] Skarlat | Optimized IoT service placement in the fog[END_REF], [START_REF] Skarlat | Towards QoS-Aware Fog Service Placement[END_REF], [START_REF] Jemaa | QoS-Aware VNF Placement Optimization in Edge-Central Carrier Cloud Architecture[END_REF], interdependant [START_REF] Doriguzzi-Corin | Dynamic and Application-Aware Provisioning of Chained Virtual Security Network Functions[END_REF], [START_REF] Xia | Combining hardware nodes and software components ordering-based heuristics for optimizing the placement of distributed IoT applications in the fog[END_REF] or all interacting [START_REF] Salaht | An overview of service placement problem in Fog and Edge Computing[END_REF], [START_REF] Donassolo | Fog Based Framework for IoT Service Provisioning[END_REF].

The last model, the deployment model, defines the conditions for deploying a service (see Section II-B2).

While the infrastructure model is linked to the network, the application and deployment models are specific to the services.

B. Formalizing a Fog Service Placement Problem

A Fog Service Placement Problem (FSPP) refers to identifying on which IoT devices security services should be deployed given specific conditions.

1) Modelling an IoT network: Finite Model Theory (FMT) allows defining the infrastructure with a structure M = (D, R) where D refers to the domain including all the elements of the structure; here it contains all the devices within the IoT network. The symbol R stands for the vocabulary, it contains the relationships and functions affecting the elements of D (devices of the IoT). R may also include constants 978-3-903176-32-4 © 2021 IFIP (zero arity functions) as well as relationships (=, < and >) comparing values to these constants.

This type of structure provides great flexibility without imposing limitations on the amount of information which can be contained. As an example, one may define the unary function W : D → R returning the processing power of a device and the binary function B : D2 → R returning the bandwidth from one device to another.

2) Modelling a service placement: The domain elements which will support each service are defined by a unary relation variable. Let S be the set of these unary variables, we assume that these choices may be defined with a first-order formula Φ over the enhanced vocabulary R ∪ S. Hence, the IoT security service placement problem can be expressed by a sentence of Monadic Existential Second-Order logic (MESO) as ∃SΦ(R, S).

a) Definition of a constraint: Let S * be an interpretation of S in D, a constraint C i is an application from a subset of (M, S * ) to {0, 1}.

b) Definition of an objective function: An objective function is expressed by a function O (M, S * ) → R. It evaluates and compares placement solutions and thereby defines a placement optimization problem. When several solutions are available, the objective function can be used to evaluate which solution is the best according to its evaluation criteria. The optimization problem can thus be solved either by minimizing or maximizing the value returned by this function.

3) Example on a toy model: Consider two security services S 1 and S 2 for which we define the following constraints: C 1 : Each service is deployed on a single device. C 2 : The S 1 service requires a power of 10 to operate. C 3 : Neither service should be hosted on devices interacting. C 4 : Each service must be reachable by a device with a minimum bandwidth capacity of 50. The numerical values given are constants of the problem and can therefore be changed according to the actual data. Similarly, the listed constraints can be changed to match the services to be deployed. In this example, the constraints are expressed as the following first-order formulas where E(x, y) is a binary relationship on D indicating whether x and y are linked together:

C 1,i : ∃x∀y S i (x) ∧ S i (y) → x = y. C 2 : ∀x S 1 (x) → W (x) ≥ 10. C 3 : ∀x∀y (S 1 (x) ∧ S 2 (y)) → ¬E(x, y). C 4,i : ∀x∃y (S i (x) ∧ E(x, y)) → B(x, y) > 50.
The placement problem so described can be expressed using the following MESO formula:

∃S 1 ∃S 2 C 1,1 ∧C 1,2 ∧C 2 ∧C 3 ∧ C 4,1 ∧ C 4,2 .
Note that MESO sentences express NP-complete problems as well as easy problems.

III. PRESENTATION OF ONTOLOGIES

Semantic structures seem appropriate to model our formal structure. They have been increasingly used in conjunction with IoT technologies for several years now. The scientific community refers to those structures as SWoT, an acronym that stands for Semantic Web of Things [START_REF] Gyrard | LOV4IoT: A Second Life for Ontology-Based Domain Knowledge to Build Semantic Web of Things Applications[END_REF] and are commonly expressed using ontologies.

An ontology is a formal and explicit description of concepts within a specific field that can model the structure of a system, i.e. the relevant entities and relationships that arise from its observation [START_REF] Guarino | What Is an Ontology?[END_REF].

A. Established IoT ontologies

The most commonly used ontology in SWoT research is the Semantic Sensor Network ontology (SSN1 ). This ontology defines all the concepts surrounding sensors and their observations. While the SSN ontology has a high level of detail in its design, it mainly covers sensor-based applications. However, the Devices class of SSN is particularly interesting as it allows a Device to be a composite of several smaller Devices.

Studies conducted by Bermudez-Edo et al. led to the development of the IoT-Lite ontology, a lightweight ontology that semantically describes the IoT [START_REF] Bermudez-Edo | IoT-Lite: A lightweight semantic model for the internet of things and its use with dynamic semantics[END_REF] with fewer axioms than other ontologies such as SSN, Fiesta-IoT [START_REF] Agarwal | Unified IoT ontology to enable interoperability and federation of testbeds[END_REF], OneM2M 2 , M3 and M3-lite [START_REF] Gyrard | Cross-Domain Internet of Things Application Development: M3 Framework and Evaluation[END_REF] or IoT-O [START_REF] Seydoux | IoT-O, a Core-Domain IoT Ontology to Represent Connected Devices Networks[END_REF]. The authors have shown that the axiom number of an ontology has a significant impact on its use in an IoT context where processing times are critical and capacity might be limited. One of the main benefits of this ontology lies in its lightweight, as well as in the defined concepts. IoT-Lite uses the definition of "Device" from SSN and describes a "Service" class with associated properties.

B. Related work for service placement based on ontologies

Allocation and, especially, the placement of services is highly dependent on the environment in which they are deployed. Having an extensive knowledge of the infrastructure is important for identifying the best location to deploy a security service. Ontologies enable the infrastructure to be described formally and therefore provide structure to the collected information.

The work of Tao et al. [START_REF] Tao | Multi-layer cloud architectural model and ontology-based security service framework for IoT-based smart homes[END_REF] illustrates the value of semantic structures for representing the infrastructure model. This study introduces a new and interesting ontology, however, it does not fit our application and has several limitations. Moreover, no other ontologies were reused in its design and it does not seem to be publicly accessible or listed on the LOV4IoT platform [START_REF] Gyrard | LOV4IoT: A Second Life for Ontology-Based Domain Knowledge to Build Semantic Web of Things Applications[END_REF] thereby impeding reusability and usage with new ontologies.

The work of Nezami et al. [START_REF] Nezami | Ontology-Based Resource Allocation for Internet of Things[END_REF] is based on IoT-O [22] associated with a second ontology, Cloud-O [START_REF] Rekik | Cloud Description Ontology for Service Discovery and Selection[END_REF] forming a third ontology. The authors suggest using ontologies to allocate resources in the edge computing paradigm. They propose a methodology with four steps, resource discovery, modelling, selection and allocation.

Among the methodologies in the literature that rely on semantic structures, the work of Petrovic and Tosic stands out particularly [START_REF] Petrovic | SMADA-Fog: Semantic model driven approach to deployment and adaptivity in Fog Computing[END_REF]. In their study, the authors define a structural component, called SMADA-fog, which automatically generates code for the management of fog computing infrastructure. The authors use meta-model which defines the structure and constraints of a family of models on two different aspects of deployment in a fog computing architecture. This modelling is the core of SMADA-fog as it allows both the annotation of the deployment model in the ontology, as well as using adaptation strategies through code generation for adapting to changes in the environment.

Smada-fog is particularly interesting for its use of semantic structures for automated code generation. This operation is similar to service deployment principles due to the customizable code that can be deployed. The work of Petrivic and Tovic [START_REF] Petrovic | SMADA-Fog: Semantic model driven approach to deployment and adaptivity in Fog Computing[END_REF] proposes a highly in-depth approach on implementation demonstrating the benefits of semantic structures for complex data representation. SMADA-fog, however, meets its limits in the definition of its semantic structure and its optimization model since we do not know if it is interoperable with other ontologies. Remark that the ontology used is not listed in LOV4IoT nor publicly available. Moreover, this structure relies solely on linear optimization and the reasoning capabilities of ontologies to perform its code deployment.

IV. SOLVING AND COMPARISON METHODOLOGY

Comparing placement solutions requires both context awareness and understanding placement problems related to these solutions.

We use our ontology to build a knowledge base and act as a cornerstone in our methodology for solving placement problems and comparisons of their solutions.

A. Knowledge base creation

Knowledge bases gather information about a specific topic under an accessible and interpretable format. For ontologies, a knowledge base is a set of individual instances of classes defined in an ontology.

1) Benefits of using an ontological knowledge base: The first benefit lies upon the formal part of the structure. As demonstrated by Shimizu, OWL (Web Ontology Language), commonly used for writing ontologies, is directly linked to descriptive logic [START_REF] Shimizu | Rendering OWL in LaTeX for Improved Readability: Extensions to the OWLAPI[END_REF]. Another benefit of using ontologies involves inferences. Inferences are especially interesting here since it allows, using logic, additional information to be retrieved from those that have been entered into the ontology.

2) Building the knowledge base: Building a knowledge base, such as designing an ontology, is an iterative process that can be expanded and refined. We build our knowledge base in three steps: a) Infrastructure modelling: This modelling gathers all information available for the addressed IoT network. It includes physical resources such as available equipment (PCs, servers, sensors...), their characteristics (CPU, RAM...) and those of the network. To model the infrastructure we define the M = (D, R) structure, introduced in Section II-B. Following the modelling, the model is annotated in the ontology using individual instances. b) IoT services definition: In this step, we define the services to be deployed within our IoT infrastructure. Defining services requires each service to be defined individually according to its deployment model. The deployment model of a service covers the positioning constraints of the service and at least one objective function.

We must also describe how these services interact, which is part of the application model. The application model defines whether the services are monolithic, meaning the services form independent blocks that do not interact, whether some are interdependent, or whether all the services are mutually dependent. The positioning of interdependent services is more complex.

The service definition step ends with the annotation of the services and related information into the knowledge base.

c) Design of resolution algorithms: Looking for solutions to service placement problems implies the use of resolution algorithms. Such algorithms can be of various natures [START_REF] Salaht | An overview of service placement problem in Fog and Edge Computing[END_REF] and therefore be annoted and linked to compatible services into the knowledge base. When annotating resolution algorithms into the knowledge base, the idea is to not only ensure that they are compatible with the infrastructure model, but also to indicate the services for which they can assist finding placement solutions.

B. Solving and comparing service placement problems

Our methodology is divided into five steps: information retrieval, solution identification, validation, evaluation and comparison.

a) Information retrieval: Using our knowledge base, we can interpret information about the network infrastructure, services data, their deployment model, the resolution algorithms and the set of relationships between entities described in the ontology.

b) Solution identification: We want to run the resolution algorithms compatible with the service we wish to deploy in order to identify placement solutions. These algorithms use the information extracted in the previous step and, more generally, the entire knowledge base (including constraints) to identify solutions. Depending on the algorithms, one or more solutions can be returned. Some algorithms may even be able to return the optimal solution directly, however, this means that these algorithms describe optimization concepts as discussed in the following steps.

c) Validation: The validation step is performed using constraint filtering on the solutions returned by the algorithms. Solutions that do not satisfy the constraints are not retained. No satisfactory solution may be found, in which case constraints may be weakened to offer looser placement conditions. Filtering can be redundant if some resolution algorithms already use constraints for identifying their solution. However, it remains necessary to validate the solutions returned by the algorithms. d) Evaluation: Each suitable solution is evaluated by an evaluator (most of the time an objective function) to assign each one a score. Solving service placement problems regularly involves optimization in the form of maximizing or minimizing scores. Examples of such evaluators could be functions that indicate the contribution of a solution on the QoS of the network [START_REF] Skarlat | Towards QoS-Aware Fog Service Placement[END_REF], the security provided on the network [START_REF] Godquin | Placement optimization of IoT security solutions for edge computing based on graph theory[END_REF] or the quality of the solution [START_REF] Godquin | Applied graph theory to security: A qualitative placement of security solutions within IoT networks[END_REF].

e) Comparison: The comparison step performs a score analysis of each solution and then selects the most appropriate one for service deployment. Using a single evaluator, the algorithm that returns the best result can be easily distinguished by tracing the origin of the solution that has obtained the best score according to the evaluation criteria. We define evaluators as functions that return a single numerical value which can compete with one another to compare placement solutions.

C. Overview of our ontology

The purpose of our ontology is to model service placement problems within the IoT, to unify the work on the service placement as well as to compare solutions. We named this ontology FSPlacementOntology in reference to the Fog Service Placement Problem (FSPP) defined by Skarlat et al. [START_REF] Skarlat | Optimized IoT service placement in the fog[END_REF].

For the ontology to best models service placement problems, it needs to include all the necessary models that Salaht et al. [START_REF] Salaht | An overview of service placement problem in Fog and Edge Computing[END_REF] mentioned, namely the infrastructure, application and deployment model.

In addition, to match our methodology proposed in the previous sections, the ontology must also provide representations of the resolution algorithms and placement solutions.

1) External concepts: Service placement problems, and more particularly FSPP, are mainly centred around the concept of service. We therefore decided to design our model around the Service class from the IoT-Lite ontology [START_REF] Bermudez-Edo | IoT-Lite: A lightweight semantic model for the internet of things and its use with dynamic semantics[END_REF].

We have modelled the rest of our application concepts around this class using a middle-out methodology [START_REF] Uschold | Ontologies: Principles, methods and applications[END_REF].

The IoT-Lite ontology models the service exposition of a device using the exposed and exposedBy properties connecting the Services to the Device class from the SSN ontology [START_REF]Semantic Sensor Network Ontology[END_REF].

In our application, services are hosted on the devices. Therefore, it is necessary to model the Device class, representative of them. As a result, the classes Device and System (superclass of Device) are imported into FSPlacementOntology, as well as the exposed and exposedBy properties.

2) Major concepts: The deployment model representation in ontology requires the definition of concepts related to constraints and objective functions. For this purpose, a Constraint class has been defined, defining the prerequisites to be met for a valid deployment solution.

In the interests of realism, when annotating a Service within the ontology, at least one placement constraint must be provided using the hasConstraint existential property.

A second class, Information connected to Constraint through the involves property has also been introduced to model the information used for the definition of a constraint. Among other information, this class lists DeviceCharacteristics and SpecificHardwareComponent.

The deployment model is completed by defining the ObjectiveFunction class linked to Service through the ObjectiveFunction existential property stating that a Service must always provide at least one ObjectiveFunction. As for the service dependencies, they are modelled using the dependsOn property.

The resolution algorithms are modelled within the ontology (ResolutionAlgorithm) and related to the solutions (PlacementSolution) using the findSolution universal property to match as closely as possible the intended application. This property allows resolution algorithms to return exclusively instances of PlacementSolution or to return nothing at all. Placement solutions include at least one instance of the Pairing class using the ContainPairs existential property. The latter exclusively matches an instance of Service with an instance of Device through the hasPairedService and hasPairedDevice properties.

Finally, the Device class is linked to numerous other concepts relating to device performance (CPU, Storage, Memory), CommunicationInterface or more generally DeviceCharacteristic.

3) Extending the ontology: A simplified ontology has been chosen to limit the number of axioms it contains, allowing greater flexibility for usage within the IoT context.

The benefits of this ontology lie in the ability to annotate service placement problems as well as in its versatility to describe placement constraints using semantic restrictions [START_REF] Horridge | A Practical Guide To Building OWL Ontologies Using Prot eg e 4 and CO-ODE Tools Edition 1[END_REF].

Restrictions may refer to concepts not currently described in our ontology. The idea is to provide a core ontology that can be used and enriched by the community to match various applications. Note that our assumption that the problem is expressible by a MESO formula limits the scope of our ontology. Modelling the deployment using a different logic may lead to changes in the ontology that could be considered in its evolution.

V. CONCLUSION

The aim of this paper is to generalize and formalize Fog Service Placement Problems with Monadic Existential Second-Order sentences. Our approach provides a common structure and vocabulary to express Fog Service Placement Problems, thus assisting in the design of our ontology devoted to such problems. Using our ontology, information related to the placement problems can be gathered, stored and used to infer additional knowledge, to form a knowledge base. Finally, based on this knowledge base, we propose a methodology for solving and comparing solutions to placement problems. While our comparison methodology remains invariant to the logic used in formalization, our ontology may not be suitable for placement problems that cannot be modelled according to Monadic Existential Second-Order logic. Therefore, it might be relevant to consider this category of problems as the ontology evolves.

http://www.w3.org/ns/ssn/

https://git.onem2m.org/MAS/BaseOntology


G Luo
email: luog@eeng.dcu.ie.

G-M Cloarec
email: cloarecg@eeng.dcu.ie.

J Ringwood
email: ringwood@eeng.dcu.ie.

Predictive Control of a Four-effect Falling-film Evaporator

Keywords: System Identification, Industrial Evaporator, Model Predictive Control, Neural Network Modeling, Multi-step Ahead Prediction

The paper demonstrates the use of predictive control techniques in conjunction with black-box modeling. The control strategy will be illustrated by an industrial application, a four-effect fallingfilm evaporator used in the dairy industry. Linear models are used to identify the multivariable system (MIMO). Particular attention is given to the multi-step ahead prediction behaviour of the process model. Model Predictive Control (MPC) is be used to control the evaporator, which allows for constrains to be set, both on control actions and output variables.

I. Introduction

T HE purpose of the paper will be to demonstrate the performance of Model Predictive Control (MPC) when applied to complex nonlinear multivariable (MIMO) systems. The control strategy will be illustrated and validated on an industrial plant, a four-effect (stage) falling-film milk evaporator. The system and the corresponding data sets [START_REF]DAISY: Database for the Identification of Systems, Department of Electrical Engineering[END_REF] are chosen since it has already been used as an illustrative example for black-box modeling [START_REF] Zhu | Comparison of Three Classes of Identification Methods[END_REF], [START_REF] Cloarec | Modelling of a Four-Effect Evaporator using Hybrid Neural Networks[END_REF] and control technique demonstration [START_REF] Van Wijck | Multivariable Supervisory Control of a Four-effect Falling-film Evaporator[END_REF].

In the dairy industry [START_REF] Bylund | Dairy Processing Handbook[END_REF], evaporators are used to concentrate products (such as milk) before the spray drying process. The basic principle of the four-effect falling-film evaporator [START_REF] Bylund | Dairy Processing Handbook[END_REF] consists of preheating the product before entering the first effect (using the heat of the extracted vapor from the previous effect). The product flows down as a film along the tubes of each effect, where a distribution plate is used at the top of each effect to direct the milk to the tubes. The film flows down on the evaporation tubes heated by condensing vapour on the outside of the tubes. A two-phase flow of concentrate and vapour reaches the separator where product and vapour are separated. To reduce energy consumption and heat damage, the evaporation is performed under partial vacuum.

In multiple effect (stage) systems, represented in Figure 1 [6], the evaporators are operated in series at increasing level of vacuum (decreasing temperatures). The concentrated product is pumped into the next distribution plate. The extracted vapour is compressed in thermocompressors and is used as a heat medium in the next evaporator.

The cascading of the evaporators leads to potential energy savings but also increases the performance of the evaporation while reducing the risk of damaging the product. Once the milk is concentrated, the product is directed to a dryer which can operate only if the product is sufficiently concentrated.

The four-effect falling-film evaporator is a complex nonlinear system. One of the causes of nonlinearity is the cascading of the evaporators and their capacity to hold quantities of product. By doing so, they can act as buffers leading to nonlinear behaviour, especially for the flow of product. The nonlinearities and the complexity of the input and outputs make the realization of optimal control more difficult. Modern control methods such as Model Predictive Control (MPC) are particularly suitable for performing multicriteria control and optimizing operation costs while ensuring safeguards for critical variables. MPC requires multi-step ahead prediction modeling techniques that can deal nonlinear multivariable industrial processes. The core of MPC (and the main challenge) is to obtain a reliable multi-step ahead prediction model. The nonlinear nature of the problem lead to use nonlinear model such as neural network models or alternatively, a linearised model which has good validity around the operating region of the system [START_REF] Cloarec | Modelling of a Four-Effect Evaporator using Hybrid Neural Networks[END_REF].

The data sets used to perform the modeling and control are those of an experimental four-effect falling-film evaporator at the Netherlands Institute for Diary Research NIZO and were obtained from the DAISY data base [START_REF]DAISY: Database for the Identification of Systems, Department of Electrical Engineering[END_REF]. The process is modeled as a three-input three-output process with disturbances. Other variables such as dry matter content of feed, temperature of feed, temperature of cooling water entering the condenser and surrounding temperatures are considered as disturbances.

The input signals are essentially pseudo-random binary sequences (PRBSs) and are used to ensure that the input space spans the range of operation modes. The sampling period is 10 seconds and the training and validation sets contain 3000 samples each. The model was resampled to a sampling period of 1000 seconds for control purposes in order to achieve a realistic prediction horizon of about 20 samples, which must cover the complete transient response of the system. This is vital since the controller calculation is iterative.

II. Model Predictive Control

The concept of Model Predictive Control (MPC) is to use the process model to predict future process outputs over a certain period of time called the prediction horizon, and the predictive controller calculates the future controller outputs such that the predicted process outputs (ŷ(k + 1), . . . , ŷ(k + N )) follow the set-points or the reference trajectories (w(k +1), . . . , w(k +N )) in the desired manner. The predictive controller is often regarded as an optimizer, as it attempts to minimize a cost function which involves the error between the desired output and the process output and the control actions, taking into account the constrains on the control actions.

The basic structure of the predictive controller is illustrated by Fig. 2. P A predictive controller uses the receding horizon approach, that is, instead of using the whole controller output sequence of (u(k), u(k + 1), •, u(k + N -1)) t as determined by the optimizer, only the first control signal u(k) is applied to the process. At the next sample, the optimization procedure is repeated but using the latest measured process outputs rather than the predicted outputs at the previous sample. The advantage of using the receding horizon approach is that it compensates for future disturbances and modeling errors. The operation of the MPC is summarized in Fig. 3. At time t k , the process output is compared to the setpoint, and a 'reference trajectory' (chosen by the designer) established along which it is desired to drive the process output. The predictive controller, using an optimization procedure, then determines a set of allowable control signal required to achieve this output response and applies the fist of these control actions. This procedure is repeated at t k+1 .

The cost function is usually a quadratic function, such as equation 1:

J(k) = Q N 2 i=N 1 (ŷ(k + i) -w(k + i)) 2 + R N 2 i=1 (∆u(k + i -1)) 2 (1)
where:

• N 1 and N 2 are the minimum and maximum prediction horizon, respectively,

• N u is the control horizon,

• Q and R are weighting factors,

• bl and bu are the lower and upper bounds on the control, respectively. The constrains on the command signal are defined by equation 2:

bl ≤ u(k + i -1) ≤ bu ∀i ∈ [1, N u ] (2)
The GPC (Generalized Predictive Control) approach assumes that after an interval of N u (N u < N 2 ), projected control increments are zero [START_REF] Clarke | Generalized Predictive Control: Part 1. The Basic Algorithm[END_REF], i.e.:

∆u(k + i -1) = 0 ∀i > N u (3)
For the constrained optimization problem, an analytical solution is not available and therefore an iterative optimization method is used. This iterative optimization is implemented using a quazi-Newton algorithm, which utilized a numerical gradient of the function to be optimized. The process responses and controller actions to positive step changes in the set-points are shown by figures 4,5 and 6. For simplicity, the set-points are used as the reference trajectories, i.e. the reference trajectories are constants before and after the step change. In such a case, the dynamic response can be shaped using the control weighting function. Alternatively, the reference trajectories can also be defined as first (or higher) order trajectories, and similar results can be obtained.

III. Results and Discussions

Figures 4(c) and 5(c) show that the step changes in the set-point for y 1 (DMC) and y 2 (FP) almost have no influence on y 3 (TP) i.e. decoupling of the system has been achieved. Figure 4(b) shows that the changes in set-point for y 1 (DMC) causes y 2 (FP) to response sharply and reversely, and a heavier control weighting is preferable to reduce such an impact on y 2 . However, Figure 5(a) suggests that minimum control weighting is desirable to reduce the impact of the change in set-point for y 2 on y 1 .

Figures 4(a), 5(b) and 6(c) show that the bigger the control weighting factor R is, the smoother the process output will be following set-point changes, and the smaller R is, the faster the process output will reach the desired value. In this way, the transient response of the system can be shaped.

Figures 6(a) and 6(b) show that the impact of a step change in the set-point for y 3 (TP) is considerable and, in some senses, undesirable. It is shown that a heavier control weighting is required. The set-point 3 changes at sampling time k = 5 from 46 to 46.5, an increase of 0.5C. When R = 100, y 1 (DMC) increases by 16% from 0.5 to ≈ 0.58, and y 2 (FP) increases from 0.05 to more than 0.15. However, bigger control weighting factors will cause the outputs of y 1 and y 2 to resume their set-points in a more relaxed fashion.

In short, the set-point changes in y 1 (DMC) and y 2 (FP) do not have severe impact on process outputs. The set-point change in y 3 (TP) has a severe impact on process outputs, and this change also causes much bigger controller actions than those caused by the changes in set-point 1 and 2 as shown by figures 6(d, e, f), figures 5(d, e, f) and figures 4(d, e, f). A bigger control weighting factor R is helpful to obtain more gentle controller actions.

IV. Conclusions

This paper only presents the simulation results of the preliminary stage of a project for predictive control of complex nonlinear MIMO systems of evaporators and spray dryers. Further investigations on predictive controller for nonlinear systems described by neural network models are still to be conducted, and the predictive control strategy is to be compared with other conventional control strategies.

These preliminary results, however, do give useful insight into the behaviour of the system and confirm that predictive control is a suitable control strategy to employ.

Fig. 1 .

 1 Fig. 1. Four-effect Falling-film Evaporator.

• input 1 :

 1 Flow of Feed (FF) • input 2 : Vapour Flow to First Effect Evaporator Stage (VF) • input 3 : Cooling Water Flow to Condenser (CWF) • output 1 : Dry Matter Content of the Product (DMC) • output 2 : Flow of Product (FP) • output 3 : Temperature of Product (TP)

Fig. 2 .

 2 Fig. 2. Block Diagram of a Predictive Controller.

Fig. 3 .

 3 Fig. 3. Operation of a Predictive Controller.

Fig. 4 .

 4 Fig. 4. Process Response and Control Actions to setpoint change in DMC.

Fig. 5 .

 5 Fig. 5. Process Response and Control Actions to setpoint change in FP.

Fig. 6 .

 6 Fig. 6. Process Response and Control Actions to setpoint change in TP.

The authors would like to thank Michiel van Wijck, Consultant Advanced Control, Honeywell Hi-Spec Solutions, The Netherlands.


Ana Kostovska

Jožef Stefan

Diederick Vermetten

Carola Doerr

Sašo Džeroski

Panče Panov

Tome Eftimov

OPTION: OPTImization Algorithm Benchmarking ONtology

Keywords: Computing methodologies → Ontology engineering, Reasoning about belief and knowledge, Continuous space search

Many platforms for benchmarking optimization algorithms offer users the possibility of sharing their experimental data with the purpose of promoting reproducible and reusable research. However, different platforms use different data models and formats, which drastically inhibits identification of relevant data sets, their interpretation, and their interoperability. Consequently, a semantically rich, ontology-based, machine-readable data model is highly desired.

We report in this paper on the development of such an ontology, which we name OPTION (OPTImization algorithm benchmarking ONtology). Our ontology provides the vocabulary needed for semantic annotation of the core entities involved in the benchmarking process, such as algorithms, problems, and evaluation measures. It also provides means for automated data integration, improved interoperability, powerful querying capabilities and reasoning, thereby enriching the value of the benchmark data. We demonstrate the utility of OPTION by annotating and querying a corpus of benchmark performance data from the BBOB workshop data -a use case which can be easily extended to cover other benchmarking data collections.

INTRODUCTION

Evolutionary computation, and black-box optimization in general, is a rapidly growing field, which has seen tremendous progress in the last few decades, where many new algorithms are developed daily, making it impossible for researchers and practitioners in the field to stay up to date with all these developments and difficult to share their results with the research community.

Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for third-party components of this work must be honored. For all other uses, contact the owner/author(s). GECCO '21 Companion, July 10-14, 2021, Lille, France © 2021 Copyright held by the owner/author(s). ACM ISBN 978-1-4503-8351-6/21/07. . . $15.00 https://doi.org/10.1145/3449726.3459579 There are unfortunately no common standards for which performance data to record, nor how to store it. Storing, sharing, and making the benchmark optimization data reusable is challenging, because of the existence of different data formats that are weakly compatible. There exist a large number of different benchmarking platforms for optimization, which have mostly been designed independently from each other.

One well-established solution to the problem of managing such complex data pools is the use of ontologies. Ontologies are specifications of a shared conceptualization of data from distributed and heterogeneous systems and databases, and as such, they enable data interoperability, efficient data management and integration, and cross-database search.

THE OPTION ONTOLOGY

In this paper, we present the OPTimization Algorithm Benchmarking ONtology (OPTION). The ontology we propose was designed with the main goal of standardizing and formalizing knowledge from the domain of benchmarking optimization algorithms, where an emphasis was put on the representation of data from the benchmark performance space. OPTION offers a comprehensive description of the domain, covering the benchmarking process, as well as the core entities involved in the process, such as optimization algorithms, benchmark problems, evaluation measures, etc. The OPTION ontology is publicly available on BioPortal1 the largest repository of ontologies. We have also constructed a semantic annotation pipeline for annotating performance data, as well as a querying endpoint for querying the knowledge base. Finally, we have annotated benchmark performance data from 32 optimization algorithms and show how these annotations can be queried.

USE CASE

To demonstrate the utility of the OPTION ontology for semantic annotation, we consider the BBOB benchmark suite as a use case. We use a subset of BBOB data that includes the algorithms from the 2015-2017 competitions (three years). This gives us a total of 32 algorithms, which we semantically annotate using the OPTION ontology. All annotations were produced with the Apache Jena RDF library in the form of RDF graphs. For each of the 32 algorithms, a separate RDF graph was generated and uploaded to Apache Jena TDB2, a native triple store that efficiently stores RDF data.

The data in the ontology provides the full information on the used problem/dimension/instance/algorithm, as well as the corresponding performance data. In addition, data provenance information has been manually collected and linked to the performance data to trace its origin. The stored data provenance information includes the digital object identifier (DOI) of the paper, the paper's title, the authors' name, and the year of the publication.

QUERYING THE KNOWLEDGE BASE

After populating the knowledge base, we can use the SPARQL query language to ask queries. To enable this functionality, we have set up an Apache Jena Fuseki2 server and implemented two services. The query service provides an endpoint for handling SPARQL queries in a RESTful manner, while the upload service enables the upload of RDF data into the Apache Jena TDB2 triple store. A non-exhaustive list of queries that are currently supported by the OPTION knowledge base is as follows: (1) For a given algorithm (e.g., CMA-ES) and fixed budget scenario (e.g., between 1000 -2000 function evaluations), return the noise-free fitness -Fopt calculated for a set of benchmark problems (e.g., f1 and f7) and a set of problem instances (e.g., the first 5 instances); (2) For a given algorithm (e.g., CMA-ES) and fixed target, return the number of function evaluations needed to reach the target; (3) For multiple algorithms and fixed budget return the noise-free fitness -Fopt calculated for a set of benchmark problems (e.g., f1-f5); (4) For a set of instances from a given benchmark problem, return all noise-free fitness -Fopt values for all (or a set of) algorithms with a fixed budget of 1000 function evaluations;

(5) For a given study, return the provenance information. SPARQL queries can become very complex and sometimes are seen as a bottleneck to the wider acceptance of Semantic Web technologies. SPARQL query construction is an error-prone and timeconsuming task, which requires expert knowledge of the whole stack of semantic technologies. To facilitate the process of querying and to avoid the need of directly writing SPARQL queries, we propose the use of a GUI (for which a mock-up is depicted in Figure 1) to alleviate the problem explained above. Such a GUI would provide end-users with a search tool through which they can parameterize the predefined queries, without knowing the exact semantic data model used for the annotations. The users would construct the query by providing the search parameters of interest via input fields or predefined drop-down menus that contain ontology terms and the SPARQL queries would be generated in the background.

CONCLUSIONS & FUTURE WORK

The development of the OPTION ontology is a big step forward in the process of making benchmark data more reusable and interoperable. By annotating a large subset of the BBOB-data, we have shown the potential of the ontology to aid with data integration, while simultaneously providing powerful querying capabilities for direct analysis of the required datasets. This significantly reduces the time needed to collect data across many functions and algorithms, while being flexible in the type of performance perspective (i.e., fixed-budget, fixed-target) to be used.

Our next steps are as follows. First, we will annotate more data sets, starting with those provided by IOHprofiler [START_REF] Doerr | Benchmarking discrete optimization heuristics with IOHprofiler[END_REF][START_REF] Vermetten | Integrated vs. sequential approaches for selecting and tuning CMA-ES variants[END_REF] and those of Nevergrad [START_REF] Rapin | Nevergrad -A gradient-free optimization platform[END_REF]. Further extensions include a web-based GUI, similar to what is shown in Figure 1. We believe this to be important for better adoption of our ontology by the optimization heuristics community, whose members are not necessarily familiar with SPARQL queries. This can then be further extended to include direct integration with existing post-processing tools, e.g., IOHanalyzer [START_REF] Wang | IOHanalyzer: Performance Analysis for Iterative Optimization Heuristic[END_REF], which would greatly simplify the process of robust algorithm comparison.

Figure 1 :

 1 Figure 1: A GUI mockup for querying the knowledge base.

OPTION atBioPortal: https://bioportal.bioontology.org/ontologies/ OPTION-ONTOLOGY

ACKNOWLEDGMENTS

The authors acknowledge the support of the Slovenian Research Agency through research core grants No. P2-0103 and P2-0098, project grants No. J2-9230 and Z2-1867, and young researcher grant No. PR-09773 to AK, as well as the EC through grant No. 952215 (TAILOR). Our work is also supported by the Paris Ile-de-France region.