title: A risk-aversion approach for the Multiobjective Stochastic Programming problem
authors: Javier León Justo Puerto Begoña Vitoriano
text: 
  
  
  
  
  
  
  1,*]Javier León
  2]Justo Puerto
  1]Begoña Vitoriano
  
  [1]HUM-LOG Research Group, Instituto de Matemática Interdisciplinar (IMI), Facultad de Ciencias Matemáticas, Universidad Complutense de Madrid, Plaza de las Ciencias 3, Madrid, 28040, Spain
  [2]Mathematical Research Institute (IMUS), University of Seville, Sevilla, Spain
  [ ]<javileon@ucm.es>, <puerto@us.es>, <bvitoriano@mat.ucm.es>
  [*]Corresponding author
  
  A risk-aversion approach for the Multiobjective Stochastic Programming problem
      [
      Received December 30, 2023; accepted December 30, 2023
  ==============================================================================
  
  
  	Multiobjective stochastic programming is a field well located to tackle problems arising in emergencies, given that uncertainty and multiple objectives are usually present in such problems. A new concept of solution is proposed in this work, especially designed for risk-aversion solutions. A linear programming model is presented to obtain such solution.
  
  
  Keywords:
  Multiobjective stochastic programming; Linear programming; Risk aversion
  
  
  
  
  
  § INTRODUCTION
  
  Vida real, múltiples criterios, aversión al riesgo
  
  
  Decision making is never easy, yet we often have to make decisions. Emergencies and disaster management are fields in which many difficulties often arise, such as high uncertainty and multiple conflicting objectives. To overcome such difficulties, risk-aversion decisions are usually sought. Risk-aversion is the attitude for which we prefer to lower uncertainty rather than gambling extreme outcomes (positive or negative).
  
  Risk-aversion, although typically studied in problems with uncertainty, can as well be considered when making decisions with multiple criteria. For instance, in the field of disaster management, solutions that are sufficiently good for all criteria are usually preferred than others that perform exceptionally good for some criteria but inadequately for the others.
  
  Multicriteria decision making (MCDM) is a field worth of consideration when studying real-world problems. Such is the case that MCDM techniques have been recently used for solving problems as varied as:
  disaster management <cit.>, engineering <cit.>, finance <cit.>, forest planning <cit.>, healthcare <cit.>, location of waste facilities <cit.>, police districting <cit.>, route planning <cit.>, train scheduling <cit.> or urban planning <cit.>.
  
  This situation, in which multiple conflicting objectives have to be optimized, has led to the definition of different solution concepts and methodologies. Depending on the problem and the type of solution considered, a specific methodology should be applied.
  
  The concept of efficiency reflects the intuition that for a solution to be acceptable, another cannot exist improving that one in every objective. Multiple notions of efficiency are available. The notation that this paper follows is the given in <cit.>.
   
  	Let f_1(x),…,f_K(x) be objective functions to be minimized, and let X be the feasible set. A feasible solution x̂∈ X is called:
  
  	
  		
    * Weakly efficient if there is no x ∈ X such that f(x) < f(x̂) i.e. f_k(x) < f_k(x̂) for all k=1,…,K.
  		
    * Efficient or Pareto optimal if there is no x ∈ X such that f_k(x) ≤ f_k(x̂) for all k=1,…,K and f_i(x) < f_i(x̂) for some i ∈{1,…,K}.
  		
    * Strictly efficient if there is no x ∈ X, x ≠x̂ such that f(x) ≤ f(x̂).	
  	
  
  
  Furthermore, the set of efficient solutions is called the efficient set, and the image under f of this set is the nondominated set. It is reasonable to assume that the solution given to any problem must lie in the efficient set. 
  
  
  
  
  
  
  
  
  
  Uncertainty is another feature present in the studied problems, in which risk-averse decisions will be preferred. The most common ways for dealing with the uncertainty are stochastic programming and robust optimization, in which fuzzy optimization is also included <cit.>.
  
  The different approaches for treating uncertainty do not respond to the desires of the modeller, instead, they reflect the nature of the uncertainty.
  If the uncertainty comes with an underlying known or estimated probability distribution, then stochastic programming is used.
  On the other hand, if uncertainty comes from a lack of precision or semantic uncertainty, then robust optimization is used.
  Robust optimization does not assume a known (or existing) distribution <cit.>.
  A recent review of robust optimization is written in <cit.>.
  For an introduction to stochastic programming the reader is referred to <cit.>.
  
  Stochastic programming is the widest used technique when there are historical data or information to infer a probability distribution. Moreover, usually discrete distributions are used, calling scenarios the different values. The concepts of value-at-risk (VaR) and conditional value-at-risk (CVaR) are widely used for quantifying risk (see for instance <cit.>).
  They are typically defined for losses distributions in finance, where the right tail of the distributions are of interest.
  
  
  
  
  
  
  
  	Given F_X(x) distribution function, and β∈ [0,1], the β-CVaR is the conditional expectancy over {x: F_X(x) ≥β}.
  
  
  
  
  
  Consider now the following problem, in which multiple objectives to be minimized and uncertainty are included simultaneously:
  
      min_x ∈ X( f_1(x,ω), …, f_K (x,ω))
  
  The above problem is typically called multiobjective stochastic programming problem (MSP), especially if ω, the uncertainty source, has a known probability distribution.
  
  
  In this paper we introduce a new solution concept in multiobjective stochastic programming based on risk-aversion preferences. Such concept is complemented with a mathematical programming model to efficiently compute it, and computational experiments are performed to assess its strengths.
  
  
  
    
  Structure of paper
  The remaining of this paper is organized as follows. Section <ref> includes the definition of a novel concept of solution for MSP problems and studies its properties. Section <ref> illustrates a basic example of how this solution can be found if the decision space is finite and small. 
  
  Section <ref> shows how to obtain such a solution with a linear programming model. An application to the multicriteria knapsack problem is developed in Section <ref>. 
  
  
  
  § MULTIOBJECTIVE STOCHASTIC PROGRAMMING
  
  
  
   §.§ Literature review
  
  <cit.> develops PROTRADE method, where utility functions are defined to aggregate objectives into a single objective stochastic problem. The resulting problem is solved with an interactive method, where the decision-maker defines an expected solution and a feasibility probability.
  <cit.> reduces the stochasticity by adding some good measures to the list of objectives, such as the mean, variance, or probability of being over/below a threshold.
  
  
  The resulting multiobjective deterministic problem is solved aggregating the objectives, but it could be solved via other techniques.
  
  <cit.> compare the stochastic approach with the multiobjective approach when using different techniques. The stochastic approach transforms the MSP on a single-objective stochastic problem, while the multiobjective approach first reduces the stochasticity transforming the MSP on a deterministic multiobjective problem. They highlight that “the multiobjective approach forgets the possible existence of stochastic dependencies between objectives.”
  <cit.> study stochastic goal programming, where the deviation of the objective functions to some goals set beforehand to stochastic values is minimized.
  
  
  
  In <cit.> a chance-constrained compromise approach is proposed, with an example presented in <cit.>.
  
  
  
  
  
  
  
  In <cit.> the INTEREST method is proposed. It is an interactive reference point method. The decision-maker gives reference levels u_i and probabilities β_i, hoping to achieve a solution x^* such that ℙ(f_i(x^*) ≤ u_i ) ≥β_i. If this is infeasible, the decision-maker should either increase the reference levels or decrease the probabilities of achievement.
  <cit.> reviews different solutions methods for the MSP problem, categorizing them as stochastic approach or multiobjective approach.
  
  Some fields where MSP models have been developed are: forest management <cit.>, multiple response optimization <cit.>, energy generation <cit.>, energy exchange <cit.>, capacity investment <cit.>, disaster management <cit.>, portfolio optimization <cit.> and cash management <cit.>, among others
  
  
  
   §.§ Definitions and dominance relationship
  
  The concept of CVaR allows to aggregate several scenarios by just looking at what happens in the worst ones. The ordered weighted averaging (OWA) operators are defined in <cit.>, and independently in the field of locational analysis <cit.> under the name of ordered median function. These concepts will allow us to aggregate different criteria by looking at the least desirable ones, as a risk-aversion measure.
  
  
  	Given a_1,…,a_n ∈ℝ, the ordered weighted averaging (OWA) operator with weights λ_1,…,λ_n is defined as:
  	
      OWA(a_1,…,a_n) = ∑_i λ_i a_(i)
  
  	where (a_(1),…,a_(n)) is the ordered vector from largest to smallest (a_1,…,a_n).
  
  
  
  	For certain weights, the OWA represents a known quantity:
  	
  		
    * If λ_i = 1/n, the resulting OWA is the average of a.
  		
    * If λ_1 = 1, and λ_j = 0 for j>1, the OWA is the maximum of a.
  		
    * If λ_n = 1, and λ_j = 0 for j<n, the OWA is the minimum of a.
  	
  
  
  <cit.> later study how to assign weights for an OWA when criteria have different importances.
  
  
  	Given a_1,…,a_n ∈ℝ with importances u_1,…,u_n such that ∑_i u_i = 1 the weights λ_j for the OWA can be calculated with f, the weight generating function in the following manner:
  	
  		
    * Sort vector a such that a_(1)≥ a_(2)≥…≥ a_(n).
  		
    * With (·) as the order induced by a, define T_j = ∑_k=1^j u_(k).
  		
    * Let f be a function such that f:[0,1]→[0,1] and f(0)=0, f(1)=1. This function is called weight generating function.
  		
    * Obtain the weights as λ_j = f(T_j) - f(T_j-1).
  	
  
  
  
  [of Definition <ref>] 
  	Consider the following weight generating function, for a given r ∈ (0,1]:
  	
      f(x) = x/r   if  x<r
      
      		1              if  x≥ r
  
  
  	Let (·) be the order such that a1≥…≥ an, uj the weight associated to a j, and also let T_j = ∑_k=1^j u_(k).
  	We shall see know how the weights are obtained from f. Let j^* be such that T_j^*-1 < r ≤ T_j^*.
  	
  		
    * λ_1 = f(T_1) = f(u1) = u1/r, assuming u1<r
  		
    * λ_2 = f(T_2) - f(T_1) = f(u1 + u2) - f(u1)= u1+u2/r - u1/r = u2/r, assuming u1 + u2<r
  		
    * …
    * λ_j^* = f(T_j^*) - f(T_j^*-1) = 1 - (u1+u2+…+uj^*-1/r), since T_j^*≥ r
  		
    * λ_j^*+1 = f(T_j^*+1) - f(T_j^*) = 1 - 1 = 0
  		
    * …
    * λ_n = f(T_n) - f(T_n-1) = 1 - 1 = 0
  	
  
  	Consequently the OWA of a_1,…,a_n with importances u_1,…,u_n is:
  	
      OWA =    u1/ra1 + u2/ra2 + … + [1 - (u1+u2+…+uj^*-1/r)] aj^*
      
      	      =    u1a1 + u2a2 + … +(r - u1 - u2 - …)aj^*/r
  
  	That is, the OWA is the average of the worst a_j, weighted by their importances, with total importance adding up to r
  
  
  
  The starting point of this paper is the recurrent idea of representing ordered weighted or ordered median operators by means of k-sums. k-sums (or k-centra in the location analysis literature) are sums of the k-largest terms of a vector <cit.>. One can trace back, at least to <cit.>, the use of k-sums to represent ordered median objectives. More recent references are for instance <cit.> and <cit.>. This last reference introduces a normalized version of k-centrum, named β-average that will be used in our paper.
  
  Through the remaining of the paper consider that f_k^j(x) are functions to be minimized within a feasible set X, with k = 1,…,K representing K different objectives with importances w_k and j = 1,…,J encoding J different scenarios with probabilities π_j.
  
   
  	
  	Given β∈ (0,1], for each criterion k it can be defined g^β_k(x) which measures the average of f on the worst scenarios (f_k^1(x),…,f_k^J(x)), with accumulated probability equal to β.
  
  
  
  [<cit.>]
  	Given a value β, if the sum of the probabilities of the worst scenarios is exactly β, then the β-average is exactly ( 1 - β)-CVaR.
  
  
  
  
  
  	Consider a point x, a fixed criterion k and 5 different scenarios with probabilities π_j and values of f_k^j given. Table <ref> shows the β-averages for different values of β, in which the scenarios have been ordered from largest value of f to smallest.
  	
  		
    * For β = 0.2, the scenario j=1 is the only one needed to obtain the worst scenario with probability 0.2, and hence g^β_k(x) = 0.2 · 10/0.2 = 0.2.
  		
    * When β equals 0.3 it is necessary to include scenario 2, obtaining a β-average of 0.2 · 10 + 0.1 · 7/0.3 = 9.
  		
    * Finally if β = 0.5 scenario 3 needs to be added as well, but only with the probability needed until reaching 0.5: g^β_k(x) = 0.2 · 10 + 0.1 · 7 + 0.2· 4/0.5 = 7.
  	
  
  	
  
  
  When using the β-average the functions f_k^j(x) were transformed into g_k^β(x), a collection of K functions not depending on the scenario. An OWA will be defined now, via its weight generating function, that will reduce the K β-averages into a scalar function.
   
  	Given x_i ∈ℝ with importance w_i (i=1,…,K, w_i ≥ 0, ∑_i w_i = 1) and r ∈ (0,1], the function O_r(x) is defined as the OWA with the following weight generating function:
  	
      f(x) = x/r   if  x<r
      
      		1              if  x≥ r
  
  
  We should be careful with the multiple meanings of w (importances for criteria and weights for OWA's). Also with f_k^j(x) as the original functions and the weight generating function
  
  	The definition of O_r(x) is made on a similar manner that the one given of the β-average (Definition <ref>), but it is done on a context with importances rather than probabilities. Example <ref> shows the similarities between both approaches.
  
  
  
  	Consider a point x and let g_k(x) be the evaluation of x under 5 different criteria with importances w_j. Table <ref> shows the r-OWAs for different values of r, in which the criteria have been ordered from largest values of g_k(x) to smallest. Consider the case r=0.5:
  	
  		
    * As g_k(x) are already ordered for largest to smallest, the values of T_k are:
  		
      T_1 = 0.2,T_2 = 0.2+0.1 = 0.3, T_3 = 0.6, T_4 = 0.85, T_5 = 1
  
  		
    * The values of T_k under f:
  		
      f(T_1) = 0.2/0.5, f(T_2) = 0.3/0.5, f(T_3) = f(T_4) = f(T_5) = 1
  
  		
    * The weights of the OWA:
  		
      λ_1 = 0.2/0.5, λ_2 = 0.3-0.2/0.5 = 0.1/0.5, λ_3 = 1-0.3/0.5 = 0.2/0.5, λ_4 = λ_5 = 0
  
  		
    * Consequently the r-OWA is:
  		
      r-OWA = 0.2x_(1)+0.1x_(2)+0.2x_(3)/0.5 = 0.2· 10+0.1· 7+0.2· 4/0.5 = 7
  
  	
  
  	
  
  
  
   
  	Given x_1,…,x_K and its associated importances w_1,…,w_K, then the λ_k of the r-OWA are determined in such a way that:
  
  	
      O_r(x) = max{λ̃_1 x_1 + … +λ̃_K x_K/r|λ̃_k ≤ w_k, ∑λ̃_k = r }  with λ_k = λ̃_k/r
  	
  
  Maybe this needs to be explicitly proven
  
  
  Given r, β∈ (0,1] and x ∈ X, let us introduce the function h^β_r(x) as the r-OWA of the β-averages. That is:
  
      h^β_r(x) = O_r(g_1^β(x),…,g_K^β(x))
  
  
  
  I think β and r could be given some names. For instance, uncertainty and multicriteria risk-seeking parameters. When β is small very few scenarios are considered, so we are very protected against uncertainty risk. Same goes with r.
  
  
  
  
  
  
  	If the importance of all criteria is the same (w_k = 1/K for all k) and r = n/K with n ∈{1,…,K}, then the h^β_r(x) is the average of the n worst β-averages. Recall that this is called n-centra <cit.>.
  
  
  
   
  	Let x and y feasible solutions (x,y∈ X) and r,β∈ (0,1]. Then x dominates y (x ≿ y) if h^β_r(x) ≤ h^β_r(y), where h^β_r(x) is the r-OWA of the β-averages.
  
  
  Definition <ref> induces a domination relationship with the following properties:
  
  
  	
    Reflexivity Given x, h^β_r(x) ≥ h^β_r(x), and then x ≿ x, so ≿ is reflexive.
  	
    Transitiveness Given x ≿ y, y ≿ z, we have h^β_r(x)≥ h^β_r(y) y h^β_r(y)≥ h^β_r(z), and then h^β_r(x) ≥ h^β_r(z), which leads to x ≿ z, and we conclude that ≿ is transitive.
  	
    Antisymmetry Given x ≿ y, y ≿ x, we have h^β_r(x)≥ h^β_r(y) and h^β_r(y)≥ h^β_r(x), but from h^β_r(x) = h^β_r(y) it cannot be guaranteed that x=y, and hence ≿ is not antisymmetric.
  
  
  
  
   §.§ Idea of solution and dominance properties
  
  Consider the multiobjective stochastic programming problem:
  
      min_x ∈ X( f_1(x,ω), …, f_K (x,ω))
  
  
  The previously defined concepts of β-average and r-OWA transform the MSP problem into a deterministic multiple objective problem, and then into a deterministic single objective problem.
  
  
  
  0.45
  	
      MSP → MOP → LP (MIP)
  
  	
      f_k^j(x)  g_k^β(x)  h^β_r(x)
  
  
  0.50
  	
  		
    * For every x ∈ X there is a function f_k^j to be minimized which depends on the scenario j and the criterion k.
  		
    * The problem is transformed into a deterministic one with multiple objectives (MOP) using the β-average concept.
  		
    * Computing the r-OWA, each x ∈ X is assigned a scalar. The problem consists of finding the x which minimizes this h^β_r(x).
  	
  
  
  
  
  The solution procedure lies into what is usually called a scalarization approach. When obtaining a minimizer of h^β_r(x) it is also desired that the optimal solution is efficient for the associated MOP problem:
  
  
      min_x ∈ X(g^β_1(x),…,g^β_K(x)) MOP
  
  
  
  
  	Given x^* minimum of h^β_r(x) the following statements hold:
  	
  		
    * x^* is not necessarily efficient of the MOP problem.
  		
    * x^* is weakly efficient of the MOP problem.
  		
    * If x^* is the only minimum of h^β_r(x), then x^* is efficient.
  		
    * Given x^* not efficient, an alternative y^* can be found on a second phase such that y^* is efficient and h^β_r(x^*) = h^β_r(y^*).
  	
  
  
  These properties are known when using scalarization techniques <cit.>. Hence only an example of the first statement will be shown.
  
  [x^* is not necessarily efficient]
  	Consider the example displayed on Table <ref>, in which there are only two feasible solutions, two equiprobable scenarios (π_1 = π_2 = 12), three equally important criteria (w_1 = w_2 = w_3 = 13), and consider the values of β = 12 and r = 23 are taken.
  
  	
  
  	The β-averages are (0.8,0.4,0.65) for the first alternative and (0.8,0.45,0.65) for the second alternative. When computing the function h^β_r, both alternatives have an objective value of 0.725. Consequently, even though the second alternative is an optimal solution of h^β_r, it is not an efficient solution of the MOP problem as its β-averages are dominated by those of the first alternative.
  
  
  
  
  
  
  
   §.§ An illustrative example
   
  The solution concept proposed will be now applied, first with a discrete (and small) case. When the solution space is discrete, and all feasible solutions can be explicitly enumerated, the steps are as follows:
  
  
  	
    Step 0 Normalize all objective functions f^j_k(x).
  	
    Step 1 Set values for β,r ∈ (0,1].
  	
    Step 2 For every x∈ X and every criterion define g_k^β(x) as:
  	
  	
      g_k^β(x) = 8cmaverage of worst scenarios for criterion k
       with probabilities adding up to β
  
  	
    Step 3 Define h_r^β(x) as:
  	
  	
      h_r^β(x) = 8cmaverage of worst g_k^β(x) values
      with importances adding up to r
  
  	
    Step 4 Search for x ∈ X minimizing h_r^β(x)
  
  
  Assume a decision space with only four alternatives, evaluated under five different scenarios with six criteria. For each of those alternatives it can be computed the value of the functions f_k^j(x) to be minimized. Table <ref> shows the values of f, evaluated on feasible point x_1, by each of the scenarios and criteria considered.
  
  
  
  The first step consists on calculating the β-averages. Let assume a value of β = 0.3:
  
  	
    * For the first criterion the worst scenario is j_5, which has probability 0.1. The second worst is j_4, with a probability of 0.25. As the sum of those probabilities exceeds the β fixed, for computing the β-average just a probability of 0.2 is considered:
  	
      g_1^β(x_1) = 0.1 · 0.86 + 0.2 · 0.76/0.3 = 0.793
  
  	
    * g_2^β(x_1) = (0.2 · 0.65 + 0.1 · 0.44  )/0.3 = 0.580
  	
    * g_3^β(x_1) = (0.3 · 0.90 )/0.3 = 0.900
  	
    * g_4^β(x_1) = 0.833, g_5^β(x_1) = 0.930, g_6^β(x_1) = 0.728
  
  
  The last step is calculating the function h_r^β(x), that is, the r-OWA of the β-averages. Table <ref> calculates the r-OWA, and shows as well the information of the previously calculated β-averages, when the value of r = 0.17 is taken. 
  
  
  
  
  
    
  Results The values of the functions for the other alternatives, as well as its β-averages and r-OWAs are shown in Tables <ref>, <ref> and <ref>, starting on Page tab:inst1alt2. A summary of the results can be seen in Table <ref>, where all the β-averages and r-OWAs are shown, determining that the optimal alternative for the values of β and r given is Alternative 1.
  
  
  
  Variations on β and r yield very different results. Figure <ref> shows which of the four alternatives has the lowest h value, depending on the values of β and r. 
  
  Figure <ref> shows the optimal objective value when varying the parameters β and r. It can be appreciated how h decreases when β and r increase. This is due to the fact that the original f_k^j functions are to be minimized, and the larger the parameters β and r are, more favourable scenarios/criteria will take part on the computation of h_r^β(x), hence decreasing its optimal value.
  
  
  
  
  
  § COMPUTING THE MINIMUM: CONTINUOUS CASE
  
  A concept of solution was proposed with Definition <ref>. When the functions f_k^j(x) to be minimized are given, a new function h_r^β(x) to be minimized is defined, with parameters β and r such that h_r^β(x) is the r-OWA of the β-averages. If the decision space is sufficiently small, the procedure shown in the above example obtains such a solution.
  
  In this section, a mathematical programming model will be developed to obtain the minimum of h_r^β(x) which allows one to obtain the proposed solution for bigger decision spaces, including continuous ones.
  
  
  
  
  
   §.§ Mathematical programming model
  
  Given k and x ∈ X we have the vector (f_k^1(x),…,f_k^J(x)). Let (f_k^(1)(x),…,f_k^(J)(x)) be the ordered vector such that f_k^(j_1)(x) ≥ f_k^(j_2)(x) when j_1 ≤ j_2.
  
  Given β∈ (0,1], let  be the ordered scenario such that:
  
      ∑_j = 1^π_(j)≥β,     ∑_j = 1^ -1π_(j) < β
  
  
  Alternatively:
  
  
      f_k^(1)(x) ≥ f_k^(2)(x) ≥…≥ f_k^()(x) ≥ f_k^(+1)(x) ≥…≥ f_k^(J)(x)
  
  
      1 = π_(1) + π_(2) + … + π_(-1)_< β + π_()_≥β + … + π_(J)
  
  
  Also let:
  
      π̂_j = π_j                                          j ∈{(1),…,(-1) }
      β - ∑_j=(1)^j=( -1)π_j     j = 
      
      	0                                                  otherwise
  
  The definition of π̂_ is made in such a way that ∑_j π̂_j = β. In this way, the average of the β worst values can be computed as 1/β∑_j=1^J π̂_j f_k^(j) (x), which coincides with the definition of β-average (Definition <ref>). This computation can be written as the following optimization problem:
  
      max_ũ_j     1/β∑_j = 1^J ũ_j · f_k^j (x) 
      s.t.     ∑_j=1^J ũ_j = β
          0 ≤ũ_j ≤π_j      j = 1,…,J
  
  
  A more natural approach would be to consider u_j = ũ_j/β. These u_j represent the proportion in which scenario j plays a part on the aggregated β-average. Introducing that change, the model is:
  
  
      max_u_j     ∑_j = 1^J u_j · f_k^j (x) 
      s.t.     ∑_j=1^J u_j = 1           
          0 ≤ u_j ≤π_j/β     j = 1,…,J
  
  
  The dual formulation is:
  
  
      min_z,y_j      z + ∑_j=1^J π_j/β y_j 
      s.t.      z + y_j ≥ f_k^j(x)      j = 1,…,J           
          z  free, y_j ≥ 0
  
  
  And hence finding the x ∈ X which minimizes the average of the worst β scenarios for a given k is:
  
  
      [ min_x ∈ X ][              rl)
      		max_ũ_j 1/β∑_j = 1^J ũ_j f_k^j (x);                       s.t.            ∑_j=1^J ũ_j = β;                              0 ≤ũ_j ≤π_j      j=1,…,J;                            ]
  
  
  Or alternatively:
  
      [ min_x ∈ X ][                   rl)
      		min_z,y_j             z + ∑_j=1^J π_j/β y_j;                              s.t. z + y_j ≥ f_k^j(x)      j = 1,…,J;                                     z  free, y_j ≥ 0      j = 1,…,J;                                   ]
  
  
  Which is equivalent to:
  
  
  	min_z,y_j,x      z + ∑_j=1^J π_j/β y_j     
  
  	s.t.       z + y_j ≥f_k^j(x)       j = 1,…,J 
  
                             z  free, y_j ≥0       j = 1,…,J 
  
  	                        x ∈X    
  
  
  
  	Models (<ref>) and (<ref>) are equivalent, as for any x ∈ X chosen in (<ref>) the values z and y_j will get as small as allowed by constraint (<ref>), as this improves the objective function (<ref>). Consequently for every x, its β-average will be computed appropriately, and thus (<ref>) obtains the x∈ X with smallest β-average, as desired on (<ref>). 
  
  
  For every k ∈{1,…,K} thanks to the problem (<ref>) the function g^β_k(x) can be defined, which measures for each x∈ X the β-average for that criterion, being:
  
  
      2
      	g^β_k(x) ≡min_z_k,y_kj    z_k + ∑_j=1^J π_j/β y_kj      
      s.t.     z_k + y_kj≥ f_k^j(x)          j = 1,…,J     
          z_k  free, y_kj≥ 0          j = 1,…,J
  
  
  The already known approach for single criterion problems ends here. Given that, the next step is finding a “good” solution for all k. That is:
  
      min_x ∈ X(g^β_1(x),…,g^β_K(x))
  
  
  
  Given r ∈ (0,1] the r-OWA of the β-averages will be now computed (in accordance with the definition given in Section <ref>). That is, the solution of the following problem is sought:
  
      max_t̃_k     1/r∑_kt̃_k · g^β_k(x) 
         ∑_k t̃_k = r         
          0 ≤t̃_k ≤ w_k      k = 1,…,K
  
  
  Or equivalently: 
  
      max_t_k     ∑_kt_k · g^β_k(x) 
         ∑_k t_k = 1         
          0 ≤ t_k ≤w_k/r     k = 1,…,K
  
  
  Its dual formulation is:
  
  	min_z,v_k       z + ∑_k w_k/rv_k     
  
  	s.t.        z + v_k ≥g^β_k(x)       k = 1,…,K
  
  	                       z  free, v_k ≥0       k = 1,…,K
  
  
  Replacing the value of g_k^β(x) given in (<ref>) the next model is obtained:
  
  
  
      min_z,v_k      z + ∑_k w_k/rv_k 
      s.t.      z + v_k ≥[l]
      						       min_z_k,y_kj    z_k + ∑_j=1^J π_j/β y_kj
        s.t.     z_k + y_kj≥ f_k^j(x)   ∀ j 
          z_k  free, y_kj≥ 0 
      	                    ∀ k 
          z  free, v_k ≥ 0     ∀ k
  
  
  
  Model (<ref>) calculates for a given x∈ X the r-OWA of its β-averages, which coincides with the notion of the function h(x) given in Section <ref>. This problem is not explicit in that it contains nested optimization problems in the constraints. For that reason, we propose a single level alternative for x ∈ X fixed.
  
  Consider the following linear programming model:
  
  
  	min_z,v_k,z_k,y_kj       z + ∑_k w_k/rv_k      
  
  	s.t.       z + v_k ≥z_k + ∑_j=1^J π_j/β y_kj       ∀k  
  
  	   	                 z_k + y_kj ≥f_k^j(x)       ∀k,j    
  
  	                       y_kj ≥0       ∀k,j    
  
  	                       z_k  free, v_k ≥0       ∀k   
  
  	                       z  free	     
  
  
  
  
  
  
  
  
  	Transformation from model (<ref>) to model (<ref>) is valid, in that their optimal solution and objective values coincide. 
  
  
  
  	Let (z^*,v_k^*,z_k^*,y_kj^*) be the optimal solution of model (<ref>). (z^*,v_k^*) is feasible of model (<ref>), and it will be shown that it is also optimal for such model.
  	Assume it exists (z',v_k') feasible of model (<ref>) with:
  	
      z' + ∑_k w_k/rv_k' < z^* + ∑_k w_k/rv_k^*
  
  	This and constraint (<ref>) implies there exists k_0 such that:
  	
      z' + v_k_0' < z_k_0^* + ∑_j=1^J π_j/β y_k_0j^*
  
  	otherwise (z',v_k',z_k^*,y_kj^*) would be optimal of model (<ref>).
  	Since z_k_0^* and y_k_0j^* are feasible of model (<ref>) they are also feasible of the model on the RHS of constraint (<ref>), and thus z' and v_k_0' violate constraint (<ref>).
  
  
  Proposition <ref> showed that the optimal solutions of models (<ref>) and (<ref>) coincide. Proposition <ref> goes further showing the connection between their feasible sets.
  
  
  	The feasible set of model (<ref>) is a projection of the feasible set of model (<ref>).	
  
  
  
  
  	
  	
  		
    * For each feasible solution (z,v_k) of model (<ref>) there is at least one feasible solution of model (<ref>) with same values (z,v_k), being so the same objective function.
  
  		Let (z^1,v^1_k) a feasible solution of model (<ref>), and (z^*_k,y^*_kj) the optimal solution where the minimum of the right-hand side of equation (<ref>) is reached for each k. Since constraints (<ref>), (<ref>), (<ref>) and (<ref>) are satisfied in model (<ref>), (z^1,v^1_k,z^*_k,y^*_kj) is a feasible solution or model (<ref>).
  
  		
    * For each feasible solution (z,v_k,z_k,y_kj) of model (<ref>), (z,v_k) is a feasible solution of model (<ref>), being so the same objective function.
  		Let (z^2,v^2_k,z^2_k,y^2_kj) a feasible solution of model (<ref>). Since constraints (<ref>), (<ref>) and (<ref>) are included in model (<ref>), (z^2_k,y^2_kj) is feasible for the model included in the RHS of constraint (<ref>) and therefore greater than or equal to the minimum of that model, verifying:
  		
      z^2 + v^2_k ≥ z^2_k + ∑_j=1^J π_j/βy^2_kj≥min{ z_k + ∑_j=1^J π_j/βy_kj}
  
  		and so, feasible for model (<ref>).
  	
  
  
  Finally after proving the validity of model (<ref>) it is possible to let x ∈ X free, with the purpose of finding the one minimizing the function h_r^β(x):
  
  
  	min_z,v_k,z_k,y_kj,x       z + ∑_k w_k/rv_k      
  
  	s.t.       z + v_k ≥z_k + ∑_j=1^J π_j/β y_kj       ∀k  
  
  	   	                 z_k + y_kj ≥f_k^j(x)       ∀k,j    
  
  	                       y_kj ≥0       ∀k,j    
  
  	                       z_k  free, v_k ≥0       ∀k   
  
  	                       z  free	      
  
  	                       x ∈X    
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § KNAPSACK PROBLEM
  
  The multiobjective stochastic knapsack problem is used to illustrate the usefulness of the previously defined concept.
  
  	Let I be a collection of objects with weights v_i, which can be selected as members of a knapsack with maximum weight V. There is a set of scenarios J, each of them with probability π_j, and a set of criteria K, with importances w_k. For every pair of scenario-criterion, there is a benefit associated with selecting object i, denoted by b^i_jk. Which objects should be taken in order to maximize benefit?
  
  
  The above problem differs with the well-known knapsack problem in that there is stochasticity and multiple objectives to be maximized.
  
  The following MSP model can be adapted to analyze the problem. Note that to preserve the sense of the optimization, rather than to maximize the benefits of the carried objects, it will be minimized the value of the objects not chosen.
  
      min_x_i     { f_k^j(𝐱):= ∑_i (1-x_i)b^i_kj}
      s.t.     ∑_i v_i x_i ≤ V      ∀ i    
          x_i ∈{0,1}     ∀ i
  
  
  When using the methodology developed in the previous sections, problem (<ref>) is transformed into the following mixed-integer linear programming model:
  
  
      MSPmin_z,v_k,z_k,y_kj,x_i      z + ∑_k w_k/rv_k     
      s.t.      z + v_k ≥ z_k + ∑_j=1^J π_j/β y_kj     ∀ k  
          z_k + y_kj≥∑_i (1-x_i)b^i_kj     ∀ k,j    
         ∑_i v_i x_i ≤ V      ∀ i    
          y_kj≥ 0      ∀ k,j    
          x_i ∈{0,1}     ∀ i    
          z_k  free, v_k ≥ 0      ∀ k   
          z  free
  
  
  For given r,β∈(0,1], model (<ref>) obtains the 𝐱^* minimizing the r-OWA of the β-averages. In order to illustrate the benefits of using model (<ref>), a naive way of solving problem (<ref>) is considered:
  
  
      MIPmin_x_i     ∑_k,j w_kπ_j ∑_i (1-x_i)b^i_kj
      s.t.     ∑_i v_i x_i ≤ V      ∀ i    
          x_i ∈{0,1}     ∀ i
  
  
  Hence model (<ref>) computes the average of the f_k^j, using the importances of the criteria and the probability of the scenarios. It is clear that for “average” criteria-scenarios , the optimal solution of model (<ref>), outperforms , the optimal solution of model (<ref>). Conversely  will improve  in unfavourable criteria-scenarios, as expected of a risk-averse solution.
  
  
  
   §.§ Computational experiments
  
  The following sections will show computational experiments, for different values of r and β and different number of objects, scenarios and criteria. Algorithm <ref> shows how the random instances are created, given a number of objects, scenarios and criteria.
  
  
  
  For each of the solved instances it will be recorded:
  
  	
    * t_MSP,t_MIP: Solution time in seconds of models (<ref>) and (<ref>). With them the following value is calculated:
  	
      Δ_time := t_MSP/t_MIP    (time penalty factor)
  
  	Δ_time, the time penalty factor, indicates the increase of computing time when solving model (<ref>) rather than model (<ref>).
  	
    * ,: Optimal values of the models. 
  	
    * f_MSP(), f_MIP(): Objective value of  in model (<ref>) and vice versa.
  	
    * To grasp the difference between the MSP and the naive approach, the following will be calculated:
  	
      Δ_avg    := 100f_MIP() -/    (deteriorating rate)
      Δ_tail    := 100f_MSP() -/f_MSP()    (improvement rate)
  
  	These quantities reflect what is the effect of making decision  instead of .
  	Large values of Δ_avg indicate high penalties for making decision  instead of  in average scenarios-criteria. Similarly, the larger Δ_tail is, the higher benefit is obtained from making decision  in tail events.
  	They will be called deteriorating rate and improvement rate
  
  
  Models are solved in GAMS 26.1.0 with solver IBM ILOG CPLEX Cplex 12.8.0.0, using a personal computer with an Intel Core i7 processor and 16Gb RAM.
  
  
  
    
  Experiment 1set5
  First experiment will consist on a full factorial design, in which the values of |I|,|J|,|K|,r,β fall in these sets:
  
  	
    * |I| ∈{50,100,200}
  	
    * |J| ∈{5,25,100}
  	
    * |K| ∈{3,6,9}
  	
    * r ∈{0.33,0.5,0.67}
  	
    * β∈{0.05,0.1,0.5}
  
  For each tuple (I,J,K) random data will be generated, using algorithm <ref>, which will then be solved for every pair (r,β). All criteria and scenarios are given same importance and probabilities. That is, w_k = 1/|K|, π_j = 1/|J|. Time limit was set in two hours by instance, in which all but three of the 3^5 = 243 configurations were solved to optimality.
  
  
  
    
  Experiment 2set6
  For the next experiment 100 random instances will be created, keeping the values of |I|,|J|,|K|,r,β constant and equal to the median value of the previous experiment. That is, |I|=100,|J|=25,|K|=6,r=0.5,β=0.1. All criteria and scenarios are given same importance and probabilities. All 100 instances were solved to optimality.
  
  
  
  
   §.§ Results
  
  
  
    
  Experiment 1
  Table <ref> shows for each of the 243 instances the solution times of the MSP and the MIP models, and the deteriorating and improvement rates of using the MSP solution instead of the MIP solutions (measured in deviation to MIP solution).
  
  Table <ref> shows the correlations between times and rates with the parameters of the instance. It can be seen how the MSP solution has a higher impact when fewer scenarios are considered. In addition to that, it can be appreciated that the MSP solution times decrease when β increase, that is, when more scenarios are included in the β-average computation.
  
  
  
  
  
  
  This appreciation is confirmed by Table <ref>, in which it can be seen that the median time penalty factor (how many more times does it take to solve the MSP model than the MIP one) is much smaller when β=0.5 than when β=0.05.
  
  
  
  
  Solution times of the MSP model are alarmingly high for some instances, due to the fact that the admissible integrality gap has been set to zero. If that is relaxed, it can be seen that all of the 243 instances reach an integrality gap smaller than 5% in under 3 seconds, 2% in under 5 seconds and 1% in under 88 seconds.
  
  Table <ref> groups instances by r and β, and shows the deteriorating and improvement rates. It can be seen that the improvement rate (in the tail) is generally higher than the deteriorating rate (in the average), especially in cases with small r and β. 
  
  
  
  
  
  This claim is also supported with Figure <ref>, where each of the 243 instances is shown according to the values of Δ_avg and Δ_tail, and grouped by the values of (r,β). Almost all of the instances ara above the imaginary line Δ_avg = Δ_tail, which shows that considering the MSP solution improves in the tail more than it loses in the average situations. In addition to that, it can be seen that the largest improvements in the tail are on instances with β = 0.05 (one of the usual values taken for CVaR), and especially with the smallest values of r. When r and β grow the differences between the MIP and MSP solutions are reduced.
  
  
  
  
  
    
  Experiment 2
  Table <ref> contains the results for each of the 100 instances, all of them with constant parameters |I|=100,|J|=25,|K|=6,r=0.5,β=0.1.
  
  Table <ref> contains a summary of the results, where it is again seen that the improvements in the tail are better than the loses in the average situations. Although single instances might take a long computing time, the median MSP solution time (3.74s) is definitely satisfactory. It is worth mentioning that the models were implemented without providing any extra bounds or known cuts that could reduce solution times.
  
  
  
  Finally, figure <ref> shows the values of f_k^j(x), where x= in blue squares and x= in orange circles, for just one of the created instances. It can be appreciated how  performs better than  in average criteria-scenarios, but  is better with unfavourable situations. 
  
  
  
  
  
  
  § CONCLUSIONS
  
  In this paper a new concept of solution has been proposed for Multiobjective Stochastic Programming problems, focused on risk-aversion. As such, this concept can be particularly useful in real-life situations where there exists a great concern with respect to unfavourable situations, such as emergency management.
  
  The solution concept is supported by an efficient way to compute it by a Mathematical Programming problem. This model is linear provided that the underlying problem can be linearly representable. Numerical experiments have been conducted for validating this approach, solving a multiobjective stochastic knapsack problem.
  
  The research has also shown that the improvements in the tail (unfavourable situations) are consistently higher than loses on average situations, especially when small values of the parameters β and r are chosen. These differences, although clearly noticeable, are not as high as one could expect. This is possibly due to the randomness of the data. It is reasonable to assume that in actual real-life problems there are choices that are more conservative for every scenario and criterion, and thus being preferable for risk-aversion attitudes.
  
  Results showed that there is a clear increase in computational time; however this is arguably acceptable as a price to pay for being risk-averse. Furthermore, this could also be due to the random nature of the data. Nevertheless, it was also shown that allowing for even rather small integrality gaps (1%) leads to drastic improvement on the computing times.
  
  
  
  
  
  § ACKNOWLEDGEMENTS
  
  This work has been supported by the UCM-Santander grant CT27/16-CT28/16, the Government of Spain grants MTM2016-74983-C02-01 and MTM2015-65803-R, and grant H2020 MSCA-RISE 691161 (GEO-SAFE).
  
  
  
  
  
  
  
  § EXTRA FIGURES AND TABLES
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  bottom=1cm,top=1cm
  empty
  
  
  
  
  
  bottom=3cm,top=3cm
  
  
  