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.


Plagiarism Detection Methods and Tools : An Overview

Plagiarism Detection Systems play an important role in revealing instances of a plagiarism act, especially in the educational sector with scientific documents and papers. The idea of plagiarism is that when any content is copied without permission or citation from the author. To detect such activities, it is necessary to have extensive information about plagiarism forms and classes. Thanks to the developed tools and methods it is possible to reveal many types of plagiarism. The development of the Information and Communication Technologies (ICT) and the availability of the online scientific documents lead to the ease of access to these documents. With the availability of many software text editors, plagiarism detections becomes a critical issue. A large number of scientific papers have already investigated in plagiarism detection, and common types of plagiarism detection datasets are being used for recognition systems, WordNet and PAN Datasets have been used since 2009. The researchers have defined the operation of verbatim plagiarism detection as a simple type of copy and paste. Then they have shed the lights on intelligent plagiarism where this process became more difficult to reveal because it may include manipulation of original text, adoption of other researchers' ideas, and translation to other languages, which will be more challenging to handle. Other researchers have expressed that the ways of plagiarism may overshadow the scientific text by replacing, removing, or inserting words, along with shuffling or modifying the original papers. This paper gives an overall definition of plagiarism and works through different papers for the most known types of plagiarism methods and tools.

Introduction
Due to the rapid advancement of the computer and network technologies, such as the Internet that enables anyone to access online contents anytime and from anywhere, academic integrity in the academic community is becoming a highly sensitive issue, especially among universities and research institutions. Plagiarism, on the other hand, is defined as a kind of academic dishonest behavior that will damage academic integrity [1]. Thus, it is needed to be resisted determinedly. However, plagiarism is not only an academic issue, but it extends to almost all industries. Occasionally, plagiarism occurs accidentally but most of the time it is the outcome of a conscious process [2]. The best definition of plagiarism might be that it is "the unacknowledged copying of documents or programs" [3]. To overcome the problem of plagiarism, large number researchers have worked on detecting plagiarism since the past decades through software detection methods [4,5]. Plagiarism was originally detected manually (by hand) or by resembling previously consulted content. Today, the great number of the available online documents make it harder to detect plagiarism manually. Therefore, there is an urgent need to produce automatic plagiarism detectors [5]. There are two main types of plagiarism, namely the verbatim/literal and the intelligent plagiarism. Plagiarism detection methods are also classified into the internal detection method, where the document is analyzed for plagiarism alone, and the external detection method, where detection is made among a collection of documents. Verbatim/literal plagiarism describes the plagiarized content as the exact copying of the source content without altering or modifying the original content. While, in intelligent plagiarism, the main content is altered/modified by different ways. Intelligent plagiarism is more difficult to reveal and includes adoption of the ideas, translation to another language, and manipulations [6,7].  [6,8] This overview paper sheds the light on the description of plagiarism. In section two, plagiarism process will be reviewed, in section three, plagiarism classification and methods will be explained in details, in section four, plagiarism tools will be reviewed, in section five, the types of datasets used in plagiarism detection will be illustrated, in section six, a discussion about the reviewed works will be summarized, and finally in section seven, a conclusion will summarize the topic of plagiarism.

Plagiarism Process
To design and produce robust and no error Plagiarism Detection Process, four main stages are required. These stages described below [8,9]: 1. Collecting the content: It is the first stage, where the plagiarism detector collects the required content from the users through a search engine which acts as the interface between the users and the detector. 2. Analyzing the similarities: After collecting the content (scientific papers, assignments, and other softcopies), the detector runs an analyzing method to search for the similarities among the documents and reveal the original copy. 3. Confirming the copy: After the analyzing stage, a function for plagiarism conformation is required to reveal the plagiarized text from an original one. Sometimes, a degree of the plagiarized text is confirmed with this process. 4. Investigation: It is the final stage, which depends highly on the interference made by the user whenever a plagiarism is confirmed. It also relies on the expert of the user to distinguish between the really plagiarized documents and the cited ones.

Plagiarism Classification
Plagiarism can be divided into two basic categories, which are the monolingual and the crosslingual. Monolingual plagiarism works with most detectors. It is about homogeneous languages, as in the case of English language setting-English language setting. Cross-Lingual plagiarism works with heterogeneous languages; for example, English language setting-Chinese language setting, and this type is quiet rare [5] [10] [11]. In the next section plagiarism, types will be discussed in details.

Plagiarism Types
Plagiarism types appear in different works, documents, scientific papers, and research article. It can be classified as in the following ways [12]: (i) pretension of others work as your work, (ii) copying others' work without mentioning the credit or citation, (iii) whether citation was mentioned or not, calming someone's contribution as your own, (iv) refereeing to others work as yours by reconstructing their work, and (v) adding a misleading acknowledgments of others as your work. Textual plagiarism and Source Code plagiarism are the two main types of plagiarism and they will be reviewed in the following [13,14]

Textual Plagiarism
In researches and scientific fields, this type of plagiarism is the most common one, where the entire text or document is taken without referring to the author or mentioning a quotation. This type of plagiarism can be further divided into seven sub-classes, as in the following [13,14]: 1. Copy-paste plagiarism: This process refers to copying the original text without any acknowledgment about the authors or the original paper as if it was your work. 2. Paraphrasing plagiarism: It is classified into two categories: (i) simple paraphrasing, where the original text is presented into different way be replacing the words into similar ones with the same meaning and, (ii) Mosaic/Hybrid/Patchwork paraphrasing, where the text is a result of combining different contributions from different papers and presented differently without referring to the original citations of the works. 3. Metaphor plagiarism: presenting other ideas in better ways. 4. Idea plagiarism: the entire solution and ideas are stolen from others and claiming that it is an original research paper. 5. Recycled plagiarism: The authors here use their previous/old works and papers for a new publication. 6. 404 Error / Illegitimate Source plagiarism: when the citation of the works is invalid. 7. Re-tweet plagiarism: In this type, the citation is referred to but it is no difference between the original work and the author's work from the point of structure, grammar, and words.

Source-Code Plagiarism
This type appears typically in educational fields, where the programming code of a specific program written originally by someone is copied, adjusted, or reused by others partially or completely. It has the following four sub-classes [13,14]: 1. Manipulation plagiarism: where the source code is altered or modified by other developers by either deleting or inserting sub-codes to an original one without referring to the citation or acknowledgment. 2. Reordering structure plagiarism: where the syntax of source code is modified by functions or statements recording without referring to the original work. 3. No-change plagiarism: where the developers do not change anything in the code but add/remove spaces or comments as it was their work. 4. Language switching plagiarism: where the source code language is rewritten by other languages and declared as original code.

Plagiarism Detection
Many papers have searched for highly accurate plagiarism detection methods using different tools, but it was always challenging to find the perfect one, due to the rapid development of the technologies, software, and data mining tools. This development has become a double-edged weapon, as the methods of plagiarism have evolved; on the other hand, the methods of detecting this theft have developed in a response to the curbing of illegal methods of copying the original work of researchers [10,11] Plagiarism detection can be performed manually or by using an automated process. The automated process is very similar to natural language processing, visual identification, and bio-metric processes. All of these have a foundation of pattern recognition. Automated process does not give 100% accuracy. Thus, the manual checking is still needed.

Internal Plagiarism Detection
Thus type involves finding plagiarized passages within a document without access to the potential original text, also called intrinsic plagiarism detection.

External Plagiarism Detection
External plagiarism detection involves comparing suspicious plagiarized documents against potential original documents.  [10,15] In the next sections, different methods and tools for plagiarism detection will be illustrated.

Plagiarism Detection Methods
Many methods have been implemented by researchers to overcome plagiarism as it has grown to form a serious issue among the academic community; researchers have used different methods to overcome these activities [5]. Therefore, a comparative study depending on the attached sources about plagiarism detection, as viewed by researchers, is illustrated in Table-1. It is quite challenging to reveal a plagiarism act because it has many types of plagiarism, like cross-lingual syntax-based methods and/or dictionarybased methods. It requires extensive in-depth knowledge for multiple languages in more than one document [36,37]. 9. Grammar Semantics Hybrid Plagiarism Detection Method: A very effective and extensively implanted approach in the field of natural language processing (NLP) . It is very accurate in revealing copy-paste plagiarism or paraphrasing plagiarism. It provides a remedy for the limitations of the semantic-based method [38]. 10. Classification and Cluster-Based Methods: These are greatly helpful methods to retrieve the information during the process of searching in any plagiarized document. Also, the comparison time is reduced during the detection process when comparing these methods with other ones [39]. 11. Citation-Based Method: This is a novel method; it mainly belongs to semantic plagiarism detection methods for the usage of the semantics in the cited document. It looks for identical pair of documents based on the citation, because these techniques use semantics contained in the citation [40,41].

Plagiarism Detection Tools
A large number of tools have been developed and utilized to detect plagiarism [5]. Table 3 shows plagiarism detection tools according to their pros and cons, covering a period of 22 years from 1994 to 2020 [5]. it takes parts of the code at a time as an input and produces HTML pages as an output to analyze the similarities between a pair of documents.

Ithenticate [43] 1996
It is a text-document based plagiarism detection tool that is presented as a web page. It compares a number of documents with the original one without the need for installation on the end-user computer, but it is limited to 25,000 words per time.
JPlag [32] 1997 Similar to the previous ones, this type is an online source-code plagiarism tool. It takes a number of programming codes and selects the identical lines among them. It works with C, C++, and Java programming languages, with less than one minute to detect hundreds of

Tool
Year Characteristics code lines.
GPSP -Glatt Plagiarism Screening Program [44] 1999 Unlike the previous tools, it works off-line. It mixes different approaches and finds the similarities among the writing styles of differed authors. It reveals plagiarism by making the author goes through a fill-in-the-blank test. Then, it counts the correctly filled blanks and the time taken to finish the test. Finally, according to the results, it takes a decision about an act of plagiarism.

Turnitin [42][37] 2000
It is provided by iParadigms as a web based tool. The user is required to upload his/her required document online, then the document will be saved to the system's database. After that, the tool checks for plagiarism by creating a document fingerprint. It accepts nearly 15,000 institutions around the world, with more than 30 million users, for its flexibility and robustness. Therefore, it is considered as the best tool. Plagiarism Checker [45] 2006 This is a free and online tool, using search engine services to detect for students' plagiarism by checking if their documents have a similar copy of another online document.

Plagiarism Scanner [46]
2008 It is an effective tool that detects throughout almost all online resources, like libraries (Questia and ProQuest), online databases, websites, and search engines. When plagiarism is detected, it produces a full report including the rate, originality, and percentage of plagiarized materials.

PlagTracker [47] 2011
It is a well-known tool for all kinds of users (teachers, websites owners, and students)that accommodates a large number of academic resources in its database and produces a detailed report whenever a plagiarism is detected.

PlagScan [48] 2015
This tool provides multiple services to companies, universities, and schools, but it is not free and the users must have a paid account to register to this tool.

Exactus Like [49] 2016
This tool is a web-based online tool that works with different formats, like HTML, Microsoft Word, and Adobe PDFs. It detects moderately disguised borrowing (word/phrase reordering, substitution of some words with synonyms) by a deep parsing function.

Grammerly [50] 2016
This is a website and a mobile application service that offers a great opportunity to the individuals to correct their documents within a realtime manner and a friendly user interface. It works online; therefore it requires an internet connection.

Grammerly [50] 2018
It is an evolved version of the previous one, representing the premium type. It targets business industries, such as teams and companies. Users reported that Grammarly helps them more professionally.
DupliCheck er [51,52] 2020 This is an absolutely perfect method, available 24/7, and ready whenever the user needs it. It is one of the most effective and free plagiarism tools on the internet. The user only requires a search engine and a connection to the world wide web to access this tool. It enables the user to either copy-paste or upload the document to check for a plagiarism.

Datasets Used in Plagiarism Detection Systems
Two main datasets are used for plagiarism detection and are illustrated below.

WordNet Dataset
WordNet is a freely and publicly available dataset that contains large lexical English language words, such as nouns, verbs, adverbs, and adjectives. It contains over 155,287 words organized in 117,659 synsets. All these words are classified into groups of cognitive synonyms. WordNet not only links word forms (strings) but also specifies their meanings. As a result, words that are next to each other in the dataset are semantically disambiguated. Also, WordNet marks the semantic relationships among groups of words in the thesaurus that do not follow any clear pattern, except for similar meanings. Lexical-words that are represented by this dataset contain synonymy between words, like the words "large" and "wide", both having a relatively similar meaning. The phrase of a noun contains substitution definitions that contain one or more substitutions. Therefore, each formal meaning pair in WordNet is unique. The current version of this dataset contains not only English language words but also different languages, such as Italian and Spanish [53].

Plagiarism analysis, Authorship identification, and Near-duplicate detection (PAN)
Another well-known plagiarism detection dataset is PAN. It refers to plagiarism analysis, authorship identification, and near-duplicate detection of different types of plagiarism. It is a series of scientific events and shared tasks on digital text forensics and stylometry. Every year, an international conference and competition called PAN@CLEF are held to connect the most advanced publications about plagiarism detection techniques [54].

Discussion
In this part, an extensive review of the most serious and frequently used plagiarism techniques around the world is illustrated. Also, the most common challenges that are facing the development of effective and robust plagiarism detection systems are reviewed.

Comparison of Common Plagiarism Types
The detection systems are available in-hand and they are growing rapidly. Therefore, a comparison among the nowadays plagiarism techniques is required. This comparison will focus on two studies that show the most serious plagiarism techniques in the universities, schools, and higher education sector. This comparison is provided by conducting several scientific areas, such as those of medicine , natural sciences, engineering, and social sciences, from 40 different countries around the world, covering the period from 2013 through 2015. Tables 4 and 5 list the most commonly practiced plagiarism types [55] [56].

Challenging Factors Facing Plagiarism Detection
Among the previously published papers that have been focusing on plagiarism, many works made a dense search on plagiarism types and techniques. Most of the plagiarism detectors available today can do the following: (i) distinguish between plagiarism in source code and/or in-text documents, with or without citation, (ii) feature extraction of semantic and/or salient syntactic, and (iii) plagiarism detection for both cross-lingual and monolingual documents [5]. Despite the availability and efficiency of these types, they are not effective enough to reveal the unattended research challenges. With the technology age that we are living in, new algorithms are produced to solve many problems, such as plagiarism, which is a problematic issue that needs to be extensively solved, especially in the scientific community. Computer science approach can address these challenges relying on the ICT advancement, some of these challenges are highlighted in the following: (i) A proof for correctness and completeness of the scientific works, i.e. whether they are ate in text documents or written as source code, is not available yet, (ii) a highly accurate framework for plagiarism detection that can reveal text segment(s), for both intrinsic and extrinsic plagiarism detection, is missing, (iii) The development of pilgrim checking systems without the need for external references and with high accuracy is considered a very challenging task, and (iv) Providing a full system for scientific works repository that can combine the works of all authors and the references to their works in one place is a difficult manner [5].

Conclusions
In this paper, an extensive literature survey about plagiarism types, methods, classification, and tools was conducted. Text plagiarism, with its seven sub-types, and source-code plagiarism, with its four sub-types were highlighted. Then, plagiarism detection methods were illustrated and summarized over nearly more than twenty years. The newly developed tools are more advanced, most of which are working online using an internet connection and a web page, some of them are delivered freely and others require subscription payment. Next, the most known datasets implemented in plagiarism detection were reviewed and a table was prepared to discuss the methods to be adopted in plagiarism . We notice that each method has its strengths and weaknesses that depend on how it is described to support two important factors: time and accuracy. We also notice that there are two mechanisms of action, namely the parallel and the series mechanisms. The parallel mechanism provides higher accuracy and less time because it performs a scanning for all the contents of a dataset. For this reason, we can say that a good algorithm is the one that covers the required conditions in terms of time, accuracy, or both. Finally, a discussion about the most frequent plagiarism types was extensively provided and the most challenging steps during the implementation of plagiarism detection were investigated.


Forecasting COVID-19 confirmed cases, deaths and recoveries: Revisiting established time series modeling through novel applications for the USA and Italy

The novel coronavirus (COVID-19) is an emergent disease that initially had no historical data to guide scientists on predicting/ forecasting its global or national impact over time. The ability to predict the progress of this pandemic has been crucial for decision making aimed at fighting this pandemic and controlling its spread. In this work we considered four different statistical/time series models that are readily available from the ‘forecast’ package in R. We performed novel applications with these models, forecasting the number of infected cases (confirmed cases and similarly the number of deaths and recovery) along with the corresponding 90% prediction interval to estimate uncertainty around pointwise forecasts. Since the future may not repeat the past for this pandemic, no prediction model is certain. However, any prediction tool with acceptable prediction performance (or prediction error) could still be very useful for public-health planning to handle spread of the pandemic, and could policy decision-making and facilitate transition to normality. These four models were applied to publicly available data of the COVID-19 pandemic for both the USA and Italy. We observed that all models reasonably predicted the future numbers of confirmed cases, deaths, and recoveries of COVID-19. However, for the majority of the analyses, the time series model with autoregressive integrated moving average (ARIMA) and cubic smoothing spline models both had smaller prediction errors and narrower prediction intervals, compared to the Holt and Trigonometric Exponential smoothing state space model with Box-Cox transformation (TBATS) models. Therefore, the former two models were preferable to the latter models. Given similarities in performance of the models in the USA and Italy, the corresponding prediction tools can be applied to other countries grappling with the COVID-19 pandemic, and to any pandemics that can occur in future.

Introduction
Cases of severe respiratory illness began to be reported across the city of Wuhan in China in December 2019. These were caused by a new type of coronavirus, SARS-CoV2, alternatively known as COVID-19 [1]. The number of COVID-19 cases quickly escalated in mid-January, and the virus soon spread beyond China's borders. As of this study, it has reached 188 countries worldwide [2]. Despite broadening COVID-19 impacts, there remain key public health questions. How many people will be infected? How does the situation change day by day? Can we predict/forecast the future numbers of people infected with COVID-19 by using daily updated data on the pandemic's trajectory? These questions could be answered by forecasting the possible futures of this pandemic through statistical/time series models. However, widely used time series models and tools do not necessarily have high forecasting/prediction accuracy, especially for medical studies [3]. As the COVID-19 pandemic evolves and additional data are amassed [4], new insights emerge on the reasonableness of the prediction models developed and tested; however, predicting its future requires transparent reporting and multiple model assessments [3]. Despite these ongoing challenges, forecasting is still critical to better understand the current situation and progression of COVID-19, to prepare for the future of this ongoing pandemic. Using appropriate models and making accurate projections consistently may enhance public health decision-making; for example, intervention-induced changes in the spread of disease [5]. In this paper, we provide statistical forecasts for the confirmed cases of COVID-19 using four different, highly cited time series models, which we subsequently introduce, and compare their performance to analyze the trajectory of cases.
In this work, we propose examining the utility of these models for forecasting the number of new COVID-19 cases in select countries based on the daily reported data. By empirically comparing multiple models in terms of their forecasting accuracy, we intend to suggest an appropriate model, which can be readily used by society, organizations, or governments to assess near futures of this outbreak. We acknowledge that this is a difficult forecasting problem, since this pandemic continues and there are many factors that we cannot presently control. Our prediction accuracy will improve with time and more data. An appealing part of these models that we considered is that they are already implemented in an accessible R package and can be applied readily by analysts. Fitting each model only takes several seconds, and the models are conducive to daily updating as newly available data are accrued. Daily updating of these models is suggested, because the accuracy of the forecasting depends on the available data used for model development, producing predictions, and generating prediction intervals (PIs) [4].

Materials and methods
We considered four time series models for our case study. The first model was the Holt model [6,7], which employs double exponential smoothing; two smoothing components adjust for level and trend in the data. The second model was the Autoregressive Integrated Moving Average (ARIMA) model, which theory dictates is the most general class of models for forecasting a time series. Although the ARIMA model is useful and powerful in time series analysis, it is somewhat difficult to choose appropriate orders for its components. For our purposes aimed at utility, it was necessary to determine the orders automatically; thus, we used the 'auto.arima' function from the forecast package in R, which returns the best ARIMA model according to model fit statistics by searching over possible models within the order constraints provided. Third was the TBATS (Trigonometric Exponential smoothing state space model with Box-Cox transformation, ARMA errors, Trend and Seasonal component) model [8], which uses a combination of Fourier terms with an exponential smoothing state space model and a Box-Cox transformation, in a completely automated manner. The last model considered was the cubic smoothing spline model based on a stochastic state space model that allows the use of a likelihood approach for estimating the smoothing parameter, thereby facilitating construction of the PI. It has been shown that this model is a special case of an ARIMA(0, 2, 2) model [9]. Its advantage over the full ARIMA model is that it provides a smooth historical trend as well as a linear forecast function. An overview of the models is provided as S1 Table in S1 File, along with specific equations.
The main outcomes of interest were the daily collected numbers of confirmed cases, deaths, and recoveries. We downloaded the data from the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University (https://github.com/CSSEGISandData/COVID-19 accessed on 04/28/2020), which consists of daily cumulative cases and covered our analysis period from February 22, 2020 until April 29, 2020. The unit of time used in modeling was day. Each outcome, beginning between March and April, exhibited a nearly exponential increase over time for both the USA and Italy (see S1 Fig in S1 File).
The forecasting performance of all these models were evaluated using mean absolute error (MAE) and mean absolute percentage error (MAPE). Additionally, the model fits were evaluated by using AIC (Akaike information criteria), reported in Table 1. We examined forecast durations in a weekly fashion, to assess model performance as extent of data collection increased over calendar time. All the aforementioned models are already implemented and available in the 'forecast' package in R, which makes implementation straightforward and possible in freely available software. The computer code that we used for our analyses is provided as S1 File.

Results
We performed the forecasting analysis for confirmed COVID-19 cases, number of deaths, and number of recoveries for the USA and Italy by using the models discussed above. For both countries, we conducted 7-day-or one-week ahead point forecasts and PIs. We updated our forecasts every 7 days under each of the four models considered. The results of our analyses for both countries are summarized with figures provided below. The accuracy metrics and fit statistics for each forecasting period, outcome, model and country are provided in Table 1.
Firstly, we considered applying these models to forecast the number of confirmed cases for the USA. The observed and forecasted confirmed COVID-19 cases in the USA with 90% PIs for each forecasting period and model are presented in Fig 1. The first forecast period covered April 2-April 8, 2020 where we used the data collected during the previous days. The forecasts and 90% PIs produced at the end of April 1 are presented in Fig 1. For the four models, the mean estimates and accompanying 90% PI for the confirmed cases 7-days ahead, expressed in rounded form in terms of thousands, were 389 (356, 422) for ARIMA, 484 (354, 644) for TBATS, 386 (347, 424) for Spline, and 389 (356, 423) for Holt. The actual number of confirmed cases was 428654 (rounded to 429 in terms of thousands) on April 8, which was captured by only the 90% prediction interval of TBATS model; this model had the lowest MAE and MAPE for this specific forecasting period (see Table 1). In this period predictions were less accurate than expected, possibly due to limited available data till April 1.
For the second forecast period, we incorporated another week of historical data by including the number of cases observed until April 8, 2020. This time the forecast period covered April 9-April 15, 2020. We again used the same models to conduct 7-day ahead predictions. The forecasts and 90% PIs produced at the end of April 7 are presented in Fig 1. The mean estimate (PI) for the confirmed cases 7-days ahead was 649 (611, 687), 701 (565, 857), 650 (616, 683), and 650 (613, 686) for the models ARIMA, TBATS, Spline, and Holt, respectively. The actual number of confirmed cases on April 14 was 636350 (rounded to 636 in terms of thousands), which was included by the 90% prediction intervals of all four models. Among all models, the Spline model had the lowest MAPE 8.1%, followed by the ARIMA model with MAPE 8.7% (see Table 1).
Furthermore, we conducted another set of forecasts and PIs for Period III which includes April 16-April 22, 2020. Similar to the analyses presented above, we used the data up until April 15, 2020. The mean estimate (PI) for 7-days ahead was 837 (797, 876), 854 (721, 1003), 834 (801, 868), and 837 (798, 876) for the models ARIMA, TBATS, Spline, and Holt, respectively. The actual number of confirmed cases on April 22 was 840351 (rounded to 840 in terms of thousands), which was well inside the prediction bands of all four models (see Fig 1). Among these models, the Spline model had the lowest MAPE which is 7.1%, but overall each method had quite small MAPE (see Table 1).
In the forecasting period IV, we considered forecasting for April 23-April 29, 2020 by using the data up until April 22, 2020. The mean estimate (PI) for 7-days ahead was 1040 (999, 1081), 1052 (929, 1183), 1041 (999, 1082), and 1040 (998, 1081) for the models ARIMA, TBATS, Spline, and Holt, respectively. The actual number of confirmed cases on April 29 was 1039909 (rounded to 1040 in terms of thousands) which is aligned with the predicted values by all models and included by their 90% PI. All the models had very low MAPE, slightly differing from each other, especially the Spline and ARIMA models which had MAPE< 7% (see Table 1).
In the final step, we used entire data history up till April 29, 2020, and we attempted forecasting for period V covering April 30 to May 06, 2020. Predicted confirmed cases and corresponding PIs for all models are presented in Fig 1. For this period, the point forecast (PI) for 7-days ahead was 1231 (1184, 1279), 1232 (1111, 1360), 1232 (1185, 1279), and 1231 (1183, 1279) for the models ARIMA, TBATS, Spline, and Holt, respectively. The data on April 29, 2020, were used to provide point forecasts and prediction intervals for the period V (April 30-May 06, 2020). At the time we performed our study, this period had not occurred; thus, we did not have the actual numbers of confirmed cases up to the concluding point of our analyses. The point forecasts and prediction intervals are provided in S2 Table in S1 File. As more data becomes available and is incorporated in the forecasting analysis, the model's predictive accuracies (e.g., MAPE) improved (see Table 1). For this period, the MAPE values are 6.1%, 8.4%, 6%, and 6.9% for ARIMA, TBATS, Spline, and Holt, respectively.
Similarly, we applied these models to forecast the confirmed cases for Italy for the same forecasting periods, and the results are presented in Fig 2. Based on our analyses, we concluded that the prediction performance of all these models were similar to what we observed in the USA data application. The values of actual observed confirmed cases were contained inside the prediction bands of all four models for all five forecasting periods. In particular, the models ARIMA, Spline and TBATS had quite small MAPE and good prediction accuracies (see  Table 1). The TBATS model had a wider PI for the first four forecasting periods but behaved similar to the PIs of the other three models for period V (see Fig 2). Additionally, the point forecasts and prediction intervals for the period V are provided in S3 Table in S1 File. We observed decreased in the forecasting accuracy metrics of all four models, compared to the previous period in both applications, which is expected as more data becomes available from forecasting periods I to V.
We next applied the same models to forecast the number of deaths due to COVID-19 in the USA for the aforementioned forecasting periods, although these models did not account for all factors that affect mortality of patients infected with COVID-19. The point forecast and PIs that we produced are presented in Fig 3. Again, for this application, the prediction performance of all these models was similar, and the values of actual observed number of deaths were within the PI bands of all four models for at least part of the forecasting periods. The ARIMA and Spline models had smaller prediction errors, compared to the other two models. The forecasting accuracy metrics and AIC values for this application are also provided in Table 1. Overall, the observed number of deaths were included by the PIs from these models except the predictions for forecasting periods I and III (see Fig 3). Additionally, the point forecasts and prediction intervals for the period V (April 30-May 06, 2020) are provided in S4 Table in S1 File. For this set of analyses, the models had higher MAPE than the analyses we produced for the confirmed cases for both the USA and Italy.
We also conducted forecasting for the number of deaths caused by COVID-19 in Italy by using the same models. Similar to previous analyses, the point forecast and prediction intervals for these analyses were produced and presented in Fig 4. In this application, we observed that the prediction performance of ARIMA, Splines, and TBATS models were better than the   performance of the Holt model, which yielded much higher prediction errors compared to others. However, the PIs of all four models included the observed number of deaths for all forecast periods. Similar to previous analyses, the forecast accuracy metrics and model fit statistics of this application are summarized in Table 1. Further, the point forecasts and PIs for period V (April 30-May 06, 2020) are provided in S5 Table in S1 File. For this set of analyses, the models yielded substantially lower MAPE than the analyses that was conducted for the number of deaths in the USA.
We now considered forecasting for the number of recoveries from COVID-19 in the USA by applying the same set of models. We obtained the point forecasts and prediction intervals for different forecast period that are presented in Fig 5. The MAE, MAPE, and AIC values of this analysis are provided in Table 1. The prediction performance of ARIMA, Splines, and Holt models were much better than the performance of the TBATS model, which yielded much higher prediction errors and quite wide PIs, especially for the first three forecasting periods in this application (for these periods the prediction intervals of TBATS model are trimmed to have better visualization in Fig 5). However, the performance of the TBATS model improved as more data became available. The PIs from the other three models did not contain the observed number of recoveries for forecast period IV: likely due to an unexpected jump/ change in the trajectory of the number of recoveries in the USA. The forecasting accuracy metrics for this application are summarized in Table 1. The point forecasts and PIs for period V (April 30-May 06, 2020) are presented in S6 Table in S1 File.
Similarly, we then applied the model to generate forecasts for the number of recoveries from COVID-19 in Italy. Point forecasts and PIs for the aforementioned forecast periods are presented in Fig 6. For this application, we found that the prediction performance of ARIMA and Holt models were slightly better than the performance of the Spline and TBATS models. The TBATS model again had poor performance compared to the other models in the analyses for the first two forecasting periods, but its performance improved for the later periods based on our visual inspection (it still has quite high MAPE) with an increased amount of data. The PIs from all four models include the observed numbers of recoveries for all forecast periods. Also, the forecasting accuracy metrics for this application are summarized in Table 1. As we do for previous applications, the point forecasts and prediction intervals for the period V (April 30-May 06, 2020) are provided in S7 Table in S1 File.

Discussion
We have evaluated four different highly cited time series models for forecasting confirmed cases, number of deaths and number of recoveries of COVID-19 for both the USA and Italy. Our application demonstrated the similarities in performance across countries for key COVID-19 outcomes but also elucidated where predictive performance may be lacking in certain models based on data accrual and sudden, sharp changes in trajectory. Although past publications have largely focused on the forecasting the number of the confirmed cases, we extended our study by considering number of deaths and recoveries. We focused on the data for two countries; however, we illustrated how our approach can be applied to data arising from any country or worldwide data. Similarly, the models considered in our work can be implemented on new data as it becomes available. The forecasting analyses should be updated as the data accumulate periodically to increase accuracy of forecasting, since this is an ongoing pandemic and additional data would improve forecasting performance of these models.
Using multiple models enabled us to test and compare their prediction accuracy and make an optimal selection. Our study showed how selecting a COVID-19 prediction model among existing, powerful options for time series models can readily assist in decision-making through projections regarding important issues such as intervention-induced changes in the spread of disease or making plans for controlling the pandemic and returning to normal life [5]. Such models could be useful in instances where knowledge is limited regarding the progress of the pandemic in a country or the world; assumptions about the future of this emergent disease (COVID-19 pandemic) can be made based on these established models and their available implementations on an open-source platform. Although we considered a breadth of models and outcomes, our work did not account for individualized factors that affect death and recovery. Incorporating such factors is expected to improve accuracy and narrow the width of PIs. Subject-level data applications are needed to address this issue. We should note that these approaches can produce forecasts for longer time periods, such as several weeks or months ahead. We examined forecast durations in a weekly fashion in this work given the nature of the forecasting goal to assess immediate windows of time.
Our study showed that all models resulted in similar overarching conclusions, provided there exists a reasonable amount of historical data available; however, each model exhibited subtle differences. The ARIMA model was the most consistent model across our different applications and could be preferable over the other models if one must pick a single model for application. ARIMA and cubic smoothing spline models had smaller prediction error and narrower PIs, compared to the Holt and TBATS models in the majority of the analyses. This is a surprising finding for the spline models, as this approach tends to perform worse in extrapolation settings. Indeed, these models were not robust in all settings. The TBATS model had quite large widths on the PIs especially when a limited amount of data was available, which is not preferable in applications with constant data accrual. Wider intervals correspond to a higher degree of uncertainty for the point forecasts. On the other hand, the Holt model always had a substantially larger AIC, which indicates poorer model fit, compared to the others for both the USA and Italy applications. TBATS and Spline models mostly resulted in similar AIC values and were often slightly smaller than the AIC of the ARIMA model (the smaller AIC, the better model fit) in the analysis for the USA data. The ARIMA model had the smallest AIC in all forecast periods for the Italy application. However, the model that has better fit (smaller AIC) does not always have a better prediction performance. One should judiciously select a model that has smaller prediction error while better reflecting the nature of the data. The forecasting accuracy metrics of all four models mostly improved with the increased amount of available data as days pass during this pandemic in both the USA and Italy. The unstable patterns in historical data on the number of deaths are more likely to worsen the accuracy of forecasts, but for the future forecasts these models are capable of accounting for major changes in the historical trajectory. In conclusion, these readily available models that can be performed efficiently at any time for any country or worldwide have reasonable accuracy for short-term forecasting of this pandemic, which should not be ignored.