Modeling and Optimizing the Multi-Objective Portfolio Optimization Problem with Trapezoidal Fuzzy Parameters

: A common issue in the Multi-Objective Portfolio Optimization Problem (MOPOP) is the presence of uncertainty that affects individual decisions, e.g., variations on resources or beneﬁts of projects. Fuzzy numbers are successful in dealing with imprecise numerical quantities, and they found numerous applications in optimization. However, so far, they have not been used to tackle uncertainty in MOPOP. Hence, this work proposes to tackle MOPOP’s uncertainty with a new optimization model based on fuzzy trapezoidal parameters. Additionally, it proposes three novel steady-state algorithms as the model’s solution process. One approach integrates the Fuzzy Adaptive Multi-objective Evolutionary (FAME) methodology; the other two apply the Non-Dominated Genetic Algorithm (NSGA-II) methodology. One steady-state algorithm uses the Spatial Spread Deviation as a density estimator to improve the Pareto fronts’ distribution. This research work’s ﬁnal contribution is developing a new defuzziﬁcation mapping that allows measuring algorithms’ performance using widely known metrics. The results show a signiﬁcant difference in performance favoring the proposed steady-state algorithm based on the FAME methodology.


Introduction
The Portfolio Optimization Problem (POP) is always present in organizations. One key issue in POP's decision process is the uncertainty caused by the variability in the project benefits and resources. The latter situation arises the necessity of a tool for describing and representing uncertainty associated with real-life decision-making situations. The POP searches a subset of projects under a predefined set of resources that maximizes the produced benefits; its formal definition is as follows.
Let A be a finite set of N projects, each characterized by estimates of its impacts and resource consumption. A portfolio is a subset of A that can be represented by a binary vector x = x 1 , x 2 , . . . , x n that assigns x i = 1 for every financed project i, and x i = 0 otherwise. Let → z (x) = z 1 (x), z 2 (x), . . . , z p (x) be the vector of impacts resulting from the linear sum of the attribute values of each financed project in x, i.e., the vector of size p representing multiple attributes related to organizational goals that describe the consequences of a portfolio x. Assume w.l.o.g. that the higher an attribute's value is, the better. Then, Problem (1) formally defines POP.
Maximize z 1 (x), z 2 (x), . . . , z p (x) , x ∈ R F (1) Table 1 summarizes the main features of the previously described works. Column 1 cites the research work and the studied POP variant. Columns 2 to 7 show the considered features in the research works: the solution algorithms it proposed, the type of instances it solved, the performance metrics it used, if it integrated preferences in the search process, if it considered a static or dynamic POP's version, the type of parameters it used, and if it used a steady-state selection scheme or not. It is worth nothing that, from the information in Table 1, only approaches based on intervals address POP's variant with uncertainty, and none of them utilized a steady-state selection scheme. The Fuzzy Adaptive Multi-objective Evolutionary solution methodology (FAME) has had great success in many optimization problems; however, there is a lack of studies about its performance on the POP. The previous situations open an area of opportunity, addressed in this work, consisting of studying optimization approaches' performance derived from fuzzy numbers and steady-state selection schemes on their search process to solve the Multi-objective POP with uncertainty (MOPOP).
Evolutionary algorithms commonly use a generational selection scheme to update each generation's population; the process creates several offspring through genetic operators and combines them with the parents to form the next generation of individuals [10,14,15]. On the other hand, an algorithm using a steaty-state selection scheme produces a single offspring during the reproduction process to combine with the parents. The efficiency of the population's update process achieved by the latter method is advantageous for any research [16]. Hence, this work proposes a new method based on FAME and fuzzy numbers to handling uncertainty and obtaining more robust solutions in MOPOP; the approach mainly uses fuzzy trapezoidal sets to reflect a magnitude's imprecision. This work's main contributions are: (1) a new mathematical model for MOPOP that considers fuzzy trapezoidal parameters; (2) a new algorithm based on FAME to solve the proposed model; (3) two novel steady-state NSGA-II to solve this MOPOP's variant; and (4) a novel strategy to measure the performance of the fuzzy multi-objective algorithms with the commonly used real metrics.
The remaining structure of this paper is as follows. Section 2 includes some elements of the fuzzy theory used in this work. Section 3 describes a new mathematical model of the Portfolio Optimization Problem with Trapezoidal Fuzzy Parameters. Sections 4 and 5 contain the proposed steady-state algorithms: T-NSGA-II and T-FAME, respectively. Section 6 describes the computational experiments done to assess the performance of the algorithms. Finally, Section 7 presents the conclusions.

Elements of Fuzzy Theory
This section contains the main concepts of fuzzy theory used in this work.

Fuzzy Sets
Let X be a collection of objects x, then a fuzzy set A defined over X is a set of ordered pairs A = {(x, µ A (x))/x єX} where µ A (x) is called the membership function or grade of membership of x in A which maps X to the real membership subspace M [17]. The range of the membership function is a subset of the nonnegative real numbers whose supremum is finite. Elements with a zero degree of membership usually are not listed.

Generalized Fuzzy Numbers
A generalized fuzzy number A is any fuzzy subset of the real line R, whose membership function µ A (x) satisfies the following conditions [18]: 1.
µ A (x) is a continuous function from R to the closed interval [0, 1] 2.
µ A (x) = 0, for β < x < ∞ where 0 < w < 1, a, b, α, β are real numbers. We denote this type of generalized fuzzy number as A = (a, b, α, β, w) LR . When w = 1, the generalized fuzzy number is denoted as A = (a, b, α, β) LR . When L(x) and R(x) are straight lines, then A is a trapezoidal fuzzy number, and denoted as A = (a, b, α, β). When b = α, then A is a triangular fuzzy number, and denoted as A = (a, b, β).
A triangular membership function definition is as: A trapezoidal membership function definition is as:

Graded Mean Integration (GMI)
Graded mean integration [19] is a defuzzification method to compare two generalized fuzzy numbers. We compare the numbers based on their defuzzified values. The number with a higher defuzzified value is larger. The formula to calculate the graded mean integration of a trapezoidal number A is given by: For a trapezoidal fuzzy number A = (a, b, α, β), there is a more straightforward expression which is P(A) = (3a + 3b + β − α)/6.

Order Relation in the Set of the Trapezoidal Fuzzy Numbers
Given the trapezoidal fuzzy numbers A 1 and A 2 , then:

Multi-Objective Portfolio Optimization Problem with Trapezoidal Fuzzy Parameters
This section presents the proposed mathematical model for MOPOP with Fuzzy Trapezoidal Parameters. It offers a detailed description of the construction of the fuzzy trapezoidal instances used in this work to assess the proposed solution algorithms' performance. It also includes a description of how the fuzzy trapezoidal parameter' values participate in evaluating objective functions and the candidate solutions' feasibility when the solution algorithms search across the solution space.

Mathematical Model
Let n be the number of projects to consider, C the total available budget, O the number of objectives, c i the cost of the project i, b ij the produced benefit with the execution of the project i in objective j, K the number of areas to consider, M the number of regions, A min k and A max k the lower and upper limits in the available budget for the area k, and R min m and R max m the lower and upper limits in the available budget for the region m. The arrays a i and b i contain the area and region assigned to the project i.x = (x 1 , x 2 , . . . . . . ., x n ) is a binary vector that specifies the selected projects included in the portfolio. If x i = 1 then the project i is selected, otherwise it is not. Now we define the MOPOP with Fuzzy Trapezoidal parameters as follows: where Subject to the following constraints: In this model, all the parameters and variables in bold and italic are trapezoidal fuzzy numbers.
The objective function tries to maximize the contributions of each objective (6). We calculate each objective by adding all the selected projects' contributions in the binary vector (7). The constraint (8) makes sure that the sum of the costs required for all the selected projects does not exceed the available budget. The set of constraints (9) makes sure that the sum of the projects' costs is in the range of the involved areas' available budget. The set of constraints (10) makes sure that the sum of the projects' costs is in the range of the available budgets for the corresponding regions. The final set of constraints (11) makes sure that the binary variables x i can only have values of 0 or 1.
We should note that the problem definition is over the space defined by the binary vectors whose size is 2 n . Then the solution algorithms must search across this space to find the Pareto optimal solutions. On the other hand, given that the well-known NP-hard Knapsack problem can be easily reduced to MOPOP, the latter is also NP-hard [21].

Strategy to Generate the Fuzzy Trapezoidal Instances
This work uses instances initially designed for the POP with interval parameters, where the fuzzy representation of the parameters of the problem uses fuzzy interval type numbers (for example, the interval [76,800, 83,200]) [10]. Fixing the values of α, β to 0.5, and adding them to any interval in the original POP's instances allowed the creation of MOPOP's instances with Trapezoidal Fuzzy Parameters. Following this way, an interval value such as [76800, 83200] would be seen as [76800, 83200, 0.5, 0.5] in the new set of instances.
To create a random fuzzy interval type instance the following real parameters are considered: budget (B), number of objectives (m), projects (p), areas (a) and regions (r), and ranges of costs (c 1 , c 2 ), and objectives (m 1 , m 2 ). Then to generate a fuzzy interval instance the following interval type values must be determined: [B, B ] ← Budget as interval [a i , a i ] ← Limits of each area I = 1, 2, . . . , a [r i , r i ] ← Limits of each region r = 1, 2, . . . , r [b ij , b ij ] ← Benefit from the objective I = 1, 2, ..., m and for each project j = 1, 2, . . . , p {C i , A i , R i } ← Real values of the cost, area and region for each project i = 1, 2, . . . , p.
Implementing MOPOP's instances generator combines the previous parameters along with Equations (12)- (24) to create random instances [10]. a i = a l + Random (a l − a l ) for i = 1,2, . . . ,a (15) a i = a u + Random (a u − a l ) for i = 1,2, . . . ,a (16) r i = r l + Random (r l − r l ) for i = 1,2, . . . ,r r I = r u + Random (r u − r l ) for i = 1,2, . . . ,r (20) b ij = 1.1*o for i = 1, 2, . . . , p and i = 1, 2, . . . , m The interval instances, built with the instances generator, have names under the following format ompn_idI, where m is the number of objectives the instance has, n is the number of projects, id is a consecutive number, and I indicate that the instance is of interval type. An example of this would be the instance o2p100_1I, meaning that it is the instance number 1 with 2 and 100 projects.
The Algorithm 1 details the structure of a fuzzy interval type instance. with a = 0.5 and b = 0.5. The Algorithm 2 shows the result of converting the fuzzy interval type instance o2p25_0I to the fuzzy trapezoidal instance o2p25_0T.

Evaluating the Solutions and Verifying the Feasibility
This section describes how to calculate the objective values of a solution and how to determine its feasibility. To explain this process, let F the trapezoidal fuzzy numbers set, and R the set of real numbers. Now it is described how to apply the map δ : F → R such that δ(A) = P(A). The map associates the GMI value to each trapezoidal fuzzy number. A remarkable property of this map is that if X ⊂ F n , then δ(X) ⊂ R n , hence, the computation of a vector solution for a MOPOP's instance with two objectives is transformed into a vector of two trapezoidal fuzzy numbers, which in turn is transformed into a vector of two real numbers. As this process is consistently applied to all the solutions, the algorithms will be performed considering that the binary vector objectives space is the real vector space. The transformation must also be applied to all the trapezoidal fuzzy numbers in the constraints to validate the solutions' feasibility in the search space process. Equations (25)-(30) shows how evaluate the solution and verify the feasibility. where Subject to the following constraints: x i {0, 1} for all. i = 1, 2, . . . . . . , n An additional benefit is that this mapping transforms the approximated Pareto front in a set of real vectors. In such a case, standard commonly used metrics can be applied to evaluate the performance of the algorithms.
Example: Consider the following simplified instance: The objectives z 1 and z 2 are the benefits generated by the projects selected in the binary vector x. The constraint verifies that the cost of that project is not higher than the available budget (C).
Given the solution x = [0, 1,0], then the fuzzy trapezoidal values of the two objectives are the following: Now the GMI is used to compare the fuzzy trapezoidal numbers. For a trapezoidal fuzzy number A = (a, b, α, β), the GMI is: As P([10, 13, 0.2, 0.5]) = 11.55 ≤ P([3,20,1,5]) = 12.166, solution x is feasible. Notice that this process was done in the fuzzy trapezoidal numbers space; only at the end the GMI is used to verify the constraint. To perform the process in the real space, the two fuzzy objectives and the fuzzy costs in the constraint are transformed into real numbers using the GMI. The evaluation of the solution is as follows: Hence, the solution x is feasible given that 11.55 ≤ 12.166. The algorithms proposed in this work use the evaluation and feasibility verification procedures described in this section. The algorithms must call such methods on every new solution generated by them.

Steady-State T-NSGA-II Algorithm
This section presents the design of all the components included in the definition of the proposed algorithm. This is an adaptation of the classic Deb algorithm NSGA-II [22] modified to work with the trapezoidal fuzzy numbers. As all the algorithms proposed in this work, T-NSGA-II updates the population, applying in each generation the steadystate approach to include in the population only one of the generated individuals. In generational algorithms, the new set of offsprings are combined with the parents to create individuals' next generation; the input to the algorithm is a MOPOP's instance. The output is an approximate Pareto front for the instance.

Representation of the Solutions
A MOPOP's solution is represented by binary vector S = {0, 1} n , where n is the number of projects. This vector is a portfolio, and each value s i = 1 represents the inclusion of project i in the portfolio. The first element in the vector is s 0 , and the last is s n-1 . Figure 1 shows an example of this representation. Hence, the solution x is feasible given that 11.55 ≤ 12.166. The algorithms proposed in this work use the evaluation and feasibility verification procedures described in this section. The algorithms must call such methods on every new solution generated by them.

Steady-State T-NSGA-II Algorithm
This section presents the design of all the components included in the definition of the proposed algorithm. This is an adaptation of the classic Deb algorithm NSGA-II [22] modified to work with the trapezoidal fuzzy numbers. As all the algorithms proposed in this work, T-NSGA-II updates the population, applying in each generation the steadystate approach to include in the population only one of the generated individuals. In generational algorithms, the new set of offsprings are combined with the parents to create individuals' next generation; the input to the algorithm is a MOPOP's instance. The output is an approximate Pareto front for the instance.

Representation of the Solutions
A MOPOP's solution is represented by binary vector = {0,1} , where n is the number of projects. This vector is a portfolio, and each value si = 1 represents the inclusion of project i in the portfolio. The first element in the vector is s0, and the last is sn-1. Figure 1 shows an example of this representation.

One-Point Crossover Operator
The one-point crossover operator generates two offsprings from two parents [23]. The process first defines a random cutting point cp in the range [0, n -1]. After this, it split each parent vector into left and right sections, where for parent i, the lefti contains its values {s0, …, scp}, and the righti contains its values {scp+1, …, sn-1}. Finally, it mixes the split sections to generate two new offsprings h1, h2, where h1 uses left1 and right2, and h2 uses left2 and right1. The parents are chosen at random. The steady-state approach only utilizes the first offspring h1. The number of crossovers that are done is a defined parameter. Figure 2 shows an example of this operator.

One-Point Crossover Operator
The one-point crossover operator generates two offsprings from two parents [23]. The process first defines a random cutting point cp in the range [0, n -1]. After this, it split each parent vector into left and right sections, where for parent i, the left i contains its values {s 0 , . . . , s cp }, and the right i contains its values {s cp+1 , . . . , s n-1 }. Finally, it mixes the split sections to generate two new offsprings h 1 , h 2 , where h 1 uses left 1 and right 2 , and h 2 uses left 2 and right 1 . The parents are chosen at random. The steady-state approach only utilizes the first offspring h 1 . The number of crossovers that are done is a defined parameter. Figure 2 shows an example of this operator.

Uniform Mutation Operator
The uniform mutation operator generates a new solution for the mutation population from given a solution vector S = {s0, s1, …, sn-1 } [24]. The process generates for each index i, for 0 ≤ i ≤ n -1, a random number u in the range [0, 1], and if u < mut then the value of si changes from 1 to 0 or vice versa, otherwise the value si remains intact. The parameter mut is the mutation probability used by the operator. Figure 3 shows an example of the use of this mutation.

Uniform Mutation Operator
The uniform mutation operator generates a new solution for the mutation population from given a solution vector S = {s 0 , s 1 , . . . , s n-1 } [24]. The process generates for each index i, for 0 ≤ i ≤ n − 1, a random number u in the range [0, 1], and if u < mut then the value of s i changes from 1 to 0 or vice versa, otherwise the value s i remains intact. The parameter mut is the mutation probability used by the operator. Figure 3 shows an example of the use of this mutation.

Uniform Mutation Operator
The uniform mutation operator generates a new solution for the mutation population from given a solution vector S = {s0, s1, …, sn-1 } [24]. The process generates for each index i, for 0 ≤ i ≤ n -1, a random number u in the range [0, 1], and if u < mut then the value of si changes from 1 to 0 or vice versa, otherwise the value si remains intact. The parameter mut is the mutation probability used by the operator. Figure 3 shows an example of the use of this mutation. Another parameter of the operator is the number of new mutated solutions that must be generated. Usually, the solutions that undergo this process come from the crossover operator's results; otherwise are randomly chosen.

Initial Population
A predefined number of randomly generated solutions are created to have an initial population. When a new random solution is generated, the objectives vector for the solution is determined and its feasibility is verified.

Population Sorting
This process consists of sorting the solutions of the population, and it is composed of two phases: (1) the elitist phase, which keeps the best solutions; and (2) the diversification phase, which ensures that there are solutions different enough to avoid local optima in the search process of the algorithm. The elitist phase is also known as non-dominated Another parameter of the operator is the number of new mutated solutions that must be generated. Usually, the solutions that undergo this process come from the crossover operator's results; otherwise are randomly chosen.

Initial Population
A predefined number of randomly generated solutions are created to have an initial population. When a new random solution is generated, the objectives vector for the solution is determined and its feasibility is verified.

Population Sorting
This process consists of sorting the solutions of the population, and it is composed of two phases: (1) the elitist phase, which keeps the best solutions; and (2) the diversification phase, which ensures that there are solutions different enough to avoid local optima in the search process of the algorithm. The elitist phase is also known as non-dominated sorting. It consists of separating the population in fronts or sets of non-dominated solutions, making sure that the best solutions are always on the first front. The diversification phase sorts the solutions of a front according to the Crowding Distance indicator. The solutions in the best fronts are included in the population, and when a front cannot be completely inserted, the solutions with the worst crowding distances are discarded. Figure 4 shows both phases.

Non-Dominated Sorting
This process has two parts, and works on a given population. The first part constructs the first front with the set of non-dominated solutions identified from the comparison of vectors of objective values among all the population' solutions. A solution is nondominated if its vector of objective values is not dominated by any other. Note that the Pareto dominance uses real value vectors in its definition.
The second part builds the remaining fronts one by one. Each new front integrates those solutions that are only dominated by solutions in previously built fronts. The process repeats until no more fronts can be made.
sorting. It consists of separating the population in fronts or sets of non-dominated solutions, making sure that the best solutions are always on the first front. The diversification phase sorts the solutions of a front according to the Crowding Distance indicator. The solutions in the best fronts are included in the population, and when a front cannot be completely inserted, the solutions with the worst crowding distances are discarded.

Non-Dominated Sorting
This process has two parts, and works on a given population. The first part constructs the first front with the set of non-dominated solutions identified from the comparison of vectors of objective values among all the population' solutions. A solution is non-dominated if its vector of objective values is not dominated by any other. Note that the Pareto dominance uses real value vectors in its definition.
The second part builds the remaining fronts one by one. Each new front integrates those solutions that are only dominated by solutions in previously built fronts. The process repeats until no more fronts can be made.

Calculating the Crowding Distance
According to [22], this process orders the solutions in a front by their Crowding Distance (CD). The distance is a measure of the separation of the solutions, and it is relative to the normalized value of the objectives. The CD identify the solutions with extreme values on the objectives and put it first on the front. After that, the solutions order are according to their accumulated degree of separation per objective, the greatest the separation the better. For each objective, the CD computes the degree of separation using the ordered array of objective values resulting from the front; the solutions with the highest and smallest objective values will have a specific Crowding Distance value d equal to infinite (∞), while the remaining solutions will be calculated by the following formula:

Calculating the Crowding Distance
According to [22], this process orders the solutions in a front by their Crowding Distance (CD). The distance is a measure of the separation of the solutions, and it is relative to the normalized value of the objectives. The CD identify the solutions with extreme values on the objectives and put it first on the front. After that, the solutions order are according to their accumulated degree of separation per objective, the greatest the separation the better. For each objective, the CD computes the degree of separation using the ordered array of objective values resulting from the front; the solutions with the highest and smallest objective values will have a specific Crowding Distance value d equal to infinite (∞), while the remaining solutions will be calculated by the following formula: where d is the Crowding Distance, I is the solution position in the whole population in general, j is the solution position after the ordering by objective m within the front, f is the objective value and m is the current objective. The accumulation of Crowding Distance value d of all the objectives results in the final value of CD for each solution I.

Calculating the Spatial Spread Deviation (SSD)
The Spatial Spread Deviation (SSD) is a density estimator used to rearrange the solutions in a front, so the spread is not by a wide margin [25]. The method calculates for each solution the SSD value using a matrix of normalized distances between the solutions in the approximated front. The solutions are sorted from the lowest to highest SSD value in order to punish solutions according to their standard deviation and their proximity to their closest k-neighbors. The next three equations show how to calculate the SSD values, in the process i is the solution in the front for which the SSD is calculated, and j take values over all the solutions in the front except i.
where D(i, j) is the distance from solution i to solution j. D max is the biggest distance between all the solutions and D min is the closest distance between all the solutions. K is the number of k neighbors closest to solution i. SSD 0 is the initial value of SSD, which is -INF if the solution is at one of the ends of the front when the normalized values of the graded mean integration of the objective values are calculated.

T-FAME Algorithm
This section presents the design of all the components of the T-FAME algorithm. The algorithm adapts the FAME algorithm to work with the trapezoidal fuzzy numbers [25]. The input to the algorithm is an instance of MOPOP. The output is the approximate Pareto front for that instance. T-FAME updates the population, applying the steady-state approach to include in the population only one of the generated individuals. The following algorithm components are the same described in Section 4: the structure used to represent the solutions, the evaluation of a solution, the construction of the initial population, the sorting of the population, the non-dominated sorting process, and the density SSD estimator. The components described in this section are those not included in the previous description or with significant differences, such as the fuzzy controller, the additional genetic operators, and the structure used to store the approximated Pareto front.

Fuzzy Controller
This section introduces an intelligent mechanism that allows an MOEA to apply different recombination operators at different search process stages. The use of different operators is dynamically adjusted according to their contribution to the search in the past.
Intuitively, the idea is to favor operators generating higher quality solutions over others. For this purpose, the fuzzy controller dynamically tunes the probability selection of the available recombination operators [25].
The fuzzy controller uses a Mamdani-Type Fuzzy Inference System (FIS) [26] to compute the probability of applying the different operators. Fuzzy sets defined by membership functions represent the linguistic values of the model's input and output variables. Regarding the inference, we use the approach originally proposed by Mamdani based on the "max min" composition: using the minimum operator for implication and maximum operator for aggregation. The aggregation of the consequents from the rules are combined into a single fuzzy set (output), to be defuzzified (mapped to a real value). A widely used defuzzification method is the centroid calculation, which returns the area's center under the curve. We use triangular-shaped membership functions in all inputs and outputs, Once the fuzzification of the inputs is done, the next process is to evaluate the antecedents of the rules R 1 and R 2 , determining the following values: In the rule evaluation, the min operator is associated with the logic operator and, and the max operator is associated to the logic operator or. Now the membership functions of the consequents of the rules must be determined. For each rule an operator of implication is applied to the antecedent value obtained in the previous process and to the consequent of the rule, to determine the membership function of the conclusion of the rule. The min operator is used to implement the implication logic operator, which truncates the membership function of the rule's consequent. For example, the truncated membership functions of the consequents are the following: Now the truncated membership functions are integrated using an aggregation operator to create a new membership function, which is the controller's fuzzy output. The aggregation operators that are frequently used are max and sum.
For the example, the max operator is used to determine the aggregated membership function, which is the following: Finally, the defuzzification of the fuzzy output obtained is done. In this step a real number is associated to the aggregated membership function, which is the output of the inference process. In the previous example, the center of the area under the curve of the aggregated membership function is used to defuzzify the output of the controller as following: Figure 5 graphically shows the fuzzy inference process for the example described.  All of the controller rules are of the type: Antecedent AND Antecedent then Consequent. The fuzzy rules were designed to have soft changes in the input variables (Stagnation and UseOp), to avoid abrupt changes in the output variable (ProbOp). The configuration was manually done by observing the surface that these three variables generated [25]. Table 2 shows the rules of the fuzzy controller.  All of the controller rules are of the type: Antecedent AND Antecedent then Consequent. The fuzzy rules were designed to have soft changes in the input variables (Stagnation and UseOp), to avoid abrupt changes in the output variable (ProbOp). The configuration was manually done by observing the surface that these three variables generated [25]. Table 2 shows the rules of the fuzzy controller. The Algorithm 4 shows the structure of the fuzzy controller used in the fuzzy controller implementation with the Java Library Fuzzy Lite 6.0. The line numbers are those that appear in the T-FAME algorithm pseudocode included in Section 6.4. Notice that in lines 37 and 38, the algorithm uses the fuzzy controller to update all the available recombination genetic operators' selection probability.
The Stagnation value is shared for all the operators, and it is an indicator of the evolution of the search in the current time window. This is a normalized value that is increased by 1.0/Window each time the generated solution cannot enter the set where the non-dominated solutions are kept and reset when the time window is over. UseOp[i] is a normalized value that is increased by 1.0/Window every time the operator i is used.

Additional Genetic Operators
Four operators are used in T-FAME to create new solutions: One-point crossover, Uniform Mutation, Fixed Mutation, and Differential Evolution. Two of these operators (One-point crossover and Uniform Mutation) are the same ones that are used on T-NSGA-II, and they are already described in the previous section.
Differential Evolution: This method was proposed by Rainer [27], and its implementation was based on [28]. It uses the four parents obtained with the tournament method. The first part of the process consists of creating a new solution called Candidate using Parent 1, Parent 2, and Parent 3, this solution is obtained by doing a binary addition of the parents. Figure 6 shows an example of how this operator works.
Once the Candidate is calculated, a binary crossover operator is done between the candidate and Parent 4 to create a new solution called Son, this binary crossover operator is different from the one-point crossover operator described previously, and it uses a parameter called crossover percentage (CP). The binary crossover operator consists of the following: For each array index, a random number between 0 and 1 is generated, if that number has a lesser value than CP, then that index receives the value of the Candidate, if this is not the case, then that index receives the value of Parent 4.
(One-point crossover and Uniform Mutation) are the same ones that are used on T-NSGA-II, and they are already described in the previous section.
Differential Evolution: This method was proposed by Rainer [27], and its implementation was based on [28]. It uses the four parents obtained with the tournament method. The first part of the process consists of creating a new solution called Candidate using Parent 1, Parent 2, and Parent 3, this solution is obtained by doing a binary addition of the parents. Figure 6 shows an example of how this operator works. Once the Candidate is calculated, a binary crossover operator is done between the candidate and Parent 4 to create a new solution called Son, this binary crossover operator is different from the one-point crossover operator described previously, and it uses a parameter called crossover percentage (CP). The binary crossover operator consists of the following: For each array index, a random number between 0 and 1 is generated, if that number has a lesser value than CP, then that index receives the value of the Candidate, if this is not the case, then that index receives the value of Parent 4.
Once the new solution Son is completed, a dominance test is done between Son and Parent 4, if the objective values of Parent 4 dominate the objective values of Son, then Parent 4 proceeds to be the new solution, but if this is not the case, then Son proceeds to be the new solution.
Fixed Mutation: This method is very similar to the uniform mutation operator that was described previously. The main difference lies in the fact that the whole process is done in a loop until n mutations are made, where n is a parameter previously defined. This operator also makes sure that no element in the solution is changed twice or more times, this is done by using a fixed array to keep track of the changed elements in the solution. Figure 7 shows an example of the Fixed Mutation operator. Fixed Mutation: This method is very similar to the uniform mutation operator that was described previously. The main difference lies in the fact that the whole process is done in a loop until n mutations are made, where n is a parameter previously defined. This operator also makes sure that no element in the solution is changed twice or more times, this is done by using a fixed array to keep track of the changed elements in the solution.

Used Structures to Store the Population and the Approximated Pareto Front
The algorithm uses the structure pop to maintain a solutions population, which contains the following information for each solution i:

Used Structures to Store the Population and the Approximated Pareto Front
The algorithm uses the structure pop to maintain a solutions population, which contains the following information for each solution i: vector binary associated to the solution i.

T-FAME Algorithm Pseudocode
This section presents the pseudocode for the algorithm T-FAME in Algorithm 5. if (RandomDouble(0,1) ≤ β) then **The parent is chosen from Front
Else **The parent is chosen from pop 10.

Experimental Results
Two experiments were done in order to evaluate the performance of the proposed algorithms. The tested steady-state algorithms were T-NSGA-II-CD, T-NSGA-II-SSD, and T-FAME. The first experiment was done to make sure the algorithms were implemented correctly, while the second experiment was done to compare the performance between them using performance metrics.
The software and hardware platforms that were used for these experiments include Intel Core i5 1.6GHz processor, RAM 4GB, and IntelliJ IDEA CE IDE.

Performance Metrics Used
In order to measure the performance of each algorithm, two metrics were used: hypervolume [28] and generalized spread [29].
Hypervolume is the n-dimensional solution space volume that is dominated by the solutions in the reference set. If this space is big, then that means that the set is close to the Pareto Front. It is desirable for the indicator to have large values. Generalized Spread calculates the average of the distances of the points in the reference set to their closest neighbor. If this indicator has small values, then that means the solutions in the reference set are well distributed.

Experimental Setup
In order to configure the algorithms used in this work, the parameter values reported in the state-of-the-art were considered. The parameter value for the maximum number of evaluations was determined after a preliminary experimental phase. The comparison of all the algorithms, under the same operation conditions, utilizes a steady-state approach, using the dominant son. Tables 3 and 4 show the values of the parameters used in the algorithms. The configuration of algorithm T-NSGA-II-SSD is the same one as T-NSGA-II-CD, however, it uses Spatial Spread Deviation instead of Crowding Distance as its density estimator. Table 3. T-NSGA-II-SSD parameters.

Parameter Value
Evaluation of the objective function 5000 Population Size 50 Crossover population % 70 Mutation population % 40 Mutation % 5 Table 4. T-FAME parameters.

Parameter Value
Evaluation

Experiment 1. Validating the Implemented Algorithms
For this experiment, an instance named o2p25_rand was used, this instance was originally created for POP with intervals, which was converted in a trapezoidal fuzzy instance by adding two parameters to the intervals. The optimum Pareto Front was obtained using an exhaustive algorithm, and approximate fronts were obtained with T-NSGA-II-CD, T-NSGA-II-SSD, and T-FAME algorithms. All algorithms solve the MOPOP with Fuzzy Parameters and use a steady-state election mechanism, creating one solution from the genetic operators' application. This adaptation from FAME has an advantage over algorithms using the classic generational approach in genetic algorithms.
The purpose of this experiment is to validate the correct operation of the implemented algorithms in the project. In the experiment, the fronts are generated, and they are compared to the optimum front, in order to determine if the algorithms are generating similar fronts. All the fronts that were generated are shown in Table 5. Each front is shown in two columns that contain the values of the two objectives that were originally Trapezoidal Fuzzy numbers, but they were converted into real numbers with the transformation based on GMI. The graph the fronts uses the GMI values obtained from the objectives. Table 5. Generated fronts of the algorithms with instance o2p25_rand.

Pareto Optimal Front
T-NSGA-II-CD T-NSGA-II-SSD T-FAME It is worth nothing that, in Figure 8, the approximated fronts are relatively close and below the optimum front. Also, observe that the T-NSGA-II-SSD and T-FAME algorithms managed to reach some optimum solutions. Finally, note that the T-FAME algorithm has a good distribution between its solutions. It is worth nothing that, in Figure 8, the approximated fronts are relatively close and below the optimum front. Also, observe that the T-NSGA-II-SSD and T-FAME algorithms managed to reach some optimum solutions. Finally, note that the T-FAME algorithm has a good distribution between its solutions.

Experiment 2. Evaluating the Performance of the Algorithms with Instances of 25 Projects
This experiment evaluates the performances of algorithms T-NSGA-II-CD, T-NSGA-II-SSD, and T-FAME, and utilizes 13 instances with 2 objectives and 25 projects. In order to compare the performance between the three algorithms, each algorithm was executed 30 times per instance. The performance metrics used were hypervolume and generalized

Experiment 2. Evaluating the Performance of the Algorithms with Instances of 25 Projects
This experiment evaluates the performances of algorithms T-NSGA-II-CD, T-NSGA-II-SSD, and T-FAME, and utilizes 13 instances with 2 objectives and 25 projects. In order to compare the performance between the three algorithms, each algorithm was executed 30 times per instance. The performance metrics used were hypervolume and generalized spread. For each instance, the reference set contains the non-dominated solutions obtained from the combination of the 30 generated fronts. The computation of the metrics uses the reference set as an approximation to the optimum Pareto Front. The computation of the median value and interquartile ranges uses the metric values of the 30 instances sorted in ascending order. With the sorted array, the median value was the average of the metric values from positions 15 and 16. At the same time, the interquartile ranges correspond to those in positions 23 and 8, corresponding to the 75% and 25% of the metrics values, respectively. The median value and the interquartile ranges are used instead of the average and the standard deviation because they are less sensitive to extreme values. The experiment performs a hypothesis test to validate the obtained results. The hypothesis was proven using the parametric t student test on those data sets that passed the normality and homoscedasticity tests and using the non-parametric Wilcoxon signed-rank test on those that do not. Both tests apply a confidence level of 95%, pairing T-FAME with each of the other two algorithms. Tables 6-9 show the results of the normality and homoscedasticity tests done for all the instances used in this work (25 and 100 projects) and the metrics of hypervolume and generalized spread. Tables 6 and 8 show in the last column pairs (i,j), which indicate that the comparison of T-NSGA-II-CD and T-FAME uses test i, and the comparison T-NSGA-II-SSD and T-FAME uses test j. The values t and W in (i, j) stand for t student test and Wilcoxon test. This work tests each instance separately. Table 6. Hypervolume normality test, the null hypothesis is that the samples follow a normal distribution which is accepted (a) when p-value < 0.05 and rejected (r) otherwise.   Table 8. Generalized Spread normality test, the null hypothesis is that the samples follow a normal distribution which is accepted (a) when p-value < 0.05 and rejected (r) otherwise.   Table 10 shows the performance results with the hypervolume metric, and Table 11 shows the results with the generalized spread metric. For the hypervolume metric, the algorithm with the largest value is considered to be the one with the best performance. For the generalized spread metric, the best algorithm is considered to be the one with the smallest value. The table's cells show the median value of the metric (M) and the interquartile range (IRQ) in the following format: M IRQ . In the result tables, for each instance the best and second-best values are marked with solid or light black, respectively. In order to indicate if the observed differences in the performance of the algorithms are significant or not, for each algorithm the symbol indicates that the performance of T-FAME is significantly better that the algorithm which it is being compared. The symbol indicates the opposite, and the symbol = indicates that the difference is not significant. These symbols are marked with an asterisk when the t student test was applied. To confirm the results obtained with the paired tests, a global evaluation is done with the three algorithms. This evaluation was done by applying a Friedman test with 95% confidence.  Table 11. Results with the generalized spread metric.

Generalized Spread
Instance The information presented in Table 10 shows that T-NSGA-II-CD stands out as the algorithm with the best performance in 12 of 13 cases. The results on Table 11 shows that T-NSGA-II-SSD positions itself as the best algorithm in 10 of 13 cases and T-FAME in 8 of 13 cases. It can also be observed that these differences are significant in all cases, this is due to the fact that when the differences are not significant between the best and second-best algorithms, then that means the algorithms are considered tied. Table 12 confirms the results observed with the t student and Wilcoxon tests. As a result of applying the Friedman test with the three algorithms, the ones with the lowest rank for the hypervolume and generalized spread metrics are T-NSGA-II-CD and T-NSGA-II-SSD, respectively.  T-NSGA-II-CD  14  T-NSGA-II-SSD  19  T-NSGA-II-SSD  25  T-FAME  20  T-FAME  39  T-NSGA-II-CD  39 6.5. Experiment 3. Evaluation of the Algorithm' Perfomances Using Instances with 100 Projects

T-NSGA-II-CD T-NSGA-II-SSD T-FAME
As indicated previously, the previous experiment was done with instances with 25 projects, for which the algorithms had to navigate in a space of binary vectors of length 25. In that case the size of the solution space was of 2 25 . For this experiment, 9 instances of 2 objectives and 100 projects were used, these instances represented a greater complexity for the algorithms because the solution space increased to 2 100 . The experiment conditions were just as in the previous one, using the same metrics but in a scenario of greater complexity scenario. For each instance, the reference set contains the non-dominated solutions obtained from the combination of the 30 generated fronts. The computation of the metrics uses the reference set as an approximation to the optimum Pareto Front. The computation of the median value and interquartile ranges uses the metric values of the 30 instances sorted in ascending order. With the sorted array, the median value was the average of the metric values from positions 15 and 16. At the same time, the interquartile ranges correspond to those in positions 23 and 8, corresponding to the 75% and 25% of the metrics values, respectively. The experiment performs a hypothesis test to validate the obtained results. The hypothesis was proven using the parametric t student test on those data sets that passed the normality and homoscedasticity tests and using the nonparametric Wilcoxon signed-rank test on those that do not. Both tests apply a confidence level of 95%, pairing T-FAME with each of the other two algorithms. Tables 6-9 shows the results of the normality homoscedasticity tests done for all the instances used in this work (25 and 100 projects) and the metrics of hypervolume and generalized spread. Table 13 shows the results with the hypervolume metric and Table 14 shows the results with the generalized spread metric. For the hypervolume metric, the algorithm with the largest value is considered to be the one with the best performance. For the generalized spread metric, the best algorithm is considered to be the one with the smallest value. The table cells show the median value of the metric (M) and the interquartile range (IRQ) in the following format: M IRq . In the result tables, for each instance the best and second best values are marked with solid or light black, respectively. In order to indicate if the observed differences in the performance of the algorithms are significant or not, for each algorithm the symbol indicates that the performance of T-FAME is significantly better that the algorithm which it is being compared. The symbol indicates the opposite, and the symbol = indicates that the difference is not significant. These symbols are marked with an asterisk where the t student test was applied. To confirm the results obtained with the paired tests, a global evaluation is done with the three algorithms. This evaluation was done by applying a Friedman test with 95% confidence.  Table 14. Results with the generalized spread metric.

Generalized Spread
Instance The information presented in Table 13 shows T-FAME stands out as the algorithm with the best performance in 7 of 9 cases and T-NSGA-II-SSD in 5 of 9 cases. The results on Table 14 show that T-FAME stands out as the best algorithm in 6 of 9 cases and T-NSGA-II-SSD in 4 of 9 cases. These differences are significant in all cases, this is due to the fact that when the differences are not significant between the best and second-best algorithms, then that means the algorithms are considered tied. Table 15 confirms the results observed with the t student and Wilcoxon tests. As a result of applying the Friedman test with the three algorithms, the one that has the lowest rank for both metrics is T-FAME.

Conclusions and Future Work
This work approaches the Multi-Objective Portfolio Optimization Problem with Trapezoidal Fuzzy Parameters. To the best of our knowledge, there are no reports of this variant of the problem. This work, for the first time, presents a mathematical model of the problem, and, additionally, contributes with a solution algorithm using the Fuzzy Adaptive Multiobjective Evolutionary (FAME) methodology and two novel steady state algorithms that apply the Non-Dominated Genetic Algorithm (NSGA-II) methodology to solve this variant of the problem. Traditionally, these kinds of algorithms use the Crowding Distance density estimator, so this work proposes substituting this estimator for the Spatial Spread Deviation to improve the distribution of the solutions in the approximated Pareto fronts. This work contributes with a defuzzification process that permits measurements on the algorithms' performances using commonly used real metrics. The computational experiments use a set of problem instances with 25 and 100 projects and hypervolume and generalized spread metrics. The results with the challenging instances of 100 projects show that the algorithm T-FAME has the evaluated algorithms' best performance. Three hypothesis tests supported these results, and this is encouraging because they confirm the feasibility of the proposed solution approach.
The main open works identified in this research are to develop algorithms for solving the problem with many objectives, preferences, and dynamic variants. Currently, we are working to change the fuzzy controller selector for a selector based on a reinforcement learning agent.

Conflicts of Interest:
The authors declare no conflict of interest.