Jinyan Liu

Wassila Ouerdane

Vincent Mousseau

A Metaheuristic Approach for Preference Learning in Multi-Criteria Ranking based on Reference Points

In this paper, we are interested in a family of multicriteria ranking methods called Ranking with Multiple reference Points (RMP). This method is based on pairwise comparisons, but instead of directly comparing any pair of alternatives, it compares rather the alternatives to a set of predefined reference points. We actually focus on a Simplified RMP model (S-RMP) in which the preference parameters include the criteria weights and the set of reference points ordered by importance. Elicitation of the parameters (from the data provided by the decision makers) leads us to the preference learning algorithms that cannot only be applied on relatively small dataset. Therefore, we propose in this work a preference learning methodology for learning S-RMP models from a large set of pairwise comparisons of alternatives. The newly proposed algorithm is a combination of an Evolutionary Algorithm and a Linear Programming approach. Empirical results and analysis are also presented.

Introduction

Decision makers (DMs) often face decision situations in which different conflicting viewpoints are to be considered. When modeling a real world decision problem using Multiple Criteria Decision Aid (MCDA) theories, several problematics can be considered: choice, sorting or ranking problem [START_REF] Roy | Multicriteria methodology for decision aiding[END_REF]. In this paper, we are interested in multi-criteria ranking problem, where the aim is to establish a preference order (or ranking) on the set of alternatives.

The field of MCDA offers a selection of methods and operational tools i.e. aggregation models that explicitly account for the diversity of the viewpoints considered. Each method constructs first a model of DM's preferences and then exploits this model in order to work out a recommendation. During the aggregation phase, different parameters are needed, such as weights, marginal value functions, thresholds, etc. depending on the method. Such parameters allow to elaborate models taking into account the DMs preferences. Therefore, it is clear that such preferences play a key role in the construction of the recommendations. In fact, they are meaningful and acceptable only if the DMs values are taken into account.

The process by which the DMs values or the parameters of an aggregation model are captured is called preference elicitation or disaggregation [START_REF] Roy | Multicriteria methodology for decision aiding[END_REF]. Preference elicitation aims at helping the analyst to appropriately elicit the DM's preferences to be represented in a decision model. The process consists of a dialogue with the DM, where the aim is to infer the values of the parameters from the holistic information given by the DM. We note that this is not an easy task, especially when the information provided contains inconsisten-cies. Moreover, the DM has in general a limited knowledge on the aggregation models and can only express his preferences in a rather intuitive and ambiguous manner.

Generally speaking, there are two paradigms of preference elicitation approaches, namely direct and indirect paradigms [START_REF] Mousseau | Eliciting information concerning the relative importance of criteria[END_REF]. In the direct aggregation paradigm, the parameter values are directly provided by the DM through an interactive communication with the analyst. The aggregation model is firstly constructed with these parameters and then applied to the alternative set to obtain the DM's preferences. Within such a paradigm, the DM should make enough effort to understand the meaning and the roles of these parameters and to associate appropriate values to them, which may be beyond his cognitive limitation. In the indirect disaggregation paradigm, the DM provides holistic preference information such as pairwise comparisons, assignment examples etc., from which a preference model is derived and then applied to contextual recommendation. In contrast with the direct elicitation, the parameter values are regressed from the DM [START_REF] Jacquet-Lagrèze | Preference disaggregation: 20 years of mcda experience[END_REF].

Nowadays, we are facing with decision problems involving large datasets. It requires adapting and improving the algorithms and technics to construct acceptable recommendations. Some of related challenges have been addressed in the scientific community of Preference Learning, which focus on the computational complexity rather than the decision problem itself [START_REF] Fürnkranz | Preference learning[END_REF]. In the perspective of preference elicitation, the emergence of applications in MCDA with the intention of dealing with large datasets provokes our research interests. In this perspective, we intend to provide an efficient algorithms to infer the parameter values of an aggregation model called S-RMP such that the DM has provided as input a large set of pairwise comparisons.

The paper is organized as follows: Section 2 introduces the RMP method and its simplified version S-RMP. We present our metaheuristic approach for S-RMP disaggregation in Section 3. The numerical analysis and the benefits of the proposed approach are provided in Section 4. At the end, we conclude the paper.

Ranking with Multiple reference Points

Recently, a ranking method called Ranking with Multiple reference Points (RMP) has been proposed in [START_REF] Rolland | Procédures d'agrégation ordinale de préférences avec points de référence pour l'aide à la décision[END_REF][START_REF] Rolland | Reference-based preferences aggregation procedures in multi-criteria decision making[END_REF]. The idea is to construct the global preference relation between two alternatives on the basis of their relative comparisons with specified reference points. This paper is concerned with learning the parameters of this method. We give in the next sections an overview of such a method as well as its simplified version named S-RMP, introduced in [START_REF] Zheng | Preference elicitation for a ranking method based on multiple reference profiles[END_REF][START_REF] Zheng | Preference Elicitation for Aggregation Models based on Reference Points: Algorithms and Procedures[END_REF]. Further work has also been presented in [START_REF] Bouyssou | Multiattribute preference models with reference points[END_REF].

General RMP method

We consider a multiple criteria ranking problem with n alternatives in the set A indexed by N = {1, 2, ..., i, ..., n} and m monotone criteria in the set F , indexed by means of a set M = {1, 2, ..., j, ..., m}. The evaluations of alternatives on a criterion j take their value in the associated evaluation scale Xj. The X denotes the evaluation space, X = ∏ j∈M Xj i.e. the Cartesian product of evalution scales. Obviously, the evaluation of any alternative a ∈ A is a vector denoted by a = (a1, a2, ..., am) ∈ X .

The RMP method is a three-step multi-criteria paradigm for ranking alternatives. It involves k reference points such that P = {1, 2, ..., h, ..., k}. The evaluation of each reference point p h , h ∈ P on a criterion j is denoted by p h j ∈ Xj. To establish a global preference relation between two alternatives, a and b, the method specifies the following three steps:

1. Compare each alternative a ∈ A (respectively, b ∈ A) to every reference points p h , h ∈ P on every criterion j, j ∈ M . 2. Aggregate the results of the step 1 considering the m criteria and deduce the preference relation between two different alternatives (a, b) which is depending on the reference point p h , also called the relative preference with respect to the reference point p h ; 3. For each pair of alternatives (a, b), aggregate the k preference relations into global preference relation.

The first step establishes the preference relation between each alternative and each reference point on each criterion j. In the second, we only consider the criteria for which a (respectively, b) is at least as good as p h . This set of criteria is denoted by C(a, p h ) (respectively, C(b, p h )) and defined as C(a, p h ) ∈ P(M ) such that C(a, p h ) = {j ∈ M | aj ⩾ p h j }. We define then the importance relation among criteria with respect to the reference point, denoted by ▶ p h , which means that a set of criteria is at least as important as another set of criteria as follows;

Definition 1. (Importance relation among criteria w.r.t. reference point)

The importance relation ▶ p h is defined on M × M such that:

1. ∀M1 ⊂ M, M1 ̸ = ∅ ⇒ M1▶ p h ∅, ∀h ∈ P 2. ∀M1 ⊂ M ⇒ M ▶ p h M1, ∀h ∈ P 3. ∀M1, M2 ⊂ M, M1 ⊂ M2 ⇒ M2▶ p h M1, ∀h ∈ P
Thus, the (relative) preference relation ≿ p h defined in Definition 2 expresses how a pair of alternatives compare with each other with respect to the reference point p h .

Definition 2. (Relative preference w.r.t. reference point)

The relative preference with respect to the reference point ≿ p h on A × A is defined by:

a ≿ p h b ⇔ C(a, p h )▶ p h C(b, p h )
In the third step, we define firstly the importance relation ⊵ as below.

Definition 3. (Importance relation among reference points)

The importance relation ⊵ is defined on P × P such that:

1. ∀P1 ⊂ P, P1 ̸ = ∅ ⇒ P1 ⊵ ∅ 2. ∀P1 ⊂ P ⇒ P ⊵ P1 3. ∀P1, P2 ⊂ P, P1 ⊂ P2 ⇒ P2 ⊵ P1
We deduce then the (global) preference relation ≿ as described in the definition 4. It means that a is at least as good as b if the coalition of reference points P (a, b) which affirms that a is at least as good as b is more important than the coalition of reference points P (b, a) which affirms that b is at least as good as a.

Definition 4. (Global preference relation)

The global preference relation on A×A on the basis of the relative preferences is:

a ≿ b ⇔ P (a, b) ⊵ P (b, a)
We note that there is no lack of generality to impose a dominance relation among reference points, as it was shows in [START_REF] Rolland | Reference-based preferences aggregation procedures in multi-criteria decision making[END_REF]. There exists always an equivalent RMP model for any RMP model with k reference points such that:

∀j ∈ M, ∀h, h ′ ∈ P, h > h ′ ⇒ p h j ⩾ p h ′ j (1)
It means also that:

∀a ∈ A, ∀h, h ′ ∈ P, h > h ′ ⇒ C(a, p h ) ⊆ C(a, p h ′ ) (2)
The two importance relations mentioned respectively in Definition 1 and Definition 3 can be rather general and built on the basis of different rules. Particularly, we present in the next section a simplified version of RMP method, namely S-RMP in which the importance relation ▶ p h is defined by a concordance rule for all reference points while a lexicography of dictatorial reference points is used as the importance relation ⊵.

Simplified RMP model

In this work, we focus on a simplified version of RMP, named S-RMP model, as considered in [START_REF] Zheng | Preference elicitation for a ranking method based on multiple reference profiles[END_REF][START_REF] Zheng | Preference Elicitation for Aggregation Models based on Reference Points: Algorithms and Procedures[END_REF].

Firstly, the importance relation ▶ p h on the criteria set is re-defined based on a concordance rule by using an additive decomposition as follow:

C(a, p h )▶ p h C(b, p h ) ⇔ ∑ j∈C(a,p h ) ωj ⩾ ∑ j∈C(b,p h ) ωj ( 3
)
where ωj represents the weight of the criterion j. Formally, the criteria weights are normalized to 1. Secondly, as shown in [START_REF] Rolland | Reference-based preferences aggregation procedures in multi-criteria decision making[END_REF], an important result derived from social choice theory indicates that the only importance relation, which aggregates the k relative preference relations (with respect to reference points) and leads to transitive relation on each possible set of alternatives, is obtained by a lexicographical order on the reference points. Therefore, a permutation σ on P is used lexicographically to determine the importance relations among the reference points so as to deduce globally the preference relations between alternatives. This is represented by [START_REF] Cailloux | Eliciting ELECTRE TRI category limits for a group of decision makers[END_REF].

a ≿ b ⇔ ∃h * ∈ P s.t. a ≿ p σ(h * ) b and ∀h ∈ {1, ..., h * -1}, a ∼ p σ(h) b a ∼ b ⇔ ∀h ∈ P, a ∼ p h b (4)
In details, the first reference point is denoted by p σ (1) , the second one by p σ (2) and so on. To compare a and b, we look at first the reference point p σ (1) . If a is strictly preferred to b with respect to p σ (1) , then we affirm globally that a is preferred to b without even considering the other reference points and similarly vice versa. However, if a and b are considered to be indifferent with respect to p σ (1) , we shall proceed with the second reference point p σ (2) , etc. Once the strict relative preference between a and b is confirmed, the global preference between a and b is affirmed. If we still cannot differentiate these two alternatives until the last reference point p σ(k) has been processed then they are considered globally as tied for conclusion.

The additive decomposition of importance relation for subsets of criteria shows that the preference relation between two alternatives (a, b) can be computed by:

a ≿ b ⇔ ∃h * ∈ P s.t. ∑ j∈C(a,p σ(h * ) ) ωj > ∑ j∈C(b,p σ(h * ) ) ωj and ∀h ∈ {1, . . . , h * -1}, ∑ j∈C(a,p σ(h) ) ωj = ∑ j∈C(b,p σ(h) ) ωj a ∼ b ⇔ ∀h ∈ P, ∑ j∈C(a,p h ) ωj = ∑ j∈C(b,p h ) ωj (5)

Learning an S-RMP model

Problem context

We are interested in learning S-RMP models in an indirect way. To account for that, we assume that the information provided by the DM takes form of pairwise comparisons of alternatives. Learning such a model amounts to setting values for the criteria weights and the reference points which are detailed as follow:

• The normalized weights of criteria, ωj, j ∈ M .

• The number of reference points, denoted by k.

• The reference points, p h , h ∈ P where p h = (p h 1 , p h 2 , ..., p h m ) is a vector in the evaluation space p h ∈ X .

• The lexicographic order of the reference points, defined as a permutation σ on the index set of the reference points P . For instance, p σ (1) is the first reference point to which alternatives are compared, ∃h ∈ P, σ(1) = h, etc.

To better understand the problem, we provide a brief review on the latest related researches.

Literature review

A first attempt to elicit indirectly the preference model from the holistic information given by the DM was presented in [START_REF] Mousseau | Inferring an electre tri model from assignment examples[END_REF] for the ELECTRE TRI method. They considered the pessimistic assignment rule, and developed a non-linear optimization formulation to infer all the parameters from a set of assignment examples. A similar approach was presented in [START_REF] Leroy | Learning the parameters of a multiple criteria sorting method[END_REF] for a simplified version of ELEC-TRE TRI named Majority Rule Sorting Model (MR-Sort). In [START_REF] Leroy | Learning the parameters of a multiple criteria sorting method[END_REF], only a few number of categories was considered and there was no veto threshold characterizing the discordance effect. Later, in [START_REF] Cailloux | Eliciting ELECTRE TRI category limits for a group of decision makers[END_REF], an extension of the previous work for ELECTRE TRI method has been presented in the group decision problems. The parameters were inferred from assignment examples, provided by multiple decision makers, based on a Mixed Integer Programming (MIP).

For the S-RMP method, [START_REF] Zheng | Preference elicitation for a ranking method based on multiple reference profiles[END_REF] was the first work that proposed an MIP to infer the parameters of this models from a given set of pairwise comparisons. They assumed the existence of an S-RMP model that is able to correctly restore all the pairwise comparisons given by the DM, and provided a linearized set of constraints with binary variables. However, the algorithm proposed in [START_REF] Zheng | Preference elicitation for a ranking method based on multiple reference profiles[END_REF] suffers from a limitation in the case where the pairwise comparisons provided by the DM contains inconsistent information. It cannot find any solution.

To overcome this problem, [START_REF] Liu | Preference elicitation from inconsistent pairwise comparisons for multi-criteria ranking with multiple reference points[END_REF] proposed a new algorithm dealing with inconsistencies in the elicitation process. The proposal is also based on an MIP with binary variables. However, the algorithm searches for an optimal solution that is compatible with as many as possible pairwise comparisons provided by the DM. Besides, it is able to identify the inconsistent information.

We remark a common issue for the two algorithms developed for S-RMP models that is the high computation time due to the introduction of the binary variables. Actually, when dealing with a very limited number of pairwise comparisons, the previous two algorithms are proved to be quite efficient. Nevertheless, the computation time increases exponentially when a large quantity of pairwise comparisons are provided.

Recently, [START_REF] Sobrie | Learning a majority rule model from large sets of assignment examples[END_REF] presented a well-adapted algorithm for learning the parameters of MR-Sort models from large datasets. It takes the advantages of a heuristic approach combined with a Linear Programming (LP). More details are provided in [START_REF] Sobrie | A majority rule sorting model to deal with monotone learning sets[END_REF]. We highlight the general idea of this work for sorting problems which in our work are ranking ones instead.

Preference elicitation for RMP models was also encouraged in [START_REF] Bouyssou | Multiattribute preference models with reference points[END_REF]. Besides, other works are generally interested in preference disaggregation. For interested readers, please refer to [START_REF] Doumpos | Preference disaggregation and statistical learning for multicriteria decision support: a review[END_REF][START_REF] Michael Doumpos | An evolutionary approach to construction of outranking models for multicriteria classification: The case of the ELEC-TRE TRI method[END_REF][START_REF] Doumpos | On the development of an outranking relation for ordinal classification problems: An experimental investigation of a new methodology[END_REF].

A metaheuristic approach

The different approaches presented in section 2.3.2 for learning S-RMP models have formulated the learning task as a Mixed Integer linear programming (MIP) optimization problem. However, a common observation is the considerable high computation time due to the introduction of binary variables [START_REF] Liu | Preference elicitation from inconsistent pairwise comparisons for multi-criteria ranking with multiple reference points[END_REF][START_REF] Zheng | Preference Elicitation for Aggregation Models based on Reference Points: Algorithms and Procedures[END_REF][START_REF] Zheng | Preference elicitation for a ranking method based on multiple reference profiles[END_REF].

In the case where the decision problem implies a large dataset, suchlike MIP consumes even more computation time and usually leads to the insolvability of the problem within the limited time. This type of problems is considered as hard optimization problems, which is defined as optimization problems that cannot be solved to optimality ( [START_REF] Boussaïd | A survey on optimization metaheuristics[END_REF]).

Thus, our proposal turns to a metaheuristic algorithm with the intention of dealing efficiently with large datasets. In contrast to the MIPs, the metaheuristic will approximate as accurately as possible the parameters of the S-RMP model. We remark that multiple S-RMP models may exist. Our objective is to infer a satisfactory S-RMP model that is compatible with most (if not all of them) of the input information. An interactive model calibration process in the form of supplementary constraints on the parameters should be invoked in real decision cases.

Overview of the algorithm

The proposed metaheuristic follows the general idea of Evolutionary Algorithms (EA) and makes use of a local optimization. It starts with an initialized population of N model S-RMP models rather than a single solution. Evaluation, Selection, Mutation and Substitution operations are iterated to adjust the parameters of each individual in the operating population. The iteration will repeats until either at least one of the individuals in the population is able to restore all the input pairwise comparisons or after a certain number of times. Actually, due to the limited computation time in practice, we interrupt the algorithm before it reach the optimality more often. In our context, the limited number of iterations is set to 100. It is summarized in Algorithm 1 Actually, on the one side, the reference points are adjusted while yielding mutants. On the other side, the criteria weights are adjusted by solving an LP. The rest of this section is structured as below: Section 3.2 begins with the initialization of the reference points and the criteria weights. Section 3.3 presents how to learn the reference points. The LP for learning the criteria weights is provided in Section 3.4.

Initialization

The first step of the algorithm consists in initializing a population of N model S-RMP models (in our context, N model = 10). In metaheuristic algorithms, it is important that the initial population spans adequately the solution space as shown in [START_REF] Kenneth | Genetic algorithms are not function optimizers[END_REF]. To do this, we provide a primitive method.

Reference points

First, we assume that the number of reference points k is fixed beforehand. From a practical perspective, the number of reference points in S-RMP models never exceeds 3, since such an S-RMP model has already a considerable capacity while restoring preference information as shown in [START_REF] Zheng | Preference Elicitation for Aggregation Models based on Reference Points: Algorithms and Procedures[END_REF].

Second, for a criterion j, k random numbers are generated from the uniform distribution in the evaluation scale Xj. Depending on the preference direction of the criterion j, they are then ranked either in ascending order such that p 1 j ⩽ p 2 j ⩽ ... ⩽ p k j or in descending order such that p k j ⩽ p k-1 j ⩽ ... ⩽ p 1 j . Thus, the generated reference points guarantee that p h ≻ p h-1 ≻ ... ≻ p 1 , where p h = (p h 1 , p h 2 , ..., p h j , ..., p h m ) Finally, for k reference points, there are k! possible lexicographic orders. We randomly choose one of them and it is fixed once initialized.

Criteria weights

Concerning the weights, we use the approach of Butler [START_REF] Butler | Simulation techniques for the sensitivity analysis of multi-criteria decision models[END_REF]. Firstly, m -1 random integer numbers di where i ∈ {1, ..., m -1} are generated uniformly between 0 and 100. Then, they are ranked such that d0 = 0 < d1 < d2 < ... < dm-1 < 100 = dm. The weight of a criterion j is defined as wj = (dj -dj-1)/100. It is accurate to 0.01 and 0 ⩽ wj ⩽ 1. This ensure that the weights sum up to 1 and are uniformly distributed.

Learning the reference points

Evolutionary algorithm

Our algorithm starts with an initial population of a certain number of random solutions, and evolves it while yielding a new generation of population at each iteration. Conventionally, the (t + 1)-th generation of solutions, denoted by Gt+1, is obtained from the n-th generation, denoted by Gt, through a procedure composed by several well-defined operations as described in [START_REF] William | An overview of evolutionary computation[END_REF]. [START_REF] Pirlot | General local search methods[END_REF] also proposed a review about Evolutionary Algorithms used in the context of local search.

In Algorithm 1, a proportion of the best S-RMP models in the population are selected at the beginning of each iteration. Then, they are submitted to the mutation operation. Mutations are performed on each selected model and yield mutants with a probability. The criteria weights of the mutants are adjusted (Section 3.4). After being re-evaluated, the mutants are considered as the "children", while the initial models are considered as the "parents". Finally, the substitution amounts to selecting the best models among both the parents and the children based on their evaluation. Otherwise, it also allows us to introduce newly initialized individuals to the next generation of population if necessary. The details are provided in the following sections.

Evaluation operation

The evaluation of the "best" individuals is done according to their fitness to the problem. In the context of learning an appropriate preference model from a set of information provided by the DM, several objective functions are applicable. In the case of learning S-RMP model, the fitness function corresponds to the Ranking Accuracy (RA) of the model, which is simply defined as the ratio of the number of pairwise comparisons restored correctly by the model to the total number of pairwise comparisons provided by the DM at the beginning of the process (( 6)). Actually, this is the most straightforward measure of ranking performance.

RA =

Number of pairwise comparisons restored Total number of pairwise comparisons provided (6)

Selection operation

The S-RMP model that gives the highest RA in the current population Gt will be selected randomly with replacement, with a probability ξ which increases with their fitness. Thereby, we define the probability ξ(fs) associated to a given S-RMP model fs that consists of a set of parameters s by:

ξ(fs) = RA(fs) -RAmin RAmax -RAmin (7)
where RAmin = min {RA(fs) | ∀fs ∈ Gt} (8)

RAmax = max {RA(fs) | ∀fs ∈ Gt} (9)

Mutation operation

Mutation operation, which amounts to adjust the reference points, is only applied to the selected models. For each selected S-RMP model, the operation is defined in the Algorithm 2. The mutations of the reference points proceed one by one based on their lexicographic order σ and then on each criterion from j = 1 to m. , are uniformly randomized. However, it should be bounded to prevent re-degrading the quality of the models after some iterations and the boundary should depend on the current number of iterations Nit. For example, we take herein:

θ σ(h) j ≤ ⌈ 50 √ Nit ⌉ (10)
As ±θ σ(h) j will be later applied to the current model with a probability µ, we denote the changed evaluation of the reference point p σ(h) on the criterion j by p ′σ(h) j where

p ′σ(h) j = p σ(h) j ± θ σ(h) j (11)
Then, we should define the probability µ that indicates if the variation ±θ σ(h) j should be applied. To do so, we identify at first the impacts provoked by ±θ σ(h) j on the judgment of preference between alternatives, which eventually improve (or worsen) the RA, through the calculations below.

For any pairwise comparison of alternatives a ≿ b provided by the DM, we define the slack quantity as follows:

s σ(h) (a,b) = m ∑ j=1 ( δ σ(h) a,j -δ σ(h) b,j ) • ωj ( 12
)
where ωj denotes the weight of the criterion j. ∀(a, b) ∈ BC, ∀h ∈ P and the lexicographic order σ, we compute δ

σ(h) a,j
and δ σ(h) b,j as follows:

δ σ(h) a,j = { 1 if aj ⩾ p σ(h) j 0 if aj < p σ(h) j , δ σ(h) b,j = { 1 if bj ⩾ p σ(h) j 0 if bj < p σ(h) j
(13) In fact, the (global) judgement of preference between (a, b) is based on the relative preference between (a, b) with respect to the reference point p σ(h) , which is deduced by

a ≿ σ(h) b ⇐⇒ s σ(h) (a,b) ≥ 0 (14) If δ σ(h) a,j
equals to 1, it means that the criterion j is contributing to the statement a ≿ σ(h) b. The criterion j is namely a contributing criterion for a in this case. Similarly, if δ σ(h) b,j equals to 1, it means that the criterion j is weakening the statement a ≿ σ(h) b. The criterion j is then namely a weakening criterion for b. Particularly, when δ σ(h) a,j -δ σ(h) b,j = 0, the criterion j is neither contributing nor weakening. We said that it is neutralized.

The variation ±θ σ(h) j

affects the calculation of s

σ(h) (a,b) through (δ σ(h) a,j , δ σ(h) b,j
). All the possible impacts provoked by ±θ σ(h) j can be summarized (as shown in Table 1) by the value change of (δ σ(h) a,j , δ σ(h) b,j ), where the changed value can be denoted and calculated by:

δ ′σ(h) a,j = { 1 if aj ⩾ p ′σ(h) j 0 if aj < p ′σ(h) j ( 15
)
δ ′σ(h) b,j = { 1 if bj ⩾ p ′σ(h) j 0 if bj < p ′σ(h) j ( 16
)
Formally, the positive and negative impacts are respectively defined as follows:

Ipos ( ±θ σ(h) j ) = { (a, b) ∈ BC (δ ′σ(h) a,j -δ ′σ(h) b,j ) -(δ σ(h) a,j -δ σ(h) b,j ) > 0 } (17)

Ineg

( ±θ (1,0) (0,0) negative contributing criterion neutralized I 6

σ(h) j ) = { (a, b) ∈ BC (δ ′σ(h) a,j -δ ′σ(h) b,j ) -(δ σ(h) a,j -δ σ(h) b,j ) < 0 } ( 18
)
(1,0) (1,1) negative contributing criterion neutralized I 7 (0,0) (0,1) negative criterion become weakening I 8

(1,1) (0,1) negative criterion become weakening I 9

(0,0) (1,1) neutral no impact on the slack I 10

(1,1) (0,0) neutral no impact on the slack I 11 (0,0) (0,0) neutral no impact on the slack I 12

(1,1) (1,1) neutral no impact on the slack I 13

(1,0) (1,0) neutral no impact on the slack I 14 (0,1) (0,1) neutral no impact on the slack

I 15 (0,1) (1,0) n.c. impossible I 16
(1,0) (0,1) n.c. impossible

For instance, we represent, in Figure 1, a couple of examples to illustrate the adjustment of reference points. We suppose that • The criteria (c1, c2, c3 and c4) are equal weighted, i.e. ∀j ∈ {1, 2, 3, 4}, ωj = 0.25. • The reference points are uniformly initialized and used in sequence, i.e. first p 1 , then p 2 and p 3 . • The DM provided a ≿ b as input.

With the initialized reference points, we derive at first that a ∼ p 1 b, b ≿ p 2 a and a ≿ p 3 b. Since the reference points are used sequentially, we derive then b ≿ a, which is inconsistent with the statement of the DM. Thus, the reference points should be adjusted by the algorithm.

In Figure 1(a), the criterion j = 3, which was originally weakening the statement a ≿ b (as δ p 2 a,3 = 0), is neutralized by applying the variation -θ p 2 3 to the reference point p 2 on the criterion j = 3, because δ ′p 2 a,3 = δ ′p 2 b,3 = 1 and δ ′p 2 a,3 -δ ′p 2 b,3 = 0. It is considered to be a (potential) positive impact, for the fact that, as shown in this example, we newly derive that a ∼ p 2 b and finally a ≿ b, which is consistent with the DM's statement. In this sense, the reference point p 2 is adjusted.

In Figure 1(b), the criterion j = 2, which was originally neither weakening nor contributing to statement a ≿ b (as δ p 2 a,2 -δ p 2 b,2 = 0), becomes contributing by applying the variation +θ p 2 2 to the reference point p 2 on the criterion j = 2, because δ ′p 2 a,2 = 1 and δ ′p 2 b,2 = 0. It is also considered to be a (potential) positive impact, for the similar fact that we can again finally derive that a ≿ b in this example with the adjusted reference point p 2 even though it was adjusted in a different way. So far, we identify four positive possibilities as well as four negative possibilities and sum up accordingly the quantity of the positive (respectively, negative) impacts by

Ipos = Ipos ( ±θ σ(h) j ) = |I1| + |I2| + |I3| + |I4| ( 19
)
Ineg = Ineg ( ±θ σ(h) j ) = |I5| + |I6| + |I7| + |I8| (20)
However, we need to make some empirical remarks such that:

1. Usually, the variation ±θ σ(h) j

provokes both the positive and negative impacts on different pairs of alternatives at the same time. Having more positive impacts means that we have more chance to improve the S-RMP model quality, but not necessarily. 2. Sometimes, we provoke as many negative impacts as the positive impacts. In such a case, we are more interested in the case where a large quantity of impacts provoked, as it allows us to introduce more diversity to the population. 3. When there are less positive impacts than the negative ones, or even no positive impacts provoked, we consider that it is still possible to improve the model potentially.

Based on these observations, we define accordingly the probability of applying the variation ±θ ) as it is associated to the variation ±θ σ(h) j , as below :

µ ( ±θ σ(h) j ) =      Ipos Ipos+Ineg if Ipos ̸ = 0 and Ipos > Ineg 1 -1

Ipos

if Ipos ̸ = 0 and Ipos = Ineg 0.5 * Gaussian if Ipos = 0 or Ipos < Ineg (21) Otherwise, as we treat the reference points one by one according to the lexicographic order, once ( 14) is affirmed, (a, b) will be appended to the list for keeping ignorable pairwise comparisons. Actually, (a, b) is considered to be ignorable, because the adjustment for the rest of the reference points will not affect the global judgement on the preference between (a, b) and the impacts provoked by ±θ

σ(h+1) j on (δ σ(h+1) a,j , δ σ(h+1) b,j
) (and so on, if exists) can be then ignored. In other words, once the comparison between any two alternatives (a, b) is restored by the previous reference point, no matter how they compare to the rest of the reference points, it will not be changed.

Substitution operation

As presented in Section 3.3.1, the mutants with the adjusted criteria weights (Section 3.4 for details) will be re-evaluated and then appended to the n-th generation of population Gt. The population that contains both the "children" and the "parents" is denoted by G ′ t . Substitution operation decides who will survive in the (t + 1)-th generation of population.

Firstly, it is about selecting the best individuals from the n-th generation of population G ′ t . Therefore, the principle is the same as the selection operation presented in Section 3.3.3. We define the probability ζ associated to an S-RMP model fs of being submitted to the (t + 1)-th generation by:

ζ(fs) = RA(fs) -RA ′ min RA ′ max -RA ′ min ( 22
)
where,

RA ′ min = min { RA(fs) | ∀fs ∈ G ′ t } ( 23
)
RA ′ max = max { RA(fs) | ∀fs ∈ G ′ t } (24)
Secondly, if there is not enough number of individuals that are good enough to be submitted, some newly initialized S-RMP models will be appended to the (t + 1)-th generation under the constraint of the initial size of the operating population.

After having defined the substitution operation, the evolutionary part of the algorithm for learning the reference points is completed. Now, we discuss how to adjust the criteria weights.

Learning the criteria weights

Assuming that the reference points are given, we consider the adjustment of criteria weights as a linear optimization problem.

We are interested in setting good values for the criteria weights ωj, j ∈ M to restore as many as possible pairwise comparisons provided by the DM. However, to do that we usually need to introduce binary variables into the program as it is described in [START_REF] Zheng | Preference Elicitation for Aggregation Models based on Reference Points: Algorithms and Procedures[END_REF]. In the case of a large dataset, using binary variables will considerably increase the computation time ( [START_REF] Zheng | Preference elicitation for a ranking method based on multiple reference profiles[END_REF][START_REF] Liu | Preference elicitation from inconsistent pairwise comparisons for multi-criteria ranking with multiple reference points[END_REF]). Therefore, we propose in this paper an alternative formulation without binary variables that has been proved to be quite efficient (Section 4 for the numerical test results).

We note at first that ωj, j ∈ M are the unknown variables to be adjusted. Considering then the slack defined by ( 12) and ( 13), we remark that it is actually an indicator that shows if the (relative) preference between two alternatives can be restored by the current S-RMP model. Hopefully, for any pairwise comparison (a, b) provided by the DM, there should be at least one of the slack variables s

σ(h) (a,b) such that ∃h * ∈ P, s σ(h * ) (a,b) > 0, s σ(1) (a,b) = s σ(2) (a,b) = ... = s σ(h * -1) (a,b) = 0 (25)
where σ is a specified lexicographic order of reference points.

However, we cannot translate (25) as part of the linear constraints of the problem, because it is not always true. It means that it is not always possible to find such an S-RMP model with the given reference points and the adjusted criteria weights that restores exactly all the provided pairwise comparisons.

In fact, our objective is, by adjusting the criteria weights, to make the model compatible with as many as possible pairwise comparisons provided by the DM. In this case, (25) should be integrated as part of the objective function but not as linear constraints.

Hence, to account for this, we define then another two auxiliary variables s

σ(h)+ (a,b) and s σ(h)- (a,b) for ∀h ∈ P, ∀(a, b) ∈ BC by: s σ(h) (a,b) -s σ(h)+ (a,b) + s σ(h)-
∑ h=1 ω σ(h) • ( s σ(h)+ (a,b) -α • s σ(h)- (a,b) ) ( 27
)
where we make use of two weighting system α and ω σ(h) . On the one hand, ω σ(h) for ∀h ∈ P weights the reference points. We remark that, by carefully choosing their values, it is approximately equivalent to using the reference points in their lexicographic order. Actually, if k ⩾ 2, whatever the value of ω σ (k) ,

ω σ(k-h) = 1 ϵ • k ∑ i=k-h+1 ω σ(i) (28)
that depends on the small value for ϵ. ϵ represent the accuracy of the criteria weights that we impose in the solution. For example, if ϵ = 10 -3 and we take ω σ(3) = 1 for an S-RMP model with 3 reference points, then we can deduce that ω σ(2) = 10 3 and ω σ(1) = 10 6 +10 3 . On the other hand, α balances the compensatory behavior between pairwise comparisons. We note that, to avoid using binary variables, such an objective function cannot garantee that the maximal number of pairwise comparisons can be correctly restored. However, by setting a big enough value for α (for example, α = 10 3 ), it is able to reduce the compensatory effects efficiently.

We summarize the linear program as below:

max ∑ (a,b)∈BC k ∑ h=1 ω σ(h) • ( s σ(h)+ (a,b) -α • s σ(h)- (a,b) ) s.t. m ∑ j=1 ωj = 1 ∀(a, b) ∈ BC, ∀h ∈ P, s σ(h) (a,b) = m ∑ j=1 ( δ σ(h) a,j -δ σ(h) b,j ) • ωj s σ(h) (a,b) -s σ(h)+ (a,b) + s σ(h)- (a,b) ≥ 0 s σ(h)+ (a,b) ≥ 0, s σ(h)- (a,b) ≥ 0

Numerical analysis

In this section, the proposed algorithm is further investigated numerically based on a large set of artificially generated data. Statistical techniques were used to demonstrate the advantages of the metaheuristic comparing with the conventional disaggregation methods for outranking models based on LP.

We are firstly concerned with the quality of the solution. For a suchlike evolutionary approach, it is represented by the "best" individual in the operating population. Secondly, there is usually a compromise between the optimality of the solution and the number of iterations it takes to reach a satisfactory result as we mentioned in Section 3.1. By simulating the different decision circumstances with a large quantity of randomly generated data, we investigate the improvement curves produced in each circumstance to better understand the behavior of the algorithm. Otherwise, we examine also the runtime characteristics of the proposed metaheuristic, since one of the interest of this work is to overcome the insolvability of the MIPbased disaggregation methods in dealing with large datasets (Section 2.3.2).

The different decision circumstances that we studied are discussed in the next section.

Experiments

In the experiments, we consider 1000 alternatives and a varying numbers of criteria (4, 6, 8 or 10). We generate randomly an initial S-RMP model to test the algorithm. It is denoted by f β , as it is defined by a set of parameters β. It simulates the preference model of a fictitious DM. Then, 500 reference pairwise comparisons among the alternatives are derived from the initial model f β based on their randomly generated evaluations on the criteria. We consider also different numbers of reference points in the initial model f β (Ini. NRP, for Initial Number of Reference Points) to simulate different complexity levels of the DM's preference system.

The experiments are divided into two groups according to the different preset numbers of reference points in the inferred model (Inf. NRP): In Group A, we infer S-RMP models with only 1 reference point whatever the initial number of reference points in β. While, in Group B, we infer S-RMP models with exactly their initial number of reference points, so as to compare these two groups of experiment and better understand the importance of carefully setting a corresponding number of reference points in S-RMP disaggregation. The different cases that we considered are distinguished by:

• Number of criteria = 4, 6, 8 or 10 • Ini. NRP = 1,2 or 3 • Inf. NRP = 1

or the initial number

The experiments are then numbered by a trio-index as shown in Table 2 and Table 3. For instance, "4.1.A" means that, we consider 4 criteria and the initial model is generated with 1 reference point and it is in the group "A". For each experiment, all the provided results are based on the average of 100 repeated trials under the same testing conditions2 .

Empirical results

Quality of the solution

The inferred model is denoted by f β * with some parameters β * . As shown in Table 2 andTable 3, the "closeness" between the initial model f β and the inferred model f β * is measured by the RA. Since we generated the input information without introducing any inconsistency, the measured value of RA is always expected to be as close as possible to 1.00. The "Starting RA" represents the quality of the "best" individual in the initialized population. The "Final RA" represents the quality of the final inferred model f β * , which is the "best" individual in the population when the algorithm terminated. However, we remind that multiple non-identical individuals may give the same value of RA, as the solution is not unique (as presented in Section 2.3.1).

In the group A of experiments, we generated the initial models with 1, 2 or 3 reference points and inferred S-RMP models involving only 1 reference point (Table 2). We notice that the best final RA, 0.997, is observed in the experiment 4.1.A, while the worst is observed in the experiment 4.3.A. It shows, in general, the capability of the proposed metaheuristic while inferring an S-RMP model as simply as possible.

However, exceptions are observed in the experiment 10.3.A. The final RA is instead even greater than 10.1.A on average. The randomized starting models are also better than in the other two cases (10.1.A and 10.2.A). This will be further discussed (Section 4.3). In the group B of experiments, we assume that we can correctly set the number of reference points in the S-RMP models to be inferred (Table 3). At this time, the worst final RA is observed in the experiment 10.3.B. It is reasonable, since it is the most complicated case that has been considered in our analysis. However, 0.950, which means 95% of the 500 input pairwise comparisons could be restored correctly in the inferred S-RMP model, is already satisfactory enough for the disaggregation. By comparing with the group A, we observe that the starting RAs as well as the final RAs are significantly enhanced in the group B of experiments by presetting correctly the number of reference points, especially in the experiments where the initial models are generated with more than one reference points (except for 10.3.A and 10.3.B).

To better understand the behavior of the proposed algorithm, the improvement curves are drawn and investigated below.

Behavior of the algorithm

The improvement curves show not only the starting and the final point, but the evolution of RA in function of the number of iterations during the whole process. As shown in Figure 2 (for Group A) and Figure 3 (for Group B), the RA is improved and converges progressively. Without any surprise, the exceptional phenomenon that we pointed out in Section 4.2.1 for 10.3.A is also observed in Figure 2(d). The dotted curve that represents the case involving 3 reference points is instead placed above the curve that represents the case involving 1 reference point (Section 4.3 for further discussions).

For a fixed number of criteria, by comparing each of the subfigures in Group A with the ones in Group B (for example, Figure 2(a) with Figure 3(a)), we can firstly visualize the intensively reduced level of RA due to the improperly set number of reference points. Moreover, we notice that the reduced level is more significant when the model involves fewer criteria.

Otherwise, we observe that, in Figure 3, the dashed line for "2p" and the dotted line for "3p" almost coincide with each other. It means that, comparing with the single reference point cases, the improving RA is reduced quasi-equally in the multiple reference points cases regardless of the number of reference points involved.

Actually, the improvement curves can also be grouped by number of reference points (as shown in Figure 4 and Figure 5) instead of by criteria number. In Figure 4(c), we inferred S-RMP models that involve only 1 reference point from the information derived from the initial models that involve 3 reference points. We notice that, not similar to the other figures, the curves (in Figure 4(c)) that represent the cases where more criteria are involved are above the others that represent the cases where fewer criteria are involved. Related discussions are provided in Section 4.3.

Moreover, the improvement curves can be divided into two stages. The first stage is the rising section, while the second stage is the flat section. The flat section shows the final RA that we can reach, (d) 10 criteria while the rising section also reveals some important advantages of the proposed metaheuristic.

We measure the first stage of the improvement curves by the rapidity of convergence. It is defined as the necessary number of iterations needed to reach a quasi-satisfactory RA. To be precise, we define the quasi-satisfactory RA as the 97% of the final RA. It shows the efficiency of the metaheuristic algorithm. The results are presented in Table 4.

We can observe that the metaheuristic approach is generally quite efficient. Table 4 shows that the inferred models can reach a quasisatisfactory state within around the first 30 iterations and it depends on the complexity of the preference model of the DM.

Computation time

The average computation time is provided in Table 5. We observed that, even for the case involving 3 reference points and 10 criteria, the computation time remains within around 3 minutes. These results should be further compared to the computation times of algorithms using MIP formulations.

Besides, not only should the runtime value be measured, but also the evolution of computation time as a function of the criteria number or the number of the reference points should be investigated. It is interesting while comparing with other MIP-based algorithms, as it demonstrates the significant runtime advantage of the metaheuristic algorithm. By drawing the runtime curves, as shown in Figure 6(a), we notice that the computation time is proportional with the criteria number. We remind that it increases exponentially in the MIP-based S-RMP disaggregation methods (Section 2.3.2). It endues the metaheuristic approach with the capacity of dealing with the decision instances that involve a large quantity of evaluating criteria.

Besides, we also notice that we consume more computation time when dealing with more reference points, but not as much as for the LP-based methods. To be clear, we draw once again the runtime curves in function of the number of reference points (as shown in Figure 6(b)). At this point of view, the obtained curves are all sublinear regardless of the criteria number while in the LP-based S-RMP disaggregation methods, they are usually exponential vs. the number of reference points.

Discussion

By checking the result of the numerical tests, we notice that the preset number of reference points should be carefully adjusted to derive a better solution, since the final RA is intensively reduced (except for 10.3.A) in Group A (i.e. learning S-RMP models with single reference point) of experiments when the initial number of reference points increases. Nevertheless, starting from an S-RMP model with single reference point shows greater interpretability to the problem. It actually expresses a very natural rule based on the distinction of two classes of evaluation on each criterion. Besides, it helps us to estimate the complexity of S-RMP models and to fix the number of reference points. It means that we should consider increasing the number of reference points when the final RA is not as satisfactory as we expected.

Moreover, we also notice that the number of criteria should be limited. On one hand, we observe that, when the number of reference points is fixed, the final RA is descending progressively when the criteria number increases. On the other hand, the exceptional phenomenon that we observed (10.3.A) shows that, when there are more criteria, it is easier to derive an S-RMP model that involves only one reference point and gives a better RA. In other words, it is easier to adjust one reference point by calibrating its valuation and the weights of criteria than to calibrating multiple reference points at the same time.

The study about the rapidity of convergence of the algorithm could help us to adjust when to interrupt the algorithm when dealing with real datasets, and to economize the total computation time in different application circumstances. For instance, in the online web ap- plications, the response time should usually be much more valuable than in the offline recommendations. To obtain the most accurate result, it is possible to run the program during a very long time, even day and night, in the offline cases. However, an acceptable solution should be worked out as quickly as possible in the online cases and then be adjusted step by step by a follow-up interactive process.

Conclusion

This paper presents an efficient metaheuristic approach to infer S-RMP models from a large set of pairwise comparisons provided by the DM. The proposed algorithm was tested with a large quantity of artificially generated data that simulates a variety of different decision circumstances. Firstly, the metaheuristic is able to deal with instances involving as many as 500 pairwise comparisons. Suchlike instances cannot be solved using MIP formulations. Secondly, the computation time is proportional with the number of criteria involved and sub-linearly increases with the number of reference points. Finally, even if the metaheuristic is not able to learn an S-RMP model which to restore all pairwise comparisons, it infers S-RMP models which restore up to at least 95% of the input information within a reasonable computation time. We remark that it constitutes a good trade-off of quality of result vs. computation time.

One of the interesting questions emerged from this work is the adjustment of the number of reference points. Moreover, the performance of the proposed metaheuristic should also be explored when we are in presence of inconsistencies, for example, in the case of group decision problems or in the case of real world applications.

Algorithm 1 : 4 5 6 Adjust the criteria weights for mutated individuals 7 Evaluate the newly adjusted individuals 8

 145678 OVERVIEW OF THE ALGORITHM 1 Initialize a population of N model S-RMP models 2 Evaluate each individual in the population 3 repeat Select the best individuals from the population with a probability ξ Adjust the reference points for selected individuals and yield mutants with a probability µ Substitute the population with a probability ζ under the constraint of the initial population size 9 until at least one of the stopping criteria is met;

Algorithm 2 : 3 for each criterion j do 4 Generate a random variation ±θ σ(h) j 5 Count possible impacts caused by this change 6 Apply the change to the current model with a probability µ 7 Append newly ignorable pairwise comparisons to the list 8

 2345678 MUTATION OPERATION1 Create an empty list for keeping ignorable pairwise comparisons 2 for each reference point p σ(h) do Yield a mutated model Firstly, both the sign and the value of the variation, denoted by ±θ σ(h) j

Figure 1 :

 1 Figure 1: Adjustment of reference points (a) Impact of -θ p 2 3 (i.e. I2) (b) Impact of +θ p 2 2 (i.e. I4)

  σ(h) jto the current S-RMP model, denoted by µ(±θ σ(h) j

  them are positive by definition. Actually, s σ(h)+ (a,b) represents the positive terms in the definition of s h (a,b) that contributes to the statement. Respectively, s σ(h)-(a,b) represents the negative terms in the definition of s h (a,b) that weaken the statement. Moreover, in order to maximize the number of pairwise comparisons correctly restored by the model, the objective function is defined as below: max ∑ (a,b)∈BC k

Figure 2 :

 2 Figure 2: Group A of experiments (num. of criteria)

Figure 3 :

 3 Figure 3: Group B of experiments (num. of criteria)

Figure 4 :

 4 Figure 4: Group A of experiments (num. of ref. pts)

Figure 5 : 1 .

 51 Figure 5: Group B of experiments (num. of ref. pts)

Figure 6 :

 6 Figure 6: Computation time (in seconds) (a)

Table 1 :

 1 Impacts provoked by ±θ

					σ(h) j	on	( δ	σ(h) a,j , δ b,j σ(h)	)
		Before After	Impact	Description
	I 1	(0,1)	(0,0)	positive	weakening criterion neutralized
	I 2	(0,1)	(1,1)	positive	weakening criterion neutralized
	I 3	(0,0)	(1,0)	positive	criterion become contributing
	I 4	(1,1)	(1,0)	positive	criterion become contributing
	I 5

Table 2 :

 2 Group A of experiments

	Exp.	Num. Cri. Ini. NRP Inf. NRP Starting RA Final RA
	4.1.A	4	1	1	0.807	0.997
	4.2.A	4	2	1	0.746	0.904
	4.3.A	4	3	1	0.735	0.880
	6.1.A	6	1	1	0.813	0.987
	6.2.A	6	2	1	0.761	0.934
	6.3.A	6	3	1	0.741	0.912
	8.1.A	8	1	1	0.810	0.969
	8.2.A	8	2	1	0.760	0.940
	8.3.A	8	3	1	0.746	0.925
	10.1.A	10	1	1	0.805	0.960
	10.2.A	10	2	1	0.761	0.931
	10.3.A	10	3	1	0.810	0.969

Table 3 :

 3 Group B of experiments

	Exp.	Num. Cri. Ini. NRP Inf. NRP Starting RA Final RA
	4.1.B	4	1	1	0.807	0.997
	4.2.B	4	2	2	0.810	0.983
	4.3.B	4	3	3	0.819	0.973
	6.1.B	6	1	1	0.813	0.987
	6.2.B	6	2	2	0.798	0.966
	6.3.B	6	3	3	0.796	0.968
	8.1.B	8	1	1	0.810	0.969
	8.2.B	8	2	2	0.780	0.961
	8.3.B	8	3	3	0.784	0.961
	10.1.B	10	1	1	0.805	0.964
	10.2.B	10	2	2	0.777	0.952
	10.3.B	10	3	3	0.775	0.950

Table 4 :

 4 Rapidity of convergence (group B, in number of iterations)

	Num. Cri. 1 ref. pt 2 ref. pts 3 ref. pts
	4	14	21	19
	6	20	21	22
	8	20	29	31
	10	23	29	32

Table 5 :

 5 Computation time (group B, in seconds)

	Num. Cri. 1 ref. pt 2 ref. pts 3 ref. pts
	4	44.55	69.58	77.98
	6	85.78	101.48	110.20
	8	119.67	136.46	146.35
	10	148.28	169.75	182.06

Laboratoire Génie Industriel, Ecole Centrale Paris, Grande Voie des Vignes, 92295 Châtenay-Malabry, e-mail: jinyan.liu@ecp.fr, wassila.ouerdane@ecp.fr, vincent.mousseau@ecp.fr

Intel Core i3-2120

3.30 GHz CPU, 4 GB RAM, Ubuntu 14.04 LTS, Eclipse IDE Kepler SR2, Cplex 12.6


Manel Mammar

Vincent Mousseau
email: vincent.mousseau@ecp.fr

Wassila Ouerdane
email: wassila.ouerdane@ecp.fr

Alexandre Aubry
email: alexandre.aubry@placedesleads.com

Manel Maamar
email: manel.maamar@ecp.fr

INTERNET PROSPECTS' FLOWS FORECASTING FOR A MULTI-PERIOD OPTIMIZATION MODEL OF OFFER/DEMAND ASSIGNMENT PROBLEM

Keywords: Assignment problem, multi-period optimization, Internet prospects forecasting Average daily turnover Forecasting error

published or not. The documents may come

INTRODUCTION

Various methods of forecasting have been proposed in the literature. The aim of such a methods is to reduce the uncertainty related to the unknowing of the future. Several approaches of forecasting, namely endogenous, exogenous and hybrid exist (see [START_REF] Box | Time Series Analysis: Forecasting and control, 4 th Edition[END_REF]). Most of these methods have been implemented and regrouped in R software [START_REF] De Micheaux | Le logiciel R, 2 nd Edition[END_REF]. On the other hand, a well-known problem is the standard Assignment Problem (AP) where the aim is to optimally resolve the problem of assigning n jobs (tasks) to n individuals (machines), such that minimum cost or maximum profit can be obtained [START_REF] Burkard | Assignment Problems[END_REF]. The classical AP has a wide range of applications that span manufacturing, personnel scheduling and transportation, to name a few [START_REF] Dessouky | Production scheduling on single-stage multiproduct batch chemical process with fixed batch sized[END_REF][START_REF] Soumis | A model for large-scale aircraft routing and scheduling problems[END_REF]. Various algorithms were developed to find solutions efficiently [START_REF] Balinski | A competitive (dual) simplex method for the assignment problem[END_REF][START_REF] Kuhn | The Hungarian method for the assignment problem[END_REF]. For a standard AP, only the cost or profit from each possible assignment is considered in the problem formulation. However, in real applications other issues can be addressed.

In this paper, we study a variant of the AP that involves assigning Internet prospects to orders over several time periods using two objectives to guide the search for solutions. The AP concerns Place des Leads (PdL), a marketplace where several clients (assurances, banks...) make orders of Internet prospects. These later are defined in a Web context as a contact ready to carry out an action (purchase, estimate...) [START_REF] Maamar | Multi-period modeling and optimization of offer/demand assignment problem: a marketplace application[END_REF]. More precisely, we aim to combine a forecasting model with an optimization model in order to enhance the PdL's assignment problem results. The majority of the research works deals with either forecasting models or optimization models. In this work, our idea is to integrate forecasting model into a multi-period assignment in order to improve assignment by integrating knowledge about Internet prospects' flows continuously evolving in the optimization model.

To account for that, we introduce, in section 2, the multi-period assignment model and its characteristics. In section 3, (i) we discuss the relationship between both forecasting model and optimization one, (ii) we propose a forecasting system which considers all characteristics of the multi-period optimization model and (iii) we show how this forecasting is integrated in the multi-period optimization model. To show the efficiency of the proposed model, numerical tests are addressed in section 4. A conclusion and some perspectives are discussed in section 5.

MULTI-PERIOD OPTIMIZATION MODEL OF OFFER/DEMAND ASSIGNMENT

In this section, we introduce a multi-period assignment model, defined in [START_REF] Maamar | Multi-period modeling and optimization of offer/demand assignment problem: a marketplace application[END_REF], where the aim is to establish an adequacy of offer and demand with consideration of several periods. In particular, the objective is to optimize turnover while: (i) considering instantly new offer and demand of Internet prospects and (ii) taking into account the eventual future system updates. The idea is to optimize assignment, at each time step, by leveraging knowledge about evolution of Internet prospects' flows as well as updating of orders' quotas, over several future periods.

The multi-period assignment problem is characterized by moving horizon, which means that time is modeled through several time steps. At each time step t, we consider a set T of periods such that T = {t': t' ∈ {1, .., g}}, where t' = 1 is the current period and t' ∈ {2, .., g} are future periods. Thus, at each time step t, we have to decide which Internet prospects will be assigned to which orders, while considering all time periods of the set T. We assume that both orders and clients are the same over all periods t'. However, it is possible to have a new Internet prospects' flows arriving in an indeterminist way over all future periods t'. Consequently, at each period t', we have m t' Internet prospects and n orders coming from k clients that can make one or more orders, represented respectively by the following sets: At each period t', Internet prospect's assignment to order is represented by the following binary variable: { Moreover, the set of all existing "order/Internet prospect" assignments at period t' is described by the set X t' = {x ijt' : i ∈ O, j ∈ P, t' ∈ T}. However, in each X t' , some pairs x ijt' are not possible. Thus, we will consider that x ijt' ∈ X t' only if an Internet prospect j corresponds to an order i, at a period t'. Moreover, some orders ask for Internet prospects in exclusivity. It means that Internet prospect is sold at most once. Otherwise, we can sell Internet prospect several times to different orders. We call this sale: non exclusivity. Accordingly, we have E t' : set of x ijt' , at period t', such that i is an exclusive order and Ē t' : set of x ijt' , at period t', such that i is a non exclusive one. Thus, X t ' = {E t' ∪ Ē t' } such that t' ∈ T.

Internet prospects' assignment problem while considering, at each time step t, g periods t' is submitted to the followings constraints.

 Constraint (1) (Constraint ( 2)) makes sure that the number of Internet prospects assigned to each order i (client k respectively), at each period t', respects remaining volume of the order, which is the difference between remaining volume of order i (client k respectively) at period 1, and delivered Internet prospects' number to order i (client k respectively) from period 1 to happening period t'. Such that:  Ro i (Rc k ) is Internet prospects' number remaining to be delivered to order i (client k respectively), at period 1, which is defined as:

Ro i = min {Go i , Ho i , Do i , Mo i } (Rc k = min {Hc k , Dc k , Mc k } respectively), where: Go i is global volume of order i; Ho i (Hc k ) is hourly quota of order i (client k respectively); Do i (Dc k ) is daily quota of order i (client k respectively) and Mo i (Mc k
) is monthly quota of order i (client k respectively).  a ik is a binary variable described as follow: {  Constraint (3) expresses the fact that, over all g periods considered in the optimization, each exclusive Internet prospect can be sold no more than once and each non exclusive one can be sold no more than b times. Such that: S j is number of previous sales before period 1 of Internet prospect j, and b is sales' number of non exclusive Internet prospects.  Constraint (4) guarantees that each client k receives each Internet prospect no more than once over all g periods considered in the optimization.

Finally, function f represents the multi-period model's objective function which seeks to maximize turnover, over all periods considered in the optimization. Such that: P i is sell price of order i.

INTERNET PROSPECTS' FLOWS FORECASTING

Why should we consider forecasts?

We are in a context where Internet prospects' flows evolve continuously in a determinist way over time. In order to optimize the turnover, we have to take into account this Internet prospects' flow in the assignment process. In other words, multi-period assignment model defined in the previous section allows optimizing assignment by leveraging knowledge about evolution of Internet prospects' flows. In the sense that, we provide to our model the set of all possible assignments of each period t' considered in the optimization. Consequently, we have to compute Internet prospects' flows arriving over these periods.

Defining a forecasting model of Internet prospects' flows is therefore necessary to improve optimization model. Therefore, we propose in what follows to articulate a forecasting model with an optimization model.

Forecasting/Optimization articulation

Building a forecasting model requires to define which data to forecast before considering any forecasting method. In other words, we have to choose the appropriate historical data, time step and horizon which are the parameters of the forecasting model. The challenge is to set the parameters that will fit with the functioning of the optimization.

Data to forecast

In order to integrate Internet prospects' flows in the optimization model, we have to forecast at each period t' considered in the optimization model: (i) the number of Internet prospects arriving at this period and (ii) orders for which these Internet prospects can be assigned. A first approach could consist in forecasting the number of Internet prospects for each order. However, this idea is unsatisfactory because a large part of Internet prospects is shared by several orders. Consequently, if we forecast Internet prospects' number for each order, then Internet prospects shared by several orders will appear in the historical data of each of these orders. Thus, they will be counted several times. In addition, we need to know for each foreseen Internet prospect, the set of orders whose characteristics are satisfied. However, forecasting of Internet prospects' number for each order does not provide information about correspondence of these foreseen Internet prospects to the other orders.

To overcome the previous problems, we propose to forecast Internet prospects' number for each orders' group sharing the same Internet prospects. Accordingly, we will forecast Internet prospects' number as well as orders for which each foreseen Internet prospect can be assigned. Thus, each group of orders represents a time series for which a forecasting study will be conducted. How to build such order's group is discussed in section 3.3.1.

Forecasting time step and horizon

We face a complex forecasting problem, where the challenge is to integrate the forecasting results into a multi-period assignment model. The first difficulty comes from the fact that the forecasting of Internet prospects should be compatible with the multi-period optimization problem described in section 2. The second difficulty is related to the fact that the assignment of Internet prospects to orders is grounded, at each time step, on a rolling horizon of the future periods. Consequently, we have to compute Internet prospects' flows over these future periods such that, on the one hand, forecasting time step must be equal to the period of the multi-period assignment model, on the other hand, horizon should represent the number of the future periods considered in the multiperiod assignment model.

The fact that forecasting time step and period of the optimization model (t') must be the same generates another difficulty. In fact, the orders' quotas are updated at each hour. Therefore, integrating such information in the optimization model requires considering periods t' of one hour [START_REF] Maamar | Multi-period modeling and optimization of offer/demand assignment problem: a marketplace application[END_REF]. However, one hour forecasting time step does not necessarily guarantee a good forecasting because Internet prospects' flows density is low on this period. Therefore, we should consider periods t' greater than one hour in order to improve forecasting. The consequences of this, only some hourly updating of orders' quotas are integrated in optimization model. For example, considering periods t' of 4 hours conducts to take into account only updating at 00 o'clock, 4 o'clock, 8 o'clock …. although we have a quotas' updating at each hour. However, forecasting with 4 hours time step will be better than forecasting with one hour time step because Internet prospects' density is more important with 4 hours time step. Accordingly, we have to choose a compromise solution.

Forecasting system

In the previous sections, we highlighted both the benefits and the difficulties of integrating forecasting into the optimization model. In this section, we present: (i) an algorithm which computes the time series and (ii) a forecasting model as well as its integration in the multi-period optimization model.

Time series' construction

The following algorithm computes the set of time series. More precisely, each time series, which is an orders' group sharing the same Internet prospects, provides the Internet prospects' number of this orders' group at each historical period.

Forecasting model and its integration in the multi-period optimization model

In order to compute the forecasting of each time series (y lt ) t∈T such that l ∈ δ, defined in the previous section, we use ETS forecasting model implemented in R which takes as input a time series and returns the appropriate smoothing method based on the time series. In our case, ETS forecasting model returns, at each period, the number of Internet prospects (IP) shared by the same orders' group.

To integrate Knowledge on new Internet prospects' flows in the g periods considered in the optimization model, we need both number of Internet prospects at each period as well as orders for which these Internet prospects can be assigned. For instance, if we consider the following example: first, the assignment is done on 3 periods, namely 0, 1 and 2. Second, we have three orders and two Internet prospects at period 0, such that the possible assignments are show in Table 1. Finally, the results of ETS forecasting are depicted in Table 2. Given these information, the question is how to build the set of all possible for period 1 and 2. We may integrate, at period 1, one IP for order 1 (see table 3) and two IPs, at period 2, for orders 2 and 3 (see table 4).

Table 1: Set of all possible assignments at period 0 Table 2: ETS forecasting

Table 3: Set of all possible assignments at period 1 Table 4: Set of all possible assignments at period 2

EXPERIMENTATION

In this section, we analyze the behavior of the proposed forecasting system. We perform a numerical investigation using real data from PdL database. Our aim is to check empirical validity of the proposed approach and validate its practical usefulness. Experiment is designed to address the ability of the forecasting to optimize the results of the multi-period assignment model in terms of turnover. To do so, we perform assignment using both multi-period assignment model which includes an exact forecasting, i.e. the current data; we call this OF system and multi-period assignment model without forecasting, called O system. We compute the turnover for both systems OF and O, and evaluate the benefit induced by the introduction of an error free forecasting.

Figure 1: Difference between average turnovers of both OF system and O system over one day

Figure 1 represents the evolution, over one day, of the difference between the average turnovers of both OF system and O system. At the beginning of the day, the OF system generates an average turnover less than O system's average turnover. This situation continues until 6 o'clock where the OF system begins to be better than the O system. At 10 o'clock, the OF system retrieves all lost sales compared to O system and continues to be better until 10h30. After that, it begins to lose compared to O system. At the end of the day, OF system generates turnover more important than O system. This observation is not surprising because O system makes, at each current period, all possible assignments. However, the OF system may do not satisfy some orders, at the current period, in order to perform sales more expensive later. In fact, the system has information about the evolution of Internet prospects' flows over the future periods. We see this phenomenon at the beginning of the day as well as at the end. Clearly, in order to fully appreciate the gain of the OF system, we have to replay the experiment more than one day, an example is detailed in Table 5.

Table 5 shows the evolution, over three days, of the difference between the average turnovers of both OF system and O system. We observe that the two systems generate the same turnover at the end of the first day. At the end of the second day, the OF system is better than O system and at the end of the third day, the difference between the average turnovers of the two systems is even more important than the end of the second day. This observation confirms the results of the previous experiment, in the sense that, the difference between the average turnovers of the two systems evolves from "+20 euros" at the end of the second day to "+60 euros" at the end of the third day which confirms that the gain is more appreciated on a long period. In other word, the OF system has visibility on the Internet prospects' flows over the future periods that reports some sales in order to perform them more expensive later. In the previous experiments, we studied the ability of the exact forecasting system to optimize the results of the multi-period assignment model. In this experiment, we address the influence of the forecasting error on the assignment model. Figure 3 shows that the integration of Internet prospects forecasting of the orders' groups having a forecasting error until 10% improves the results in terms of turnover. However, orders' groups with forecasting error more than 10% deteriorate the results.

CONCLUSION

In this paper, we presented a forecasting system that allows improving the results of a multi-period optimization model. Thus, the challenge behind this work was to determine the right data to forecast (not only a number of internet prospects) in order to integrate them into the optimization model. In other word, the proposed forecasting system returns the number of Internet prospects arriving at each period as well as their characteristics in order to determine orders for which these Internet prospects can be assigned. In addition, our forecasting system ensures the fact that forecasting's time step and optimization model's period must be the same in order to articulate both forecasting model and optimization model. Finally, in order to assess the benefit of using such an approach, empirical investigation was conducted on real data from PdL data base. Numerical tests proved that forecasting model is effective for realistic data amount. Hence, our contribution can be directly integrated into PdL platform as well as other real-world applications which are similar, for example, call centers and service companies.



  Set of Internet prospects: P = {j: j ∈ {1, .., m t' }}.  Set of orders: O = {i: i ∈ {1, .., n}}.  Set of clients: C = {k: k ∈ {1, .., r}}. The model P represents a deterministic linear mathematical program which optimizes, at each time step t, Internet prospects' assignment to orders by considering several future periods, whose components are presented in what follows.

  For i ∈ O do { if x ijt = 1 then { l := l ∪ {i}}; } return l; End

Figure 3 :

 3 Figure 3: influence of the forecasting error rate on the average daily turnover

Time_Series (): Data:

 T = {t: t ∈ {1,..., d}} is set of historical periods.  P = {j: j ∈ {1, ..., m t }} is set of Internet prospects.  O = {i: i ∈ {1, .., n}} is set of orders.  δ = {l: l ∈ {1, .., e}} is set of orders' groups.  A = {y lt : l ∈ δ, t ∈ T} such that y lt is Internet prospects' number shared by orders' group l at period t.  X = {x ijt : i ∈ O, j ∈ P, t ∈ T} such that x ijt' = 1 if Internet prospect j corresponds to the characteristics of order i at period t, 0 otherwise. The function Time_Series() generates the set of the time series. For each orders' group l, we have the time series (y lt ) t∈T which is the Internet prospects' number of the orders' group l at each historical period t. The orders' group sharing the Internet prospect j at the period t is computed by the function Order_Group(j , t , O , X).

Begin

Order_Group (j , t , O , X):