II Theory

L I S T O F T A B L E S

Table Deterministic exact approaches to mean-variance portfolio selection problem (see also [START_REF] Mansini | Twenty Years of Linear Programming Based Portfolio Optimization[END_REF] and references within). References are sorted in chronological order; papers published in the same year are sorted according to alphabetic order of the last name of the corresponding first author.

A B S T R A C T

This thesis presents a new algorithm for Mixed Integer NonLinear Programming, inspired by the Multiplicative Weights Update framework and relying on a new class of reformulations, called the pointwise reformulations. The thesis is divided in three main parts: a foreword consisting in Chapter 1, a theoretical foundation of the new algorithm in Chapter 2, and the application of this new methodology to two real-world optimization problems, namely the Mean-Variance Portfolio Selection in Chapter 3, and the Multiple NonLinear Knapsack Problem in Chapter 4.

Mixed Integer NonLinear Programming is a hard and fascinating topic in Mathematical Optimization both from a theoretical and a computational viewpoint. These problems are characterized by nonlinear objective function and constraints, and continuous and integer decision variables. Many real-world problems can be cast this general scheme and, usually, are quite challenging in terms of efficiency and solution accuracy with respect to the solving procedures. Another very important tool in Mathematical Optimization is represented by formulations and reformulations: in particular, we introduce a new family of reformulations, namely pointwise reformulations, depending on a given parameter, which are easier to solve than the original formulation. A remarkable characteristic we look for in the pointwise reformulation is exactness, i.e., the existence of a given value of the parameter such that a global optimum of the original problem is also a global optimum for the reformulation. The basic idea to heuristically solve Mixed Integer NonLinear Problems consists in finding the optimum of the (easier) exact pointwise reformulation, which immediately yields the corresponding global optimum of the original problem. We employ the Multiplicative Weights Update algorithm in order to identify the correct value of the parameter for the pointwise reformulation.

In Chapter 1 we give an overview of the mathematical concepts and entities we use in the rest of the thesis. Chapter 2 is devoted to illustrate the general scheme of the new algorithm, the Multiplicative Weights Update for Mixed Integer NonLinear Programming, and its main theoretical properties. Moreover, in this chapter we define several automatic building procedures to determine the pointwise reformulation of a given Mixed Integer NonLinear Problem for specific, but broad, classes of optimization problems.

In the rest of the thesis we deal with two real-world challenging optimization problems: Mean-Variance Portfolio Selection and Multiple NonLinear Knapsack Problems. In Chapter 3 we give a survey on the models, formulations and reformulations, and exact methods for the single-objective singleix period Mean-Variance Portfolio Selection problem. Among all the versions for the portfolio problems proposed in the specialized literature, we choose the most addressed class, namely cardinality constrained portfolio selection with semi-continuous variables. We consider also a possibly non-convex nonconcave transaction cost function, with the only hypothesis of separability, which is quite natural in the context of financial markets. The Multiplicative Weights Update for Mixed Integer NonLinear Problems behaves sufficiently better than the benchmarks with respect to the quality of the solution produced and the number of assets which compose the optimal portfolios. In general, in fact, minimizing the number of assets in the optimal portfolio is a second goal for the optimal selection procedures.

Then, in Chapter 4 we consider the Multiple NonLinear Knapsack Problem, addressed here for the first time in its entirely. We adapt the Multiplicative Weights Updated framework to this problem, proposing a new pointwise reformulation. Unfortunately, extensive computational experiments show this algorithmic approach is not well suitable to solve challenging instances of this knapsack problem. Hence, we illustrate a different heuristic method based on the discretization of the solution space and on the surrogate relaxation. The method consists of three phases: we propose a constructive greedy procedure, and two procedures for the feasibility recovering of the surrogate solution. A local search post-procedure is also implemented in order to improve the overall quality of the solution produced by the heuristics. Computational experiments indicate that this method prevails over the benchmarks both in terms of quality of the solution and of total computational elapsed time.

R É S U M É

L'objectif de cette thèse consiste à présenter un nouvel algorithme pour la programmation non linéaire en nombres entiers, inspirée par la méthode Multiplicative Weights Update et qui compte sur une nouvelle classe de reformulations, appelées les reformulations ponctuelles. La thèse est divisée en trois parties principales: une introduction composée par le Chapitre 1, une définition théorique du nouvel algorithme dans le Chapitre 2 et l'application de cette nouvelle méthodologie à deux problèmes concrets d'optimisation, tels que la sélection optimale du portefeuille avec le critère moyennevariance dans le Chapitre 3 et le problème du sac à dos multiple non linéaire dans le Chapitre 4.

La programmation non linéaire en nombres entiers est un sujet très difficile et fascinant dans le domaine de l'optimisation mathématique à la fois d'un point de vue théorique et computationnel. Ces problèmes sont caractérisés par une fonction objective et des contraintes non linéaires, ainsi xi que des variables de décision continues et entières. Il est possible de formuler de nombreux problèmes dans ce schéma général et, habituellement, ils posent de réels défis en termes d'efficacité et de précision de la solution obtenue quant aux procédures de résolution. Un autre outil très important dans l'optimisation mathématique est représenté par les formulations et reformulations. En particulier, nous introduisons une nouvelle famille de reformulations, appelées reformulations ponctuelles, en fonction d'un paramètre donné. Elles sont plus simples à résoudre que la formulation originale. Une caractéristique remarquable recherchée dans la reformulation ponctuelle est l'exactitude, c'est-à-dire l'existence d'une valeur donnée du paramètre telle que un optimum global du problème d'origine est aussi un optimum globale pour la reformulation. L'idée de base pour résoudre heuristiquement les problèmes non linéaires en nombres entiers consiste à trouver l'optimum des reformulations ponctuelles exactes (les plus faciles), qui produit immédiatement l'optimum global correspondant du problème d'origine. Nous employons alors l'algorithme Multiplicative Weights Update afin d'identifier la valeur correcte du paramètre pour la reformulation ponctuelle.

Dans le Chapitre 1, nous définissons les concepts et les objets mathématiques utilisés dans le corps de la thèse. Le Chapitre 2 est consacré à illustrer le cadre général du nouvel algorithme, le Multiplicative Weights Update pour la programmation non linéaire en nombres entiers et ses principales propriétés théoriques. En outre, dans ce chapitre, nous définissons plusieurs procédures de construction automatique pour déterminer la reformulation ponctuelle d'un problème non linéaire en nombres entiers pour des classes spécifiques, mais larges, de problèmes d'optimisation.

Dans le corps de la thèse, nous nous occupons de deux problèmes d'optimisation difficiles: la sélection du portefeuille moyenne-variance et le problème du sac à dos multiple non linéaire. Dans le Chapitre 3, nous donnons un résumé des modèles, formulations et reformulations, ainsi que les méthodes spécifiques orientées sur le problème de la sélection du portefeuille moyenne-variance avec un seul objectif et sur une période unique. Parmi toutes les versions du problème de portefeuille proposé dans la littérature spécialisée, nous avons choisi la catégorie la plus abordée, appelée la sélection de portefeuille avec contrainte de cardinalité avec variables semicontinues. Nous considérons également une fonction de coût de transaction non-concave et non-convexe, avec la seule hypothèse de séparabilité, ce qui est tout à fait naturel dans le contexte des marchés financiers. Le Multiplicative Weights Update pour les problèmes non linéaires en nombres entiers fait preuve d'un meilleur comportement par rapport aux autres méthodes de résolution, notamment en termes de qualité de la solution produite et du nombre d'actions qui composent le portefeuille optimal. En général, en fait, la minimisation du nombre d'actions dans le portefeuille optimal est un deuxième objectif pour les procédures optimales de sélection.

Par conséquent, dans le Chapitre 4, nous considérons le problème du sac à dos multiple non linéaire, abordé ici pour la première fois dans son intégralité. Nous adaptons la structure du Multiplicative Weights Update à ce problème, proposant une nouvelle reformulation ponctuelle. Malheureusement, des expériences computationnelles poussées montrent que cette approche algorithmique n'est pas bien adaptée pour résoudre les cas difficiles de ce problème de sac à dos. Ensuite, nous illustrons une méthode heuristique différente basée sur la discrétisation de l'espace de solutions et sur sa relaxation agrégée. La méthode consiste en trois phases: nous proposons une procédure gloutonne constructive et deux algorithmes pour la récupération de la faisabilité de la relaxation agrégée. Une post-procédure de recherche locale est également exécutée afin d'améliorer la qualité globale de la solution produite par les heuristiques. Les expériences de calcul indiquent que cette méthode prévaut sur les autres tant en termes de qualité de la solution que de temps total de calcul.

I N T R O D U C T I O N

This thesis is devoted to introduce and analyze a new methodology to solve optimization problems, i.e., a broad class of problems, in which we want to find the minimum of a single deterministic objective function, in a finite or infinite set of feasible points, described by a finite number of inequalities and possibly implicit constraints such as integrality or membership in given polyhedra. In particular, we consider Mixed Integer NonLinear Problems (MINLPs), i.e., optimization problems involving continuous and discrete decision variables and possibly nonlinear terms in the objective function and the constraints.

The methodology we employ in this thesis involves three main ingredients: MultiStart algorithm, Multiplicative Weights Update algorithm, and reformulations. The MultiStart algorithm is a simple random (heuristic) procedure to globally solve optimization problems. The Multiplicative Weights Update algorithm can be explained as a stochastic (heuristic) method for a decision maker to iteratively take a decision among different choices, by observing the prediction of a finite number of advisors. Reformulations are a fundamental tool in optimization both from a theoretical and applicative viewpoint. The formulation describes the structure of a Mathematical Programming (MP). Reformulations change the symbolic structure of a MP formulation while keeping some of its mathematical properties invariant.

This chapter constitutes a foreword to the rest of the thesis: we introduce the notation and the mathematical entities we will use in the other chapters. The rest of this chapter is organized as follows. In Section 1.1 we give the mathematical formal definition of Mixed Integer NonLinear Programming (MINLP) and we remark its general properties. In Section 1.2 we describe a general classification for the algorithms for MINLP. In Sections 1.3 and 1.4 we describe the MultiStart and the Multiplicative Weights Update algorithms, respectively. In Section 1.5 we present a definition of formulations and reformulations of a given optimization problem, and we propose a general classification, based on the relationship between the original problem and its reformulated versions. Finally, in Section 1.6 the thesis structure is drawn with all the complete references from which the others chapters are sourced.

.

Let x ∈ R n be an n-dimensional vector of continuous decision variables and y ∈ Z p be a p-dimensional vector of integer variables. The general MINLP is defined as follows:

min f(x, y) (1.1a) s.t. g(x, y) 0 (1.1b) x ∈ X (1.1c) y ∈ Y ∩ Z p , (1.1d)
where f(x, y) : R n+p → R and g(x, y) : R n+p → R m represent the objective function and the constraints, respectively. The sets X ⊆ R n and Y ⊆ R p are two polyhedra of suitable dimensions. Let F be the feasible set of the MINLP, i.e.,

F := {(x, y) : g(x, y) 0, x ∈ X, y ∈ Y ∩ Z p }. (1.2)
Moreover, we define projections of the feasible set over the continuous and discrete variables, respectively:

F X :={x ∈ X : ∃y ∈ Y ∩ Z p such that g(x, y) 0} (1.3a) F Y :={y ∈ Y ∩ Z p : ∃x ∈ X such that g(x, y) 0}. (1.3b) Definition 1.1.1.
With a slight abuse of notation, we say that a MINLP is convex if f(x, y) and g (x, y) are convex for all ∈ {1, . . . , m}, otherwise the MINLP is non-convex.

Remark 1.1.2. We emphasize that the continuous relaxation of a convex MINLP, i.e., the nonlinear problem obtained by removing the integrality requirements of the variable y, has a convex feasible set.

Definition 1.1.3. With a slight abuse of notation, we say that a MINLP is strictly convex if f(x, y) and g (x, y) are strictly convex for all ∈ {1, . . . , m}.

If p = 0, the MINLP reduces to NonLinear Problem (NLP). If the objective function and the constraints are linear, the MINLP reduces to Mixed Integer Linear Problem (MILP). Finally, if p = 0 and the objective function and the constraints are linear, the MINLP reduces to Linear Problem (LP).

Then we introduce the definitions of neighborhoods of a given point

(x, ỹ) ∈ R n × Z p . Definition 1.1.4. Let (x, ỹ) ∈ R n × Z p , then B((x, ỹ); ε) := {(x, y) ∈ R n × Z p : y = ỹ, ||x -x|| q ε} (1.4a) N(x, ỹ) := {(x, y) ∈ R n × Z p : x = x, ||y -ỹ|| 0 1}, (1.4b)
where q ∈ N and ε > 0.

Henceforth, solving a MINLP means finding at least one global solution for the optimization problem, i.e., a feasible point whose objective value is the minimum among all the points in the feasible set.

Definition 1.1.5. A point (x * , y * ) ∈ F is a global solution of a MINLP if f(x * , y * ) f(x, y) for all (x, y) ∈ F.
We introduce also the definition of local optimum of a MINLP.

Definition 1.1.6. A point (x * , y * ) ∈ F is a local solution of a MINLP if, for some ε > 0, f(x * , y * ) f(x, y) for all x ∈ B((x * , y * ); ε) ∩ F X and y ∈ N(x * , y * ) ∩ F Y .
A global solution is given by the best solution in all the feasible set: this definition takes into account the global behavior of the objective function and the entire feasible set. A local solution, on the contrary, is a solution restricted to a small subset of the feasible set, with respect to which the local behavior of the objective function is considered.

While a strictly convex NLP has at most one global solution, the same property does not necessarily hold for strictly convex MINLPs [START_REF] Leyffer | Deterministic Methods for Mixed Integer Nonlinear Programming[END_REF][START_REF] Loh | A Sequential Linearization Approach for Solving Mixed-Discrete Nonlinear Design Optimization Problems[END_REF].

MINLP is NP-hard because it includes Mixed Integer Linear Programming (MILP) [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF][START_REF] Kannan | On the Computational Complexity of Integer Programming Problems[END_REF] and Mixed Integer Quadratic Programming (MIQP) [START_REF] Del Pia | Mixed-Integer Quadratic Programming is in NP[END_REF] as special cases, when constraints are affine and objective function is linear or quadratic, respectively. MINLP is, in general, undecidable [START_REF] Jeroslow | There Cannot be any Algorithm for Integer Programming with Quadratic Constraints[END_REF], even for p 10, when the objective function is linear and the constraints are polynomial [START_REF] De Loera | Integer Polynomial Optimization in Fixed Dimension[END_REF]. For a survey on computational complexity of MINLP we refer the interested reader to Hemmecke et al. [START_REF] Hemmecke | Nonlinear Integer Programming[END_REF] and Köppe [START_REF] Köppe | On the Complexity of Nonlinear Mixed-Integer Optimization[END_REF]. However, if the MINLP is convex or polyhedra X and Y are bounded, the undecidability is fortunately avoided.

Usually, solving MINLPs in practice could be extremely difficult, as the following example shows.

Example 1.1.7. Fermat's Last Theorem, formulated by Pierre de Fermat in 1637 on the margin of his copy of Arithmetica by Diophantus of Alexandria and published by his son Samuel de Fermat in Tolosa in 1670 (see [64, p. 62]) states: "Cubum autem in duos cubos, aut quadratoquadratum in duos quadratoquadratos et generaliter nullam in infinitum ultra quadratum potestatem in duos eiusdem nominis fas est dimidere cuius rei demonstrationem mirabilem sane detexi. Hanc marginis exiguitas non caperet.". 1This conjecture, finally proved by Wiles [START_REF] Wiles | Modular Elliptic Curves and Fermat's Last Theorem[END_REF], asserts the Diophantine equation x n + y n = z n has no integer solution when n ∈ Z ∩ [3, +∞) and x, y, z ∈ Z + ∩ [1, +∞). Now, Fermat's conjecture is false if and only if the optimum value of the following MINLP is zero [START_REF] Murty | Some NP-complete Problems in Quadratic and Nonlinear Programming[END_REF]:

min x ,y ,z ,n (x n + y n -z n ) 2 (1.5) s.t. x ∈ Z + ∩ [1, +∞) (1.6) y ∈ Z + ∩ [1, +∞) (1.7) z ∈ Z + ∩ [1, +∞) (1.8) n ∈ Z ∩ [3, +∞).
(1.9)

For simplicity, from now on, we assume all the optimization problems we introduce are characterized by lower and upper bounds on the decision variables, i.e., by constraints such as x ∈ [x, x] and y ∈ [y, y], where the underline symbol indicates the lower bounds and the upperline symbol indicates the upper bounds.

. The algorithms proposed in the literature to solve optimization problems can be subdivided into two different classes, with respect to the type of solution produced:

• local algorithms produce local solutions;

• global algorithms produce global solutions.

Optimization algorithms can also be subjected to a further classification:

• exact algorithms;

• heuristic algorithms.

Exact algorithms produce a solution which meets a given optimization criterion which characterized a local or a global solution [START_REF] Bonnans | Numerical Optimization: Theoretical and Practical Aspects[END_REF][START_REF] Fletcher | Practical Methods of Optimization[END_REF][START_REF] Nocedal | Numerical Optimization[END_REF], while heuristic algorithms produce, in general quickly, a good feasible solution for the problem. Global algorithms consider the entire feasible set and they try to explore it in order to find a global optimum, exactly or heuristically. According to the way they search the feasible set, they can be:

• deterministic methods;

• stochastic methods.

Deterministic methods are such that, if applied several times to the same instance of the optimization problem, they produce the same output every time; instead, stochastic methods are characterized by probabilistic procedures, which may produce different solutions each time the algorithm is executed on the same input.

One of the first general heuristic methods for global optimization problems, i.e., mathematical programs in which we want to find the global minimum, is the Pure Random Search (PRS) algorithm. The PRS (see Algorithm 1) is based on a simple restart procedure: at each iteration a new (feasible) point is generated (Step 3) and it is valued with respect to the objective function, if the current point has a better objective function value, the current returned point is updated (Steps 4-8).

Algorithm 1 Pure Random Search 1: Initially set t=1, (x * 0 , y * 0 ) := (x 0 , y 0 ) and U := f(x 0 , y 0 ) 2: while termination condition is not met do t := t+1 ) and upper bound U.

Step 3 can be implemented in many different ways from pure systematic generation to pure randomization procedure [START_REF] Martí | Multi-Start Methods[END_REF]. The termination condition is, in general, given by the maximum number of iterations or, in the case the global minimum value is known, the distance between the current function value and the global one [START_REF] Schoen | Stocastic Global Optimization: Stopping Rules[END_REF].

Remark 1.2.1. Although Algorithm 1 is quite simple, the implementation of this method cannot be a trivial task, since generating points in a feasible set is generally a non easy operation. One possible implementation is to generate a point in a box containing the feasible set and accept the point only if it is feasible. Assumption 1.2.2. We assume the feasible set is defined only by lower and upper bounds constraints: in this case Step 3 can be implemented in a pure randomization fashion, i.e., the current point is generated according to a uniform distribution over the feasible set. Proposition 1.2.3. Let Assumption 1.2.2 hold, and {(x t , y t )} be a sequence of uniformly random distributed points on the feasible set F of the optimization problem. Let Z t = {(x 1 , y 1 ), . . . , (x t , y t )}. Then, for all subset A ⊂ F with strictly positive Lebesgue measure,

lim t→∞ P (Z t ∩ A = ∅) = 1.
(1.10)

Proof. Let meas(A) denote the Lebesgue measure of set A. p = meas(A)/ meas(F) ∈ (0, 1) is the probability that a point sampled in F belongs to A. Therefore, we have:

P (Z t ∩ A = ∅) = 1 -(1 -p) t . (1.11)
The statement follows.

Proposition 1.2.3 guarantees the entire feasible set is covered by the sequence of random points: no part of the feasible set remains unexplored. The following proposition states that the Algorithm 1 converges to the global optimum in probability as the number of iterations goes to infinity.

lim t→∞ P ((x * t , y * t ) ∈ {(x, y) ∈ F : f(x, y) f * + ε}) = 1, (1.12)
where f * is the optimal value of the optimization problem.

Proof.

Let A = {(x, y) ∈ F : f(x, y) f * + ε}. The statement follows by invoking Proposition 1.2.3.
The serious drawback of Algorithm 1 consists in the fact that, in order to reach a point in the neighborhood of the global solution of the optimization problem, many iterations may be necessary: in the worst case the number of iterations is infinite. The probability of finding a point (x t , y t ) such that f(x t , y t ) f * + ε is ε/meas(F) and the probability of finding such point in N iterations is .13) In order to produce a point in the neighborhood of the optimal solution of the problem with a confidence level α, i.e., such that

1 -1 - ε meas(F) N . ( 1
1 -1 - ε meas(F) N = α, we need N = log (1 -α) log(1 -ε/meas(F)) (1.14)
iterations. Therefore, as pointed out by Locatelli and Schoen [START_REF] Locatelli | Global Optimization: Theory, Algorithms, and Applications[END_REF], if the optimization problem is pure continuous, i.e., p = 0, the feasible set F is a unit box, and the neighborhood of the global minimum is a box of edge length , the number of iterations needed is .15) .

log (1 -α) log(1 -n ) = O 1 n . ( 1
In order to speed up Algorithm 1 a local search phase can be also implemented, defining the MultiStart (MS) algorithm: at Step 4 in Algorithm 2 a local algorithm is applied to the randomly generated point. Obviously, in practical implementation, this local phase is crucial for the computational behavior of the algorithm since it is, in general, much more expensive than the global one [START_REF] Liberti | Writing Global Optimization Software[END_REF].

Remark 1.3.1. Note that, when the mathematical program is a non-convex NLP, Step 4 is NP-hard [START_REF] Floudas | Quadratic Optimization[END_REF][START_REF] Pardalos | Checking Local Optimality in Constrained Quadratic Programming is NP-hard[END_REF], however it is practically fast. Moreover, we can consider local NLP optimization as tractable as long as we limit ourselves to constraints involving only factorable functions over a simple operator alphabet such as +, -, * , /, log, exp.

Algorithm 2 MultiStart Algorithm

1: Initially set t := 1, (x * 0 , y * 0 ) := (x 0 , y 0 ) and U := f(x 0 , y 0 ) 2: while termination condition is not met do apply a local algorithm from (x t , y t ), obtaining point (x t , y t )

5: if f(x t , y t ) < U then 6: set x * t := x t , y * t := y t , U := f(x * t , y * t ) 7: else 8: set x * t := x * t-1 , y * t := y * t-1 9:
end if In order to avoid too many local minimization procedures, variants of Algorithm 2 have been proposed: in particular, the local phase could be performed only if the randomization phase produces a better point in terms of the objective function [START_REF] Liberti | Writing Global Optimization Software[END_REF].

Moreover, a population of points could be randomly generated in the randomization phase and a certain subset of points could be selected by means of clustering techniques (see [START_REF] Boender | Stochastic Methods[END_REF][START_REF] Liberti | Writing Global Optimization Software[END_REF][START_REF] Schoen | Two-Phase Methods for Global Optimization[END_REF]) and local minimizations can be performed, considering the selected points as starting points [START_REF] Liberti | Writing Global Optimization Software[END_REF]. In other versions of the MS, an escaping strategy is also implemented in order to avoid local minima: for instance, Simulated Annealing [START_REF] Kirkpatrick | Optimization by Simulated Annealing[END_REF] (see also [START_REF] Metropolis | Equation of State Calculations by Fast Computing Machines[END_REF][START_REF] Cerny | Thermodynamical Approach to the Travelling Salesman Problem: An Efficient Simulation Algorithm[END_REF][START_REF] Kirkpatrick | Optimization by Simulated Annealing: Quantitative Studies[END_REF][START_REF] Dekker | Global Optimization and Simulated Annealing[END_REF][START_REF] Locatelli | Simulated Annealing Algorithms for Global Optimization[END_REF]) is considered in the context of non-convex MINLP by Cardoso et al. [START_REF] Cardoso | A Simulated Annealing Approach to the Solution of MINLP Problems[END_REF] and Tabu Search [START_REF] Glover | Tabu Search, Part 1[END_REF][START_REF] Glover | Tabu Search, Part 2[END_REF][START_REF] Glover | Tabu Search[END_REF] is applied to non-convex MINLP by Munawar et al. [START_REF] Munawar | Solving Extremely Difficult MINLP Problems Using Adaptive Resolution Micro-GA with Tabu Search[END_REF] (see also [START_REF] Lin | Using Tabu Search to Solve MINLP Problems for PSE[END_REF][START_REF] Munawar | arGA: Adaptive Resolution Micro-Genetic Algorithm with Tabu Search to Solve MINLP Problems Using GPU[END_REF]).

Finally, the generation of the random points can be implemented by means of non uniform distributions [START_REF] Solis | Minimization by Random Search Techniques[END_REF], such as in Simulated Annealing: in this case the following proposition guarantees the convergence of the algorithm in probability. [START_REF] Solis | Minimization by Random Search Techniques[END_REF]) Let meas(A) be the Lebesgue measure of set A. Let {m t (•)} be a sequence of probability measures such that, for all (Borel) subset A ⊂ F with meas(A) > 0,

Proposition 1.3.2. (Solis and Wets

∞ t=1 (1 -m t (A)) = 0. (1.16)
Then, if f is a Lebesgue-measurable function and F ⊆ R n is a Lebesgue-measurable set, the sequence of random point (x * t , y * t ) generated by the algorithm is such that

lim t→∞ P ((x * t , y * t ) ∈ {(x, y) ∈ F : f(x, y) f * + ε}) = 1, (1.17)
where ε > 0.

. The Multiplicative Weights Update (MWU) algorithm is a "meta-algorithm", i.e., it could be adapted to many different settings, with a broad application in Optimization, Machine Learning, and Game Theory. In particular, we point out the Plotkin-Shmoys-Tardos Algorithm [START_REF] Plotkin | Fast Approximation Algorithm for Fractional Packing and Covering Problems[END_REF] for fractional packing and covering LPs, which can be derived by the MWU framework. The Plotkin-Schmoys-Tardos Algorithm is connected with the Lagrangian relaxation approach, in which several "complicated constraints" are relaxed. Let A ∈ R n×n , x, b ∈ R n and P be a convex subset of suitable dimension. The Plotkin-Schmoys-Tardos Algorithm solves the following feasibility problem:

∃?x ∈ P : Ax b (1.18)
in which x ∈ P represent the "easy constraints", while Ax b indicate the "complicated constraints". The Algorithm calls repetitively a sub-procedure (Oracle) which deals with the following feasibility problem:

∃?x ∈ P : p Ax p b, (1.19)
where p is a vector of suitable dimensions. Arora et al. [3] show that, if the sub-procedure is an ( , ρ)-bounded Oracle, then there exists an algorithm which solves the problem (1.18) up to an additive error ε or derives that the problem (1.18) is infeasible. The algorithm calls the Oracle only O( ρ log(m)/

2 ) times.

For a general explanation of MWU we refer the reader to the excellent survey by Arora et al. [3], from which we borrow the following example.

In a very simple stock market with just one stock characterized by two possible daily price movements (up and down), an investor has the possibility to observe the prediction of q experts. The investor wants to gain as much as possible in terms of overall return, according to the prediction of the best expert which, however, is of course not known a priori. The first (trivial) algorithm one could think about consists in operating in the market, selling or buying the stock, according to the majority opinion of the experts. Nevertheless, the experts could be correlated or even not really experts in finance, so the majority opinion could be systematically wrong.

The MWU algorithm (see Algorithm 3) corrects the trivial one: the investor chooses his/her financial strategy according to the weighted majority opinion of the experts. Initially all the experts have the same weights w i (i q) (Step 1), but, as the time goes on, the algorithm gives a gain to the experts who made the correct prediction, and gives a cost to ones who made a wrong prediction (Step 4). For technical reasons the costs/gains must be in [-1, 1]. If this is not the case, a step to suitably scale the costs/gains has to be implemented. The weight of each expert is updated in a multiplicative fashion (Step 5). At each iteration, the investor chooses on the basis of the predictions in a random way according to the probability distribution induced by the weights (Step 3).

Algorithm 3 Multiple Weights Update Algorithm

1: Initially set t:=1, η 1 2 and w i,t := 1 for all i q 2: while termination condition is not met do 3:

sample i q from the distribution p t ∼ (w i,t : i q).

4:

each decision incurs a cost/gain ψ t ∈ [-1, 1] q .

5:

update weights w i,t+1 := w i,t (1 -ηψ i,t ).

6:

t := t + 1.

7: end while

Even though the MWU is a quite simple strategy, it is possible to show an upper bound for the overall expected cost t T p t ψ t . Theorem 1.4.1 (Arora et al. [3]). Let T be the number of iterations, ψ i,t ∈ [-1, 1], for all i ∈ q and t T , be the cost/gain associated to expert i at time t and η 1 2 . The MWU Algorithm guarantees an overall expected cost

E MWU := t T ψ t p t min i q   t T ψ i,t + η t T |ψ i,t |   + ln q η (1.20)
.

Intuitively a formulation is a way to write down a given optimization problem. Specific types of formulation are, for instance, required by the implemented solvers in order to address a specific type of problems, such as LPs, NLPs, MILPs, and MINLPs. Formulations can be:

• flat formulations;

• structured formulations.

In the objective and constraints of structured formulations, quantifiers, such as ∀, , , appear; in flat formulations no quantifier is present. When a formulation P is cast into another formulation Q, we say that Q is a reformulation of P.

Several definitions for reformulation have been proposed (see, e.g., Audet et al. [START_REF] Audet | Links betwwn Linear Bilevel and Mixed 0-1 Programming[END_REF] and Sherali [START_REF] Sherali | Personal Communication[END_REF]). Generally, reformulations are defined in such a way that several properties of the original formulation are preserved, such as the set of the optimal solutions or the set of the feasible points. A systematic theory for reformulations is presented in Liberti [START_REF] Liberti | Reformulations in Mathematical Programming: Definitions and Systematics[END_REF], which proposes the following classification. Definition 1.5.1. A reformulation Q of a formulation P is a relaxation if its feasible set contains the feasible set of P. Definition 1.5.2. A reformulation Q of a formulation P is exact if it shares all the optimization properties (local optima, global optima, feasible set) with P.

Moreover, in the next chapters we will consider also bounding reformulations [START_REF] Mencarelli | A Multiplicative Weights Update Algorithm for MINLP[END_REF].

Definition 1.5.3. A reformulation Q of a formulation P is bounding if, when solved to optimality, produces a lower bound for P and its feasible set contains the feasible set of P. Remark 1.5.4. Note that all the relaxations are also bounding reformulations.

. This thesis is based on several published papers, namely [START_REF] Mencarelli | A Multiplicative Weights Update Algorithm for MINLP[END_REF][START_REF] Mencarelli | Heuristics for the General Multiple Non-linear Knapsack Problem[END_REF], and other working papers [START_REF] Mencarelli | Complex Portfolio Selection via Convex Mixed-Integer Quadratic Programming: A Survey[END_REF][START_REF] Mencarelli | Relaxations and Heuristics for the General Multiple NonLinear Knapsack[END_REF] submitted to the international refereed journals Computers & Operations Research and International Transactions in Operational Research, respectively. In particular, Chapter 2 is sourced from [START_REF] Mencarelli | A Multiplicative Weights Update Algorithm for MINLP[END_REF], Chapter 3 from [START_REF] Mencarelli | Complex Portfolio Selection via Convex Mixed-Integer Quadratic Programming: A Survey[END_REF], and Chapter 4 from [START_REF] Mencarelli | Relaxations and Heuristics for the General Multiple NonLinear Knapsack[END_REF][START_REF] Mencarelli | Heuristics for the General Multiple Non-linear Knapsack Problem[END_REF].

Beside the introductory section, the thesis is structured in two main parts: a theoretical part (Chapter 2), and an applicative one (Chapters 3 and 4). In the theoretical part we introduce a new algorithm to solve MINLP, and we explain its fundamental steps and its theoretical properties. In the second part we apply the methodology to two different optimization problems, namely the Mean-Variance Portfolio Selection and the Multiple NonLinear Knapsack Problems, both modeled as MINLPs and difficult to solve in practice. We will see that for the first problem the algorithm performs quite well with respect to the benchmarks, while this is not the case for the second problem.

Hence, for the sake of completeness, we describe a constructive heuristic procedure to find good solutions in a reasonable amount of computational time.

Part II.

Theory

T H E M W U A L G O R I T H M F O R M I N L P .

The MWU algorithm can be described as a special case of the MS algorithm to solve MINLPs. The MS algorithm is composed of two main steps: the choice of a random starting point, and a local optimization procedure, meaning we solve the problem with a local solver (for instance, we can heuristically solve non-convex problems by means of a solver which solves convex problems exactly). In the MWU algorithm for MINLPs we follow the same structure of the MS, but we introduce a strategy for generating promising points for the objective function and for the constraints. It is quite clear what "promising point" means in terms of objective function: if we consider a minimization problem and two feasible points, the one with the lowest objective function value is more promising than the other. Intuitively, being a "promising point" for the constraints corresponds to a sufficiently little violation of the constraints in terms of the difference between the value of the left hand side of the constraint and the value of the right hand side. Moreover, since in the MWU algorithm for MINLPs we generate "promising points" according to the problem formulation, in order to guarantee low objective function values and low violation of the constraints, we say that MWU algorithm for MINLPs is a matheuristic [START_REF] Maniezzo | Hybridizing Metaheuristics and Mathematical Programming[END_REF], i.e., a heuristic algorithm based on the MP formulation. In other words, we generate the promising point by solving an auxiliary formulation of the original problem, where several terms are fixed to the values given by the MWU, we call a pointwise reformulation. The solution we obtain is the starting point for a local optimization procedure. Then, we iterate the two steps (solving the pointwise reformulation and applying a local procedure) for a given number of iterations or until a given optimality criterion is satisfied. If the pointwise reformulation is built in such a way that there exists a value of the parameter for which the reformulation has the same optimum as the original one, i.e., the pointwise reformulation is exact, then we only need to guess the correct value of the parameter and solve the simpler reformulation. If the pointwise reformulation is bounding, at each iteration we obtain a lower bound for the optimal value of the original formulation, by simply solving the reformulated problem.

The rest of this chapter is organized as follows. In Section 2.2 we formally introduce the pointwise reformulation, describing its theoretical properties and characterizations. In Section 2. Given a MINLP P as in (1.1), a term of P is a symbolic expression in its expression tree, i.e., it is represented by a set of adjacent nodes in its expression tree. The expression tree of a MINLP is a way to computationally represent the problem in a tree, having as leaves the variables and the constants of the problem and as the other nodes the logical operators (+, ×, /, sin, etc.) linking variables and constants [START_REF] Cohen | Computer Algebra and Symbolic Computation: Mathematical Methods[END_REF].

Both the original terms t and the substituting ones t are functions of the decision variables x and y; furthermore, the substituting terms t are also functions of the parameters θ: hence, we extensively indicate the substituted terms as t(x, y) and the substituting terms as t (x, y; θ). Definition 2.2.2. Given a MINLP P and R θ = ptw t←t'(θ) (P), both defined on a vector

x of continuous decision variables in R n and a vector y of discrete decision variables in Z p : (a) For every replaced term t s (for s r) in Definition 2.2.1, let D s be the range of t s (x, y), where the term is interpreted as a function of the decision variables x and y of P ranging in the respective domains.

(b) For every replacement term t' s (for s r) in Definition 2.2.1, let D s (θ s ) be the range of t' s (x, y; θ s ) when the term is interpreted as a function of the decision variables x and y of P ranging in the respective domains.

(c) For s r, let Θ s be the range of the corresponding parameter θ s , and let Θ = (Θ s | s r).

Finally, given a parameter θ ∈ Θ and a function φ, we indicate with φ θ the function obtained by replacing the terms t s (x, y) with t s (x, y; θ) (s r). Let X θ and Y θ be the sets obtained by replacing the terms t s (x, y) into the inequalities defining sets X and Y, respectively. With the previous notation, we rewrite the pointwise reformulation as follows:

min f θ (x, y) (2.1a) s.t. g θ (x, y) 0 (2.1b) x ∈ X θ (2.1c) y ∈ Y θ ∩ Z p . ( 2
D (θ s ) ∧ t' s (x, y; θ s ) = t s (x, y) (2.2)
(b) R θ is exact if, for each globally optimal solution (x * , y * ) of P, there is at least one vector θ ∈ Θ such that (x * , y * ) is also an optimal solution of R θ ;

(c) R θ is efficient if there is a polynomial-time algorithm for approximately solving R θ (for θ ∈ Θ) to within a given ε > 0 approximation factor.

Lemma 2.2.5. Let P be a MNLP cast in form (1.1) and R θ = ptw t←t (P) be a spanning reformulation, then we have:

feas(P) ⊆ θ∈Θ feas(R θ ). (2.3)
where feas(P) is the feasible set of problem P.

Proof. Let (x , y ) ∈ feas(P). Given that R θ is spanning, there exists a parameter ξ ∈ Θ such that t s (x , y ) = t s (x , y ; ξ s ) for each s r. Hence it follows g ξ (x * , y * ) = g (x * , y * ) for all m. Note (x , y ) ∈ feas(P), g (x , y ) 0 for all m, therefore (x , y ) is also feasible for R ξ . Since for all ξ ∈ Θ feas(R ξ ) is a subset of θ∈Θ feas(R θ ), the statement follows.

Lemma 2.2.6. Let P be a given formulation and R θ be a spanning pointwise reformulation, then there exists ξ ∈ Θ such that R ξ is a bounding reformulation of P.

Proof. In the proof of Lemma 2.2.5 we saw that, if ξ ∈ Θ is such that t s (x * , y * ) = t s (x * , y * ; ξ s ) for all s r, then (x * , y * ) ∈ feas(R ξ ). Analogously, it is possible to prove f ξ (x * , y * ) = f(x * , y * ), and therefore:

val(R ξ ) f ξ (x * , y * ) = f(x * , y * ) = val(P), (2.4)
where val(P) indicates the value of the optimal solution of P. The statement of the theorem follows.

.

. . A First-Order Reformulation

In this section we derive a pointwise linear reformulation for a general non-convex MINLP. We assume functions f and g ( m) are at least once continuously differentiable. Moreover, we suppose without loss of generality that the objective functions of the original problem P and of its pointwise reformulation R θ are the same, i.e., f θ ≡ f. If this is not the case, we add a dummy variable γ ∈ R to P, and we consider the following exact reformulation, in the sense of [START_REF] Liberti | Reformulations in Mathematical Programming: Definitions and Systematics[END_REF], of the original MINLP problem: min γ (2.5a)

s.t. f(x, y) γ (2.5b) g(x, y) 0 (2.5c) x ∈ X (2.5d) y ∈ Y ∩ Z p . ( 2

.5e)

We assume all the inequalities describing sets X and Y are included on inequalities g(x, y) 0 with g(x, y) : R n+p → R m . We replace each nonconvex multivariate function in (2.5) with an affine approximation, i.e., with a first-order approximation at a given point (x, ỹ): min γ (2.6a)

ν 00 + ν 10 x - x y -ỹ γ (2.6b) ν 0 + ν 1 x - x y -ỹ 0 ∀ m (2.6c) y ∈ Z p . (2.6d) Remark 2.3.1.
In the pointwise reformulation (2.6), the parameter θ is the matrix ((x ỹ) , ν 0k , ν 1k ) (k ∈ {0, . . . , m}), whose dimensions are (n + p) × (2m + 3).

For simplicity of notation, in the rest of this subsection, we set g 0 (x, y) := f(x, y). Lemma 2.3.2. The pointwise reformulation (2.6) is spanning.

Proof. For each k ∈ {0, . . . , m}, we note that the replacement term ν 0k + ν 1k (xx yỹ) and the replaced term g k (x, y) have the same value for all values of θ such that (x, ỹ) = (x, y) and ν 0k = g k (x, y).

The previous approach to generate a pointwise reformulation has been successfully applied to a special class of the Hydro Unit Commitment (HUC) problems arising in energy industry in the paper [START_REF] Mencarelli | A Multiplicative Weights Update Algorithm for MINLP[END_REF].

. . Polynomial MINLPs

In this section, we consider polynomial MINLPs, i.e., formulations (1.1) where functions f and g ( m) are polynomials of (x, y). For simplicity in this section we assume p = 0. All the results, except where explicitly indicated, are valid also in case discrete decision variables are present.

A Trivial Pointwise MILP Reformulation

A trivial method to derive a pointwise MILP reformulation for any polynomial MINLP P consists in turning every variable but one to a parameter in every monomial of P. For two integer sequences h = (h

a µ = θ d -1 h j = θ d j h j .
Since every monomial is reduced to a linear term, but the rest of the constraints remains unchanged, it is immediate to show that the formulation obtained by the above procedure is a pointwise MILP reformulation of P.

We point out that, if, with respect to the monomial µ d h , there exists an index k t, such that d k = 1, because of the rule (a), we do not introduce θ h k in the pointwise reformulation keeping the number of the parameters as low as possible.

Remark 2.3.3. By construction, the previous pointwise reformulation is spanning since, for each feasible x d , the replaced term µ is always equal to the replacement term x h θ d -1 h j = θ d j h j for θ h j = x h j (j t).

. . Bilinear MINLPs

The trivial pointwise MILP reformulation appears particularly well suited to the case of bilinear MINLPs where the incident graph corresponding to the matrix of the quadratic form is bipartite, i.e., polynomial MINLPs of degree 2 where the variables can be partitioned into two sets I and J such that every monomial x i x j has i ∈ I and j ∈ J. Then the procedure reduces to choose I (resp. J) and replace x i with θ i (resp. x j with θ j ) for all i ∈ I (resp. j ∈ J).

. . Quadratic MINLPs

Consider MINLPs involving quadratic terms only and assume without loss of generality p = 0:

min x Q 0 x + q 0 x (2.7a) s.t. x Q x + q x b ∀ m (2.7b)
Let L be the set of indices in {0, . . . , m} such that Q is indefinite. For each ∈ L, we look for two matrices A , B such that:

(i) A is positive semidefinite;

(ii) B is sparse;

(iii) A + B = Q .
Next, we rewrite x Q x + q x as x A x + x B x + q x, introduce as many parameters θ as there are variables x, and define the pointwise reformulation by replacing every term x B x by θ B x. This yields a pointwise convex Mixed Integer Quadratic Problem (MIQP) reformulation, since all the nonconvex terms x B x have been turned into linear terms. Thereby, given the indefinite matrix Q , we aim the following optimization problem to recover A 0 and B :

min B 0 (2.8a) s.t. A 0 (2.8b) A + B = Q . ( 2
sample θ t from the distribution p t ∼ (w i,t | i q) 4: solve R θ t = ptw t←t'(θ t ) (P), get solution (x t , y t ) 5:
optionally refine (x t , y t ) (e.g. using local descent)

6: if (x t , y t ) is better than the incumbent, set (x * , y * ) = (x t , y t ) 7: (x t , y t ) yields decision costs ψ t ∈ [-1, 1] q 8:
update weights for next iteration:

w i,t+1 ← w i,t (1 - ψ i,t 2 ) 9:
increase t.

10: end while 11: return (x * , y * )

. . Sampling

In MWU, a weight is maintained for each advisors. In this case we want to estimate parameters θ: hence, we have to take q decisions at each step, not only one. In Step 3 we randomly select r values for the parameters according to the distribution of the weights: we have one weight, and therefore one expert, for each decision.

In general, the number r of the parameters and the dimension q of the weights w could be different. Techniques of aggregation, if r < q, or disaggregation, if, on the contrary, r > q, must be implemented before the sampling step. For instance, in [START_REF] Mencarelli | A Multiplicative Weights Update Algorithm for MINLP[END_REF] an example of disaggregation methodology is given for the Euclidian Distance Geometry Problem (DGP).

. . Solution and Refinement

Usually, solving the pointwise reformulation is not an easy task if the reformulation is not efficient. However, if the reformulation is good and easier to solve than the original formulation, the proposed methodology can be successfully exploited. This is the case, for instance, when the original prob-lem P is a non-convex MINLP and the reformulation is a convex MINLP or even a MILP. Remark 2.4.1. In order to obtain an upper bound, we do not need to solve the pointwise reformulation exactly, but we can use any heuristic we like. Several termination criteria can be considered, such as the number of iterations or the optimality gap.

The refinement procedure (Step 5) is not expressly needed to guarantee the upper bound on the cumulative error of MWU algorithm. Nevertheless, computational experiments show it is necessary in practice to speed up the algorithm. Moreover, if we consider Algorithm 4 as a special case of the two-phase MS algorithm, the refinement step coincides with a local phase in which we locally improve the initial solution.

. . Computing MWU Costs/Gains

This is the most critical step of the MWU algorithm, since it can influence the performance guarantee. It has basically two requirements: (a) each cost vector ψ t should have components in [-1, 1]; (b) if θ t replaces several terms in the original formulation, ψ t reflects the contribution of its value to optimality and feasibility. Requirement (a) is necessary for the MWU performance guarantee to hold, while requirement (b) tries to correlate E MWU (see (1.20)) to suboptimality and infeasibility.

Hence, ψ i,t consists of a scaled contribution α t to the error of (x i , y i ) from the suboptimality of (x t , y t ) in P and of β t from its infeasibility. Since it is not necessary to reward "better feasibility", we consider β t belonging to [0,1].

Let f θ (x, y) and g θ (x, y) 0 ( m) be respectively the objective function and the constraints of R θ .

After Steps 4-5, we can evaluate the current solution (x t , y t ) in the original formulation P by computing f(x t , y t ) and g (x t , y t ) for each . We let α t be proportional to f θ (x t , y t ) -f(x t , y t ) and β t = (β ,t : m) to (max(g θ (x t ), 0) : m). These quantities have, however, to be scaled in order to guarantee ψ t = [-1, 1] q . For simplicity, in the rest of the Section, we assume the original problem is continuous, i.e., p = 0. Formally, here is how ψ t is computed:

1. for t = 1, let α t = sgn(f θ (x t ) -f(x t )) and β ,t = sgn(max(g (x t ), 0))
for all ;

2. for each t > 1 let:

α t = f θ (x t ) -f(x t ) max s t |f θ (x s ) -f(x s )| and β ,t = max(g (x t ), 0) max s t (max(g (x s ), 0)) for all ; 3. Ψ t = 1 m+1 (α t + m β ,t ); 4. for each s r let ψ s,t = Ψ t |t s (x t , y t ) -t s (x t , y t ; θ t )| max(|t s (x t , y t )|, |t s (x t , y t ; θ t |)
.

Remark 2.4.2. Ψ t could be the result of many different ways to combine of α t and the β ,t depending from the specific application.

Remark 2.4.3. Since Ψ t ∈ [-1, 1] and |t i (x t , y t ) -t i (x t , y t ; θ t )| max(|t i (x t , y t )|, |t i (x t , y t ; θ t |) ∈ [0, 1],
we have

ψ t ∈ [-1, 1].
.

In this chapter we have given a theoretical insight with regard to the MWU algorithm for MINLPs. In particular we have introduced the pointwise reformulation of a given MP formulation and we have stated its main properties. Moreover, we have drawn several building procedures for the pointwise reformulations with respect to several classes of MINLPs. Finally, we have discussed in detail the main steps of the MWU algorithm.

Part III.

Applications

3 M E A N -V A R I A N C E P O R T F O L I O S E L E C T I O N P R O B L E M .
The first milestone in modern single-period portfolio selection theory is undoubtedly Harry Markowitz's 1952 seminal paper [START_REF] Markowitz | Portfolio Selection[END_REF] (for an historic perspective, see also [START_REF] Rubinstein | Markowitz's Portfolio Selection: A Fifty-Year Retrospective[END_REF] and surveys [START_REF] Wallingford | A Survey and Comparison of Porfolio Selection Models[END_REF][START_REF] Constantinides | Portfolio Theory[END_REF][START_REF] Elton | Modern Porfolio Thory, 1950 to Date[END_REF]), in which the meanvariance portfolio optimization model was proposed for the first time. Although several ideas and results are already introduced by de Finetti [START_REF] De Finetti | Il Problema dei Pieni (in Italian)[END_REF] (see [START_REF] Barone | Bruno de Finetti. The Problem of Full-Risk Insurances[END_REF] for the English translation of the first chapter "The problem in a single accounting period"), the contribution of the Italian mathematician was discovered only recently by the financial international community (see [START_REF] Rubinstein | Bruno de Finetti and Mean-Variance Portfolio Selection[END_REF][START_REF] Rubinstein | A History of the Theory of Investments: My Annotated Bibliography[END_REF][START_REF] Pressacco | The Origins of the Mean-Variance Approach in Fnance: Revisiting de Finetti 65 Years Later[END_REF]) and was acknowledged by Harry Markowitz himself [START_REF] Markowitz | de Finetti Scoops Markowitz[END_REF].

The mean-variance approach is based on the fundamental observation that, according to what Markowitz states in [START_REF] Markowitz | Portfolio Selection[END_REF], the investors should try to increase their portfolio return and contemporaneously to decrease, as much as possible, its volatility or its risk (see also [START_REF] Markowitz | Portfolio Selection: Efficient Diversification of Investimen[END_REF][START_REF] Markowitz | Mean-Variance Analysis in Portfolio Choice and Capital Markets[END_REF]). The portfolio variance is the most widely used measurement of the portfolio volatility; other possible risk measurements are reported, for instance, in [START_REF] Chang | Portfolio Optimization Problems in Different Risk Measures Using Genetic Algorithm[END_REF][START_REF] Hoe | An Empirical Comparison of Different Risk Measures in Portfolio Optimization[END_REF][START_REF] Mansini | Twenty Years of Linear Programming Based Portfolio Optimization[END_REF]. If the expected returns of the assets follow a Gaussian distribution [START_REF] Hanoch | The Efficiency Analysis of Choices Involving Risk[END_REF] or the investor's utility function is quadratic [START_REF] Castellani | Manuale di Finanza[END_REF], then the mean-variance criterion is theoretically compatible with the expected utility hypothesis originally introduced by Bernoulli [START_REF] Bernoulli | Specimen Theoriae Novae de Mensura Sortis[END_REF] (see [START_REF] Sommer | Daniel Bernoulli. Exposition of a New Theory on the Measurement of Risk[END_REF] for the English translation and see also [START_REF] Fishburn | Utility Theory for Decision Making[END_REF][START_REF] Huang | Foundations for Financial Economics[END_REF][START_REF] Neumann | Theory of Games and Economic Behaviour[END_REF] for more modern approaches). As pointed out by Markowitz in [START_REF] Markowitz | Mean-Variance Approximations to Expected Utility[END_REF][START_REF] Markowitz | Portfolio Theory: As I Still See It[END_REF], the previous assumptions are sufficient, but not necessary conditions. However, the assumption of Gaussian asset returns might be unrealistic: the probability distribution for the expected returns is generally leptokurtic [START_REF] Mills | Stylized Facts on the Temporal and Distributional Properties of Daily FT-SE Returns[END_REF].

The resulting mathematical program might not represent completely the problem solved nowadays by the practitioners, but it can be enriched with various constraints to model the different characteristics of the modern financial markets. Moreover, the mean-variance approach considers only the firstand second-order moments of the probability distribution of the returns: consequently, in specific situations, this approach might lead to counterintuitive or even paradoxical solutions [START_REF] Copeland | Financial Theory and Corporate Policy[END_REF].

Kallberg and Ziemba [START_REF] Kallberg | Comparision of Alternative Utility Functions in Portfolio Selection Problems[END_REF] compare the effects of different utility functions with respect to the optimal portfolios when the distribution of the expected returns is Gaussian, and show empirically that utility functions with similar absolute risk aversion indices -defined by Arrow [START_REF] Arrow | Essays in the Theory of Risk-Bearing[END_REF] and Pratt [START_REF] Pratt | Risk Adversion in the Small and in the Large[END_REF], but 29 originally introduced by de Finetti [START_REF] De Finetti | Sulla Preferibilità (in Italian)[END_REF] (see also [START_REF] Castellani | Manuale di Finanza[END_REF][START_REF] Montesano | de Finetti and the Arrow-Pratt Measure of Risk Aversion[END_REF]) -give rise to similar optimal portfolios. The reminder of the chapter is organized into two parts. In the first part we survey convex MIQP approaches to solve the Mean-Variance Portfolio Selection (MVPS) problem: in the next section we present the mathematical formulation of the portfolio problem with quadratic risk measure according to Markowitz [START_REF] Markowitz | Portfolio Selection[END_REF] and we analyze its main disadvantages; Sections 3.3, 3.4, and 3.5 provide several ways to enrich the original formulation in order to overcome its drawbacks; several equivalent mathematical reformulations for the mean-variance probabilistic portfolio problem are described in Section 3.6; finally, in Section 3.7 we summarize exact methods proposed in the literature to solve MVPS problems.

In the second part of the chapter, in Section 3.8, we choose a specific family of portfolio problems with cardinality constraint and transaction costs and we apply the MWU algorithm to this class of uncountable many MVPS problems. The aim of this part consists in a computational study of the performances of MWU, with a particular focus on the behavior of MWU depending on the degree of nonlinearity of the cost function. We propose a promising pointwise reformulation for this class of problems and a procedure to compute costs/gains in the MWU framework. Computational experiments show that the algorithm outperforms the benchmarks with respect to the quality of the solution produced. A summary of the main notation used throughout the chapter is reported in the Appendix A. Sections 3.3.2 and 3.4.1-3.4.3 are entirely based on the papers by Bonami and Lejeune [START_REF] Bonami | An Exact Solution Approach for Portfolio Optimization Problems under Stochastic and Integer Constraints[END_REF] and Lejeune [START_REF] Lejeune | A VaR Black-Litterman Model for the Construction of Absolute Return Fund-of-Funds[END_REF]. Section 3.6.1 is completely based on [START_REF] Bonami | An Exact Solution Approach for Portfolio Optimization Problems under Stochastic and Integer Constraints[END_REF], while Sections 3.6.2-3.6.3 on [START_REF] Filomena | Stochastic Portfolio Optimization with Proportional Transaction Costs: Convex Reformulations and Computational Experiments[END_REF].

. We consider r possibly risky assets characterized by a mean return vector µ ∈ R r . Let x ∈ R r + be the vector, whose generic entry j (j = 1, . . . , r) represents the fraction of the portfolio value invested in asset j. For the moment, following Markowitz [START_REF] Markowitz | Portfolio Selection[END_REF], we assume that the entries of µ and of the covariance return matrix Σ ∈ R r×r are known precisely. In the meanvariance approach, we aim to minimize the portfolio variance x Σx under the constraint that the portfolio return is at least equal to a given level R > 0. Therefore, the problem we aim to solve can be stated as follows:

min x Σx (3.1a) s.t. µ x R (3.1b) 1 r x = 1 (3.1c) x 0, (3.1d)
where 1 r ∈ R r is the all-one vector. Then, by repeatedly solving problem (3.1) for different values of return R, we can compute the efficient frontier, i.e., the set of the non-dominated portfolios in the sense of Pareto optimality.

In several cases (see, e.g., [START_REF] Bonami | An Exact Solution Approach for Portfolio Optimization Problems under Stochastic and Integer Constraints[END_REF]) an additional non-risky asset with mean µ 0 and zero variance is also considered, in order to algorithmically derive the efficient frontier (see Section 3.4.8). Several papers (see, e.g., [START_REF] Crama | Simulated Annealing for Complex Portfolio Selection Problems[END_REF]) consider an equality version for the return constraint (3.1b), namely µ x = R.

Constraint (3.1c) ensures that the whole capital available is invested in the portfolio and in several papers (see, e.g., [START_REF] Burdakov | Mathematical Programs with Cardinality Constraints: Reformulation by Complementarity-type Conditions and a Regularization Method[END_REF]) is substituted by

1 r x 1. (3.2)
Constraint (3.1d) prevents short selling, i.e., the possibility for the investor to sell financial assets not already in his/her portfolio. This financial operation is generally performed with speculative intents when the investor expects a bearish trend in the financial stock market. In case short selling is allowed, (3.1c) can be replaced by the constraint 1 r x = 0, which defines the so-called dollar neutral portfolio, by requiring the exposure on long part of the portfolio to be equal to the one on the short part. Several authors consider the case where the decision variables x represent the absolute amount invested per asset so that the inequality (3.1c) becomes 1 r x B, where B is the investor's total initial budget.

In [START_REF] Buchheim | A Frank-Wolfe Based Branch-and-Bound Algorithm for Mixed-Integer Portfolio Optimization[END_REF], Buchheim et al. introduce the budget constraint:

v x B, (3.3)
where the decision variables x are the units of financial asset held in the investor's portfolio and v ∈ R r + is the vector of the costs per unit of corresponding asset. Problem (3.1) is a convex continuous linearly constrained quadratic program, because, by definition, matrix Σ is symmetric and positive semidefinite; hence, we have a computationally tractable problem. However, the main drawback of this model consists in the sensitivity of the optimal solutions with respect to the input parameters (expected returns and covariance matrix), which are clearly unknown in real-world applications (see, e.g., [START_REF] Britten-Jones | The Sampling Error in Estimates of Mean-Variance Efficient Portfolio Weights[END_REF][START_REF] Broadie | Computing Efficient Frontiers using Estimated Parameters[END_REF][START_REF] Ceria | Incorporating Estimation Errors into Portfolio Selection: Portfolio Construction[END_REF][START_REF] Chen | Sensitivity to Estimation Errors in Mean-Variance Models[END_REF][START_REF] Fastrich | Robust Portfolio Optimization with a Hybrid Heuristic Algorithm[END_REF][START_REF] Jobson | Confidence Regions for the Mean-Variance Efficient Set: An Alternative Approach to Estimation Risk[END_REF][START_REF] Jorion | Portfolio Optimization in Practice[END_REF][START_REF] Jobson | Putting Markowitz Theory to Work[END_REF][START_REF] Michaud | The Markowitz Optimization Enigma: is "Optimized" Optimal?[END_REF]). Furthermore, when Σ is estimated starting from empirical measurements, it might happen that semidefiniteness is not directly satisfied and some ad-hoc procedures are required (see [START_REF] Cournéjols | Optimization Methods in Finance[END_REF][START_REF] Higham | Computing the Nearest Correlation Matrix: A Problem from Finance[END_REF][START_REF] Schoen | Two-Phase Methods for Global Optimization[END_REF]).

Chopra [START_REF] Chopra | Mean-Variance Revisited: Near-Optimal Portfolios and Sensitivity to Input Variations[END_REF] empirically analyzes the effects of slight differences in the estimate. Best and Grauer [START_REF] Best | On the Sensitivity of Mean-Variance-Efficient Portfolios to Changes in Asset Means: Some Analytical and Computational Results[END_REF] conduct a theoretically rigorous analysis with computational results about the sensitivity of mean-variance efficient portfolios with respect to possible changes in asset means.

Several papers study the instability and ill-conditioning of problem (3.1): for instance, Kallberg and Ziemba [START_REF] Kallberg | Mis-Specifications in Portfolio Selection Problems[END_REF] consider estimation errors in the investor's utility function and the mean vector and covariance matrix of the return distribution for normally distributed portfolio selection problems and observe that errors in mean vector give rise to significant problems. Chopra and Ziemba [START_REF] Chopra | The Effect of Errors in Means, Variances, and Covariances on Optimal Portfolio Choice[END_REF] show that the estimating errors with respect to the expected return means is generally one order of magnitude larger than the one corresponding to estimating errors in asset variances or covariances, assuming negative exponential utility function with joint normal distribution of returns.

.

. . Robust Approaches

The robust version of the mean-variance problem (3.1) has been considered in quite recent works (see the surveys [START_REF] Fabozzi | Robust Portfolios: Contributions from Operations Research and Finance[END_REF][START_REF] Gabriel | Recent Advances in Robust Optimization: An Overview[END_REF]). It consists in assuming that the expected returns are uncertain and their expected values and variances belong to a given set. By supposing that the unknown input parameters belong to a given uncertainty set (see [START_REF] Ben-Tal | Robust Convex Optimization[END_REF][START_REF] Ben-Tal | Robust Solutions of Uncertain Linear Programs[END_REF][START_REF] Ghaoui | Robust Solutions to Least-Squares Problems with Uncertain Data[END_REF][START_REF] Ghaoui | Robust Solutions to Uncertain Semidefinite Programs[END_REF]), it is possible to show some theoretical results: for instance, Goldfarb and Iyengar [START_REF] Goldfarb | Robust Portfolio Selection Problems[END_REF] and Tütüncü and Koenig [START_REF] Tütüncü | Robust Asset Allocation[END_REF] established the robust portfolio selection problem can be formulated as a Second-Order Cone Problem (SOCPs) (see [START_REF] Alizadeh | Second-Order Cone Programming[END_REF][START_REF] Ben-Tal | Lectures on Modern Convex Optimization: Analysis, Algorithms, and Engineering Applications[END_REF][START_REF] Boyd | Convex Optimization[END_REF][START_REF] Lobo | Applications of Second-Order Cone Programming[END_REF]) for ellipsoidal and box uncertain sets, respectively.

Under the assumption that the return mean belongs to a convex polytope, whose vertices are known, Costa and Paiva [START_REF] Costa | Robust Portfolio Selection Using Linear-Matrix Inequalities[END_REF] prove that program (3.1) can be formulated as a Linear Matrix Inequalities (LMI) problem (see [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF][START_REF] Nesterov | Interior-Point Polynomial Algorithms in Convex Programming[END_REF]). Moreover, El Ghaoui et al. [START_REF] Ghaoui | Worst-Case Value-at-Risk and Robust Portfolio Optimization: A Conic Programming Approach[END_REF] show that, when the mean and the covariance are unknown, but bounded, the worst-case mean-variance portfolio selection problem can be reformulated as a Semidefinite Program (SDP) (see [START_REF] Ben-Tal | Lectures on Modern Convex Optimization: Analysis, Algorithms, and Engineering Applications[END_REF][START_REF] Nesterov | Interior-Point Polynomial Algorithms in Convex Programming[END_REF][START_REF] Saigal | Handbook of Semidefinite Programming and Applications[END_REF][START_REF] Vandenberghe | Semidefinite Programming[END_REF]).

Finally, Ye et al. [START_REF] Ye | Robust Portfolio Optimization: A Conic Programming Approach[END_REF] introduce uncertain sets both for the mean vector and the second moment matrix of the returns, showing the connection between the fully robust portfolio selection problem with box uncertain set for the mean and ellipsoid uncertain set for the second moment of the returns and SOCP, SDP, and Semi-Infinite Programming (see [START_REF] Žaković | Semi-Infinite Programming and Applications to Minimax Problems[END_REF]).

. . Probabilistic Approach

Bonami and Lejeune [START_REF] Bonami | An Exact Solution Approach for Portfolio Optimization Problems under Stochastic and Integer Constraints[END_REF] take into account the uncertainty in the expected assets returns by dealing with a probabilistic problem and by introducing a probabilistic constraint, which imposes that the expected return of the op-timal portfolio should be not less than a given return level R with a high probability p > 0.

Let ξ be the random vector representing the expected returns of the r risky assets. We assume that the random vector ξ admits a probability density function and the density function of ξ x is strictly positive. Moreover, let µ ∈ R r with µ = E [ξ] and Σ = E (ξ -µ)(ξ -µ) be the mean and the covariance matrix for the r-variate distribution of ξ, respectively. Formulation

min x Σx (3.4a) s.t. P ξ x R p (3.4b) 1 r x = 1 (3.4c) x 0 (3.4d)
is usually referred to as the probabilistic Markowitz formulation and its deterministic equivalent defines a NLP (see [START_REF] Bonami | An Exact Solution Approach for Portfolio Optimization Problems under Stochastic and Integer Constraints[END_REF][START_REF] Filomena | Stochastic Portfolio Optimization with Proportional Transaction Costs: Convex Reformulations and Computational Experiments[END_REF][START_REF] Filomena | Warm-Start Heuristic for Stochastic Portfolio Optimization with Fixed and Proportional Transaction Costs[END_REF][START_REF] Lejeune | Portfolio Optimization with Combinatorial and Downside Return Constraints[END_REF]). Let ψ = (ξ x -µ x)/ √

x Σx be the standardized random variable representing the normalized portfolio return. Equation (3.4b) can be equivalently rewritten, as follows:

P ξ x R = P ψ R -µ x √ x Σx = 1 -F (x) R -µ x √ x Σx , (3.5)
where F (x) (•) is the cumulative distribution of the standardized portfolio return. We assume that F (x) (•) is a continue strictly increasing function. Moreover, we point out that the analytic form of the probability distribution F depends on the portfolio weights x. It follows that the probabilistic constraint (3.5) becomes

1 -F (x) R -µ x √ x Σx p ⇐⇒ 1 -p F (x) R -µ x √ x Σx ⇐⇒ µ x + F -1 (x) (1 -p) √ x Σx R, (3.6)
where

F -1 (x) (•) is the inverse of the cumulative distribution F (x) (•) and F -1 (x) (1 -p) is the (1 -p)-quantile of F (x) (•).
Therefore, the deterministic equivalent of optimization problem (3.4) corresponds to the following NLP [START_REF] Kataoka | A Stochastic Programming Model[END_REF]:

min x Σx (3.7a) s.t. µ x + F -1 (x) (1 -p) √ x Σx R (3.7b) 1 r x = 1 (3.7c) x 0. (3.7d)
In the following, we survey for which classes of probability distributions the problem can be reformulated as a SOCP. We thus recall several definitions in probability theory and convex optimization. Definition 3.3.1. (Serfling [225]) Let ξ ∈ R r be a random variable, whose probability density function is f : R r → R. If f(ξ -θ) = f(θ -ξ), then ξ has a distribution that is centrally symmetric about θ ∈ R r . Definition 3.3.2. (Boyd and Vandenberghe [START_REF] Boyd | Convex Optimization[END_REF]) Let x ∈ R n be the decision variables and

A i ∈ R (n i -1)×r , H ∈ R p×r and h, c i ∈ R r , β i ∈ R n i -1 , g ∈ R p , d i ∈ R (∀i ∈ {1, .
. . , n}) be the parameters of a given convex continuous optimization problem and • 2 indicates the Euclidean norm. If a convex continuous optimization problem can be (equivalently) rewritten as follows,

min h x s.t. A i x + β i 2 c i , x + d i , i = 1, . . . , n Hx = g, (SOCP)
then it is an SOCP. A constraint is a Second-Order Cone Constraint (SOCC) of dimension n i , if it can be equivalently rewritten as

A i x + β i 2 c i x + d i . (SOCC)
Remark 3.3.3. (Boyd and Vandenberghe [START_REF] Boyd | Convex Optimization[END_REF]) Observe that SOCP generalizes Quadratically Constrained Programming, i.e., the case when c i = 0 for all i = 1, . . . , n: a Quadratically Constrained Problem can be obtained by squaring the constraints. SOCP generalizes also Liner Programming (LP), i.e., the case when instead A i = 0 (n i -1)×r for all i = 1, . . . , n, where 0 (n i -1)×r is the zero matrix of suitable dimensions. Definition 3.3.4. (Lobo et al. [START_REF] Lobo | Applications of Second-Order Cone Programming[END_REF]) A convex set C ⊆ R r is SOC-representable if it is equivalent to the intersection of a finite number of SOCC, i.e., there exist parameters

A i ∈ R (n i -1)×(r+m) , β i ∈ R n i -1 , c i ∈ R r+m , and d i ∈ R such that x ∈ C ⇐⇒ ∃y ∈ R m : A i x y + β i 2 c i x y + d i i = 1, . . . , n.
Moreover, a given function

f : R r → R is SOC-representable if the set {(x, t) : f(x) t} is SOC-representable.
Ultimately, given an objective function f : R r → R and a feasible convex set C ⊆ R r , which are SOC-representable, then the corresponding convex optimization program, i.e., min f(x)

s.t. x ∈ C, (SOC)
can be dealt with as a SOCP by means of efficient interior point methods (see, for instance, [START_REF] Freund | Interior Point Methods: Current Status and Future Directions[END_REF][START_REF] Nemirovskii | Interior-Point Methods for Optimization[END_REF][START_REF] Nesterov | Interior-Point Polynomial Algorithms in Convex Programming[END_REF][START_REF] Potra | Interior-Point Methods[END_REF][START_REF] Wright | The Interior-Point Revolution in Optimization: History, Recent Developments, and Lasting Consequences[END_REF]), characterized by polynomial time computational complexity [START_REF] Gondzio | Interior Point Methods 25 Years Later[END_REF][START_REF] Nesterov | Interior-Point Polynomial Algorithms in Convex Programming[END_REF]. Bonami and Lejeune [START_REF] Bonami | An Exact Solution Approach for Portfolio Optimization Problems under Stochastic and Integer Constraints[END_REF] showed some convexity results for problem (3.7) that we briefly discuss in the following for the sake of completeness.

Theorem 3.3.5. (Bonami and Lejeune [START_REF] Bonami | An Exact Solution Approach for Portfolio Optimization Problems under Stochastic and Integer Constraints[END_REF]) Let p ∈ [0.5, 1). If the probability distribution of ξ x is centrally symmetric, then the deterministic constraint (3.7b), equivalent to (3.5), is a SOCC.

Therefore, optimization problem (3.4) is an SOCP because its objective function is convex quadratic and its feasible region is described by the intersection of a second-order cone and several linear constraints. Constraint

µ x R -F -1 (x) (1 -p) √
x Σx ensures that the expected portfolio return is greater than the given return plus a penalty term, which is function of the portfolio variance and is increasing with the confidence level p [START_REF] Filomena | Stochastic Portfolio Optimization with Proportional Transaction Costs: Convex Reformulations and Computational Experiments[END_REF].

We recall also the definition of the skewness of a multi-variate distribution of a real-valued random variable ξ with mean µ and standard deviation σ [START_REF] Avèrous | Skewness for Multivariate Distributions: Two Approaches[END_REF]:

skew(ξ) = E [ξ -µ] 3 σ 3 , ( 3.8)
The skewness is basically an asymmetry index of the distribution: perfectly symmetric distributions have zero skewness.

Definition 3.3.6. (Bonami and Lejeune [START_REF] Bonami | An Exact Solution Approach for Portfolio Optimization Problems under Stochastic and Integer Constraints[END_REF]) A probability distribution of an r-variate real-valued random vector ξ with mean µ and median m has positive skewness if The exact value of the (1 -p)-quantile, F (x) (1 -p), is known only for few probability distributions. If we assume, for example, that the distribution of the expected returns is Gaussian, which is a quite restrictive assumption (see, for example, [START_REF] Fama | Mandelbrot and the Stable Paretian Hypothesis[END_REF][START_REF] Fama | The Behaviour of Stock-Market Prices[END_REF][START_REF] Mandelbrot | The Variation of Certain Speculative Prices[END_REF][START_REF] Rachev | An Empirical Examination of Daily Stock Return Distributions for U.S. Stocks[END_REF]), but rather common in several theoretical frameworks (see [START_REF] Hanoch | The Efficiency Analysis of Choices Involving Risk[END_REF][START_REF] Kallberg | Comparision of Alternative Utility Functions in Portfolio Selection Problems[END_REF][START_REF] Montesano | de Finetti and the Arrow-Pratt Measure of Risk Aversion[END_REF]), then the numerical values of quantiles F -1 (x) (1 -p) of the normalized portfolio return ψ are computationally known.

P (0 ψ) P (m ψ) ⇐⇒ F -1 (x) (α) 0, α 0.5, where E [ψ] = E [ξ -µ] = 0 and F (x) (m) = P (m ψ) = 0.5.
. Beyond the ill-conditioning of problem (3.1), the other serious drawback of Markowitz's original proposal is represented by the mismatch with the problems faced by practitioners in real-world applications (see, e.g., [START_REF] Cournéjols | Optimization Methods in Finance[END_REF][START_REF] Schoen | Stocastic Global Optimization: Stopping Rules[END_REF]). Nevertheless, we can consider additional constraints to problems (3.1) or (3.4), which describe the most common restrictions observed in real-world financial markets (see [START_REF] Fabozzi | Robust Portfolio Optimization and Management[END_REF][START_REF] Kolm | 60 Years of Portfolio Optimization: Practical Challenges and Current Trends[END_REF][START_REF] Lejeune | Portfolio Optimization with Combinatorial and Downside Return Constraints[END_REF][START_REF] Lobo | Portfolio Optimization with Linear and Fixed Transaction Costs[END_REF][START_REF] Mansini | Twenty Years of Linear Programming Based Portfolio Optimization[END_REF]). However, this kind of constraints could make the efficient frontier discontinuous and more challenging to compute [START_REF] Jobst | Computational Aspects of Alternative Portfolio Selection Models in the Presence of Discrete Asset Choice Constraints[END_REF].

. . Buy-in Thresholds

Generally, investors avoid extremely small long positions in their portfolios, because, on one side, they have a limited impact on the return value of the portfolio and, on the other side, they could be quite expensive with respect to finance fees and monitoring costs [START_REF] Scherer | Introduction to Modern Portfolio Optimization[END_REF]. Long positions not belonging to a given range [x j , x j ] ⊂ [0, 1] of the total initial budget B can be prevented by the simple range constraint x j x j x j , j = 1, . . . , r.

(3.9)

Several authors (see, e.g., [START_REF] Cesarone | A New Method for Mean-Variance Portfolio Optimization with Cardinality Constraints[END_REF][START_REF] Chang | Heuristics for cardinality constrained portfolio optimization[END_REF][START_REF] Frangioni | Perspective Cuts for a Class of Convex 0-1 Mixed Integer Programs[END_REF][START_REF] Jobst | Computational Aspects of Alternative Portfolio Selection Models in the Presence of Discrete Asset Choice Constraints[END_REF]) require x to be a semi-continuous variable [START_REF] Sun | Recent Advances in Mathematical Programming with Semi-continuous Variables and Cardinality Constraint[END_REF], i.e., they require x j ∈ [x j , x j ] ∪ {0} for all j = 1, . . . , r: they introduce extra binary variables δ ∈ {0, 1} r such that, for all j = 1, . . . , r, δ j = 1 if the investor holds the asset j, i.e. if x j > 0, and add the following constraints avoiding too small or huge holding positions:

x j δ j x j x j δ j , j = 1, . . . , r. (3.10)
Note that constraints (3.10) directly imply 0 x j δ j for all j = 1, . . . , r.

. . Round Lot Purchasing

Usually, investors manage only given lots of shares and other financial agreements, because of the facility in monitoring and purchasing/selling operations. Furthermore, for small private investors, splitting a large lot involves a premium, that has to be paid to the broker. The higher the cost of splitting large batch in single shares, the greater the impact of this kind of cost with respect to the optimal portfolios. Round lot purchasing constraints prescribe that investors hold, for the risky asset j (j = 1, . . . , r), batches or lots of S j stocks.

Let us define γ ∈ Z r + as a vector of general integer variables. We require that the number of the shares of asset j (j = 1, . . . , r), namely η j ∈ Z + , is an integer multiple of the lot-size S j :

η j = γ j S j , j = 1, . . . , r. (3.11)
Let q j be the market value of asset j (j = 1, . . . , r) held in portfolio, then we have η j = x j B/q j and constraint (3.12) can be equivalently rewritten as follows,

x j = q j γ j S j B , j = 1, . . . , r. (3.12)
The reader is referred to [START_REF] Bonami | An Exact Solution Approach for Portfolio Optimization Problems under Stochastic and Integer Constraints[END_REF] for further discussion.

Mansini and Speranza [START_REF] Mansini | Heuristic Algorithms for the Portfolio Selection Problem with Minimum Transaction Lots[END_REF] have shown that finding a feasible solution of problem (3.1) with round lot constraints (3.11), upper bound on γ j , i.e., the number (j = 1, . . . , r) of minimum lots, and bound constraints with respect to the total portfolio expenditure is NP-complete.

. . Sector Diversification

Generally, either there exist law limitations about the risk exposure (this is the case, for instance, of pension funds) or investors try to hold a representative portion of their portfolio in a prescribed number of asset categories or industrial sectors. However, in general, optimal portfolios for problem (3.1) are not well-diversified [START_REF] Green | When Will Mean-Variance Efficient Portfolios be well Diversified?[END_REF]. Usually, given are lower bound on the fraction of portfolio value held in specific sets of shares. For classical empirical analysis about financial benefits of a well-diversified portfolio, we refer the reader to [START_REF] Copeland | Financial Theory and Corporate Policy[END_REF][START_REF] Fama | Foundations of Finance: Portfolio Decisions and Securities Prices[END_REF][START_REF] Solnick | The Advantages of Domestic and International Diversification[END_REF].

Let us assume that every asset can be allocated to a specific financial category and let C k (k = 1, . . . , n) be the index set of all risky assets connected with the category k. Moreover, we suppose that sets C k define a partition of {1, . . . , r}. We introduce a binary variable ζ k ∈ {0, 1} for each financial category, such that ζ k = 1 if and only if the investment in financial category k (k = 1, . . . , n) is above a prescribed minimum level s:

sζ k j∈C k x j s + (1 -s)ζ k . (3.13)
Moreover, we have to consider an additional constraint in order to satisfy the diversification prescription [START_REF] Bonami | An Exact Solution Approach for Portfolio Optimization Problems under Stochastic and Integer Constraints[END_REF], which requires to hold portions of assets in at least n > 0 categories:

n k=1 ζ k n. (3.14)
. .

Cardinality Constraints

Beyond diversification requirements, asset managers (for instance in index tracking funds) wish to replicate as accurately as possible a market index with a limited number of financial agreements, namely K > 0. This can be modeled through the following cardinality constraint:

x 0 = r j=1 sgn(|x j |) K. (3.15)
By introducing additional decision variables δ j , already presented for constraints (3.10), we can straightforwardly reformulate the previous constraint in the following equivalent form [START_REF] Lejeune | Portfolio Optimization with Combinatorial and Downside Return Constraints[END_REF]:

r j=1 δ j K. (3.16)
Bienstock [START_REF] Bienstock | Computational Study of a Family of Mixed-Integer Quadratic Programming Problems[END_REF] (see also [START_REF] Shaw | Lagrangian Relaxation Procedure for Cardinality-Constrained Portfolio Optimization[END_REF]) shows that problem (3.1) with cardinality constraint (3.16) is NP-hard, even when r = 3. Several authors (see, e.g., [START_REF] Chang | Heuristics for cardinality constrained portfolio optimization[END_REF][START_REF] Deng | Ant Colony Optimization for Markowitz Mean-Variance Portfolio Model[END_REF][START_REF] Fernández | Portfolio Selection Using Neural Networks[END_REF][START_REF] Soleimani | Markowitz-based Portfolio Selection with Minimum Transaction Lots, Cardinality Constraints and Regarding Sector Capitalization using Genetic Algorithm[END_REF][START_REF] Woodside-Oriakhi | Heuristic Algorithms for the Cardinality Constrained Efficient Frontier[END_REF]) consider an equality version for cardinality constraint (3.16) and propose mainly heuristic methods to solve the corresponding problem:

r j=1 δ j = K.
(3.17)

Moreover, finding the K assets that should be included in the optimal portfolio is, in general, an NP-hard problem [START_REF] Moral-Escudero | Selection of Optimal Investment Portfolios with Cardinality Constraints[END_REF].

Using the theoretical results in [START_REF] Scozzari | A Clique Algorithm for Standard Quadratic Programming[END_REF][START_REF] Tardella | Connections between Continuous and Combinatorial Optimization Problems through an Extension of the Fundamental Theorem of Linear Programming[END_REF] and extending [START_REF] Cesarone | Efficient Algorithms for Mean-Variance Portfolio Optimization with Hard Real-Word Constraints[END_REF][START_REF] Cesarone | Efficient Algorithms for Mean-Variance Portfolio Optimization with Hard Real-Word Constraints[END_REF][START_REF] Cesarone | Portfolio Selection Problems in Practice: A Comparison Between Linear and Quadratic Optimization Models[END_REF], Cesarone et al. [START_REF] Cesarone | A New Method for Mean-Variance Portfolio Optimization with Cardinality Constraints[END_REF] have shown that the problem (3.1) with cardinality constraints (3.16) has the same optimal solution of problem (3.1) with equality cardinality constraints (3.17) and reduce this kind of programs to Standard Quadratic Programming Problem (see [START_REF] Bomze | On Standard Quadratic Optimization Problems[END_REF][START_REF] Bomze | New and Old Bounds for Standard Quadratic Optimization: Dominance, Equivalence and Incomparability[END_REF]), avoiding to explicitly introduce binary variables and considering an exact tailored solving procedure, called Increasing Set Algorithm. The Standard Quadratic Programming Problem is an NP-hard problem when the Hessian matrix of the objective function is indefinite, i.e., if the Hessian matrix of the objective function is neither positive nor negative semidefinite [START_REF] Bomze | On Standard Quadratic Optimization Problems[END_REF].

Di Gaspero et al. [START_REF] Panigrahi | Hybrid Local Search for Constrained Financial Portfolio Selection Problem[END_REF] consider an "interval" version for the cardinality constraint (3.16): K r j=1 δ j K, where K and K are such that 1 K K r. Cardinality constraints are closely related to buy-in threshold constraints [START_REF] Jobst | Computational Aspects of Alternative Portfolio Selection Models in the Presence of Discrete Asset Choice Constraints[END_REF]. Finally, in several papers (see, e.g., [START_REF] Cesarone | A New Method for Mean-Variance Portfolio Optimization with Cardinality Constraints[END_REF][START_REF] Chang | Heuristics for cardinality constrained portfolio optimization[END_REF][START_REF] Jobst | Computational Aspects of Alternative Portfolio Selection Models in the Presence of Discrete Asset Choice Constraints[END_REF]) it is observed that the problem (3.1), with cardinality constraints (3.16) and with minimum and maximum buy-in thresholds (3.10) can be straightforwardly reformulated as a convex MIQP.

. . Sector Capitalization

Sector capitalization constraints are introduced by Soleimani et al. [START_REF] Soleimani | Markowitz-based Portfolio Selection with Minimum Transaction Lots, Cardinality Constraints and Regarding Sector Capitalization using Genetic Algorithm[END_REF], in order to mathematically formulate the behavior of investors generally inclined to hold assets in financial sectors with higher capitalization value to reduce the total portfolio risk.

Let be the number of economic sectors and suppose, without loss of generality, that they are sorted in non-increasing way according to their capitalization value. Define L l as the set of assets for economic sector l (l ∈ {1, . . . , }). We introduce additional binary variables y l such that

1 M j∈L l δ j y l M j∈L l δ j l ∈ {1, . . . , } (3.18a) j∈L l µ j + (1 -y l ) j∈L l+1 µ j l ∈ {1, . . . , -1}, (3.18b)
where M ∈ R + is a sufficiently large positive number. The "big-M" constraints (3.18) ensure that the assets belonging to the sectors with higher capitalization values have basically higher probability to be in the optimal portfolios than the ones belonging to sectors with less capitalization values.

. . Turnover and Trading

Frequently, investors already hold a portfolio x (0) and, because of mutations in the financial market or others, they want to change their portfolio, by considering the new financial environment and by limiting, however, the variations with respect to the portfolio already held [START_REF] Perold | Large-Scale Portfolio Optimization[END_REF].

Crama and Schyns [START_REF] Crama | Simulated Annealing for Complex Portfolio Selection Problems[END_REF] propose to introduce restrictions on purchasing and selling variations. In particular, let P j and S j be respectively the maximum purchasing and selling levels for asset j (j = 1, . . . , r), turnover constraints can be stated as follows:

max x j -x (0) j , 0 P j j = 1, . . . , r (3.19a) max x (0) j -x j , 0 S j j = 1, . . . , r. (3.19b)
Because of fixed transaction costs (see Section 3.4.1), additional constraints are, generally, introduced in order to prevent small variations between portfolios. Let P j and S j be respectively the minimum purchasing and selling levels for asset j, trading disjunctive constraints can be stated as follows,

x j = x (0) j ∨ x j x (0) j + P j ∨ x j x (0)
j -S j for all j = 1, . . . , r .

. . Benchmark Constraints

Often, investors want to obtain a portfolio which is as close as possible to a benchmark (or target) portfolio x B [START_REF] Bertsimas | Portfolio Construction Through Mixed-Integer Programming at Grantham, Mayo, Van Otterloo and Company[END_REF]. With respect to economic sector diversified investments, Bertsimas and Shioda [START_REF] Bertsimas | Algorithm for Cardinality-Constrained Quadratic Optimization[END_REF] introduce the following additional constraints in order to bound variances between the optimal and the target portfolios:

j∈S l (x j -x B j ) ε l l = 1, . . . , . (3.20)
. .

Collateral Constraints

Di Gaspero et al. [START_REF] Gaspero | Local Search for Constrained Financial Portfolio Selection Problems with Short Sellings[END_REF] (see also [START_REF] Jacobs | Portfolio Optimization with Factors, Scenarios, and Realistic Short Positions[END_REF]) discuss the following legal constraints for short selling portfolios imposed by US Regulation T, a set of US laws concerning the margin requirements for the collateral agreement. The complete text of the regulation is available at https://www.ecfr.gov/cgi-bin/ text-idx?tpl=/ecfrbrowse/Title12/12cfr220 _ main _ 02.tpl. In particular, they introduce a free-risk asset with mean µ 0 and zero variance, the so-called collateral agreement, such that

x 0 -a r j=1 min{0, x j } (3.21a) r j=0 |x j | 2 (3.21b)
where a ∈ N + is the security level for the collateral agreement. In this case the decision variables x are not constrained to be positive, since short-selling is allowed, and variables δ defined in (3.10) are replaced by ternary variables z ∈ {-1, 0, 1} r , such that, for each j (j = 1, . . . , r), z j = 1 if the investor bought the asset j, i.e., if x j > 0, z j = -1 if the investor sold the asset j, i.e., if x j < 0, and z j = 0 if the investor does not hold asset j. Therefore, cardinality constraint (3.16) becomes r j=1 |z j | K.

. Besides (3.1a), several different objective functions have been proposed in the literature in order to make problems (3.1) and (3.4) simpler with respect to computational tractability or to better model real behaviors of investors and money savers. We consider only objective functions involving quadratic risk measure, namely portfolio variance (for an exhaustive survey on approaches proposed for portfolio selection problem with linear risk measures we refer the interested reader to the paper [START_REF] Mansini | Twenty Years of Linear Programming Based Portfolio Optimization[END_REF] and to the recent book [START_REF] Mansini | Linear and Mixed Integer Programming for Portfolio Optimization[END_REF]).

. . Penalty Functions

In order to define an unconstrained NLP, Bartholomew-Biggs and Kane [START_REF] Bartholomew-Biggs | A Global Optimization Problem in Portfolio Selection[END_REF] introduce the following penalty function for problem (3.1) with minimum buy-in threshold constraints (3.9) with x j := x and x j := 1 for all j ∈ {1, . . . , r},

φ(x j ) = 4 x j (x j -x) x 2 , j = 1, . . . , r (3.22)
which is non negative when x j 0 or x j x. Moreover, -1 φ(x j ) < 0 when x j ∈ (0, x), so that additional constraint (3.9) can be replaced by the following continuous one:

φ(x j ) 0, j = 1, . . . , r. (3.23)
Therefore, an unconstrained NLP can be easily defined, by introducing additional continuous variables s ∈ R r , such that x j := s 2 j for all j = 1, . . . , r and considering the resulting objective function, adjoint with three penalty terms, one replacing each set of constraints:

x Σx + ρ(1 -1 r x) 2 + ρ µ x R -1 2 +τ r j=1 κ j (x j ) 2 , ( 3.24)
where ρ and τ are suitable positive parameters and κ j (x j ) := min{0, φ(x j )} for all j = 1, . . . , r.

A similar approach is stated also for the round lot purchasing constraints (3.11) or (3.12), that can be replaced by the following constraints:

κ j (x j ) = Bx j q j - Bx j q j 1 - Bx j q j - Bx j q j = 0, j = 1, . . . , r (3.25
) where v denotes the integer part of v ∈ R.

However, round lot purchasing constraints (3.11) might make impossible satisfy at the same time request (3.1c): consequently, the following new quadratic risk measure [START_REF] Mitchell | Rebalancing an Investment Portfolio in the Presence of Convex Transaction Costs and Market Impact Costs[END_REF] is considered:

x Σx (1 r x) 2 , ( 3.26)
leading to an alternative definition of (3.24):

x Σx

(1 r x) 2 + ρ min{0, 1 -1 r x} 2 + ρ µ x R -1 2 +τ r j=1 κ j (x j ) 2 . (3.27)
Bartholomew-Biggs and Kane [START_REF] Bartholomew-Biggs | A Global Optimization Problem in Portfolio Selection[END_REF] apply a DIRECT (DIviding RECTangles) type global algorithm (see [START_REF] Finkel | DIRECT Optimization Algorithm User Guide[END_REF][START_REF] Gablonsky | DIRECT Version 2.0[END_REF][START_REF] Gablonsky | Modifications of the DIRECT Algorithm[END_REF][START_REF] Jones | The DIRECT Global Optimization Algorithm[END_REF][START_REF] Jones | Lipschitzian Optimization without the Lipschitz Constant[END_REF]) to the previous unconstrained NLPs (3.24) and (3.27).

. . Balanced Objective Functions

Mean-variance portfolio selection problems (3.1) and (3.4) are naturally multi-objective optimization programs since usually investors want to gain the maximum profit at the minimum risk: these are, of course, conflicting targets, that have to be considered at the same time.

Several authors (see, e.g., [START_REF] Chang | Heuristics for cardinality constrained portfolio optimization[END_REF]) use standard (linear) scalarization techniques such as the Weighted Sum approach (see, for example, [START_REF] Ehrgott | Multicriteria Optimization[END_REF]). Namely, they consider the "balanced" objective function

λ(x Σx) -(1 -λ)(µ x), (3.28)
where λ ∈ [0, 1] is a parameter which represents investor's risk aversion. Let θ 1 , θ 2 ∈ R + be two parameters, a more general variant is proposed by Schaerf [START_REF] Schaerf | Local Search Techniques for Constrained Portfolio Selection Problems[END_REF]:

θ 1 (x Σx) + θ 2 max{0, µ x -R}. (3.29)
Bertsimas and Shioda [START_REF] Bertsimas | Algorithm for Cardinality-Constrained Quadratic Optimization[END_REF] (see also [START_REF] Bertsimas | Portfolio Construction Through Mixed-Integer Programming at Grantham, Mayo, Van Otterloo and Company[END_REF]) introduce an extended "balanced" objective function, considering also trading and turnover requirements with respect to a given initial portfolio x (0) :

1 2 x Σx -µ x + r j=1 ι j x j -x (0) j 2 , ( 3.30)
where ι j > 0 is a coefficient for asset j which models the symmetric purchasing/selling impact with respect to the stock price. Finally, Tadonki and Vial [START_REF] Tadonki | Portfolio Selection with Cardinality and Bound Constraints[END_REF] and Shaw et al. [START_REF] Shaw | Lagrangian Relaxation Procedure for Cardinality-Constrained Portfolio Optimization[END_REF] consider respectively constant and linear transaction costs embedded in a quadratic "balanced" objective function, namely respectively

λ 1 (x Σx) -λ 2 (µ x) + c δ (3.31) λ 1 (x Σx) -λ 2 (µ x) + c x, (3.32)
where λ 1 ∈ R + and λ 2 ∈ R + are two positive scalars, c ∈ R r + is a vector, whose entries represent the transaction costs for the portfolio assets and δ ∈ R r is the binary vector defined in constraints (3.10).

.

In this section we present several different possible reformulations and approximations for the mean-variance portfolio optimization problem.

. . SOCC Inner Approximations

As observed in Section 3.3.2, given a probability distribution on the portfolio returns, it is not always possible to write the problem (3.7) in a closed form: the exact value for the quantile F -1 (x) (1 -p) is known only for special distributions (e.g., normal distribution, Student distribution, uniform distribution on an ellipsoid). However, if the probability distribution of the expected returns is only partially known, the value of its quantiles can be approximately computed using several well-known probability inequalities [START_REF] Lejeune | Portfolio Optimization with Combinatorial and Downside Return Constraints[END_REF], such as, e.g., Cantelli [START_REF] Bonami | An Exact Solution Approach for Portfolio Optimization Problems under Stochastic and Integer Constraints[END_REF], Chebyshev [START_REF] Bonami | An Exact Solution Approach for Portfolio Optimization Problems under Stochastic and Integer Constraints[END_REF], and Camp-Meidell [START_REF] Lejeune | A VaR Black-Litterman Model for the Construction of Absolute Return Fund-of-Funds[END_REF] inequalities (see also [START_REF] Hardy | Inequalities[END_REF][START_REF] Ion | Nonparametric Statistical Process Control[END_REF][START_REF] Lin | Probability Inequalities[END_REF][START_REF] Marshall | Inequalities: Theory of Majorization and Its Applications[END_REF]). Theorem 3.6.1. (Bonami and Lejeune [START_REF] Bonami | An Exact Solution Approach for Portfolio Optimization Problems under Stochastic and Integer Constraints[END_REF]) Assume the first and the second moment of the probability distribution of the portfolio return are finite. The SOCC

µ x - p 1 -p √ x Σx R (3.33)
is an approximation of the chance constraint (3.4b).

Theorem 3.6.2. (Bonami and Lejeune [START_REF] Bonami | An Exact Solution Approach for Portfolio Optimization Problems under Stochastic and Integer Constraints[END_REF]) Assume the first and the second moment of the probability distribution of the portfolio return are finite and the distribution is symmetric. The SOCC

µ x - 1 2 (1 -p) √ x Σx R (3.34)
is an approximation of the chance constraint (3.4b).

Theorem 3.6.3. (Lejeune [START_REF] Lejeune | A VaR Black-Litterman Model for the Construction of Absolute Return Fund-of-Funds[END_REF]) Assume the first and the second moment of the probability distribution of the portfolio return are finite and the distribution is symmetric and unimodal. The SOCC

µ x - 1 9 (1 -p) √ x Σx R (3.35)
is an approximation of the chance constraint (3.4b).

Remark 3.6.4. (Bonami and Lejeune [START_REF] Bonami | An Exact Solution Approach for Portfolio Optimization Problems under Stochastic and Integer Constraints[END_REF]) The approximation given by Theorem 3.6.2 for a symmetric probability distribution is tighter than the one given by Theorem 3.6.1 and that the approximation given by Theorem 3.6.3 for a symmetric unimodal probability distribution is tighter than the one given by Theorem 3.6.2.

. . Variance Reformulation

Given the symmetric positive definite matrix Σ, we consider its Cholesky decomposition Σ = CC , where C ∈ R r×r is a lower triangular matrix. From a computational viewpoint, the Cholesky decomposition is twice faster and more stable than LU factorization or Gauss elimination method (see [START_REF] Kielbasinski | A Note on Rounding-Error Analysis of Cholesky Factorization[END_REF][START_REF] Meinguet | Refined Error Analyses of Cholesky Factorization[END_REF][START_REF] Sun | Rounding-Error and Perturbation Bounds for the Cholesky and LDL T Factorizations[END_REF]) and it is implemented in High Performance Computing numerical software libraries (see [START_REF] Anderson | LAPACK User's Guide[END_REF][START_REF] Blackford | ScaLAPACK User's Guide[END_REF][START_REF] Dongarra | LINPACK User's Guide[END_REF]).

Note that Cholesky decomposition exists and is unique if matrix Σ is positive definite (see [START_REF] Golub | Matrix Computation. Third. Johns Hopkins Studies in the Mathematical Science[END_REF][START_REF] Heath | Diophantus of Alexandria[END_REF]) and this property is verified by variancecovariance matrix, if we exclude the case of exact collinearity of the random variables, i.e., we assume that none of the risky asset can be exactly replicated by a linear combination of the other ones. The hypotheses to apply Cholesky decomposition to positive semidefinitive matrices are identified in [START_REF] Higham | Analysis of the Cholesky Decomposition of a Semi-Definite Matrix[END_REF][START_REF] Heath | Diophantus of Alexandria[END_REF][START_REF] Moler | On the Householder-Fox Algorithm for Decomposing a Projection[END_REF] and error analysis is instead formally stated in [START_REF] Moler | On the Householder-Fox Algorithm for Decomposing a Projection[END_REF] for idempotent matrices and in [START_REF] Higham | Analysis of the Cholesky Decomposition of a Semi-Definite Matrix[END_REF] for the general case.

By assuming positive definiteness for covariance matrix Σ and introducing non negative decision variable h ∈ R + , we obtain the following problem, equivalent to (3.7):

min x,h C x 2 2 (3.36a) s.t. µ x -R h (3.36b) F -1 (x) (1 -p) C x 2 -h (3.36c) 1 r x = 1 (3.36d)
x 0, h 0. (3.37)

. .

Period-separable Reformulation

As pointed out by Filomena and Lejeune [START_REF] Filomena | Stochastic Portfolio Optimization with Proportional Transaction Costs: Convex Reformulations and Computational Experiments[END_REF], the variance of the portfolio can be reformulated as the Euclidean norm of a vector, whose number of components T corresponds to the number of data points, by using the following preliminary result: Theorem 3.6.6. (Konno and Suzuki [START_REF] Konno | A Fast Algorithm for Solving Large Scale Mean-Variance Models by Compact Factorization of Covariance Matrices[END_REF]) Let ν jt be the (observed) return of asset j at time t and introduce the extra variables b t = r j=1 (ν jt -µ j )x j (t = 1, . . . , T ). The variance of the portfolio return can be rewritten as

x Σx = 1 T b 2 2 .
Therefore, the probabilistic Markowitz portfolio model (3.4) can be reformulated as the following convex NLP:

min x,h,b 1 T b 2 2 (3.38a) s.t. µ x -R h (3.38b) F -1 (x) (1 -p) √ T b 2 -h (3.38c) b t - r j=1 (ν jt -µ j ) x j = 0, t = 1, . . . , T (3.38d
)

1 r x = 1 (3.38e)
x 0, h 0.

Furthermore, Filomena and Lejeune [START_REF] Filomena | Stochastic Portfolio Optimization with Proportional Transaction Costs: Convex Reformulations and Computational Experiments[END_REF] observe that in order to mathematically compute the variance in problems (3.4) and (3.37) the estimate of only r(r + 1)/2 covariance terms is needed: this situation can lead to several coherence problems for the covariance matrix (see Section 3.1). Moreover, the approach described in this section does not require any assumption on matrix Σ. Finally, we can consider the corresponding equivalent epigraph formulation of problem (3. . Mean-variance portfolio selection problem with the constraints introduced in Section 3.4 gives rise to a convex MIQP, which is at least as difficult as NP-hard, because it includes MILP as special case [START_REF] Kannan | On the Computational Complexity of Integer Programming Problems[END_REF][START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF]. Nowadays, MIQPs can be solved via commercial and open-source solvers (see, e.g., [START_REF] Bonami | Algorithms and Software for Convex Mixed Integer Nonlinear Programs[END_REF][START_REF] Lodi | Mixed Integer Nonlinear Programming Tools: An Updated Practical Overview[END_REF][START_REF] Vigerske | MINLP Solver Software[END_REF]). In this section we overview specialized and more efficient computational procedures recently proposed in literature.

In [START_REF] Bienstock | Computational Study of a Family of Mixed-Integer Quadratic Programming Problems[END_REF] Bienstock proposes a tailored Branch-and-Cut (BC) procedure to solve the cardinality constrained portfolio problem, where (3.16) is replaced with the "surrogate" constraint r j=1

x j x j K.

(3.39)

Several types of cutting planes, namely mixed-integer rounding inequalities, knapsack cuts, intersection cuts, and disjunctive cuts are also considered in the same paper.

Bertsimas and Shioda [START_REF] Bertsimas | Algorithm for Cardinality-Constrained Quadratic Optimization[END_REF] develop a BC algorithm where at each node of Branch-and-Bound (BB) tree the convex continuous relaxation of problem (3.1) with cardinality (3.16) and buy-in (3.10) constraints is solved by means of Lemke's method [START_REF] Cottle | The Linear Complementarity Problem[END_REF]. The portfolio problem with objective function (3.31) and cardinality (3.16) and buy-in (3.10) constraints was solved by Tadonki and Vial [START_REF] Tadonki | Portfolio Selection with Cardinality and Bound Constraints[END_REF] with BB techniques together with a Bender decomposition approach.

Lee and Mitchell [START_REF] Lee | Computational Experience of An Interior-Point SQP Algorithm in a Parallel Branch-and-Bound Framework[END_REF] describe a parallel BB framework for the cardinality constrained portfolio selection problem, in which each node is approximated by means of Sequential Quadratic Programming (SQP) and each quadratic subproblem is solved via interior-point method (see, e.g., [START_REF] Nemirovskii | Interior-Point Methods for Optimization[END_REF][START_REF] Nesterov | Interior-Point Polynomial Algorithms in Convex Programming[END_REF]).

Frangioni and Gentile [START_REF] Frangioni | Perspective Cuts for a Class of Convex 0-1 Mixed Integer Programs[END_REF] solve problem (3.1) with minimum and maximum buy-in thresholds additional constraints (3.10) with a BC method improved by using Perspective Cuts (see also [START_REF] Frangioni | Approximated Perspective Relaxations: A Project and Lift Approach[END_REF]), a family of valid inequalities, related to the perspective function (see [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms[END_REF][START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms[END_REF]) and to the convex envelope of the objective function (see [START_REF] Günlük | Perspective Reformulation and Applications[END_REF]).

Zheng et al. [START_REF] Zheng | Successive Convex Approximations to Cardinalityconstrained Convex Programs: A Piecewise-linear DC Approach[END_REF] propose a difference of convex functions approach to the cardinality constrained quadratic program, by replacing the cardinality constraint (3.15) with the following piecewise linear approximation:

1 ω x 1 - r j=1 max{x j -ω, 0} + max{-x j -ω, 0} 0, (3.40)
where ω > 0 is a given parameter. Non-smooth approximation (3.1) with constraint (3.40) is solved by means of Successive Convex Approximation method. This algorithm determines a Karush-Kuhn-Tucker (KKT) point or defines a sequence of points converging to a KKT point for the ω-parametrized approximation. Moreover, the authors show that, letting ω → 0 + , the optimal value of the approximate problem approaches the optimal value of the original problem.

Shaw et al. [START_REF] Shaw | Lagrangian Relaxation Procedure for Cardinality-Constrained Portfolio Optimization[END_REF] solve cardinality constrained portfolio problem under the assumption that vector µ of assets returns can be decomposed according to a multiple factor model [START_REF] Castellani | Manuale di Finanza[END_REF], i.e., µ = Ξf + u, where r represents the number of different factors, Ξ ∈ R r×r is the sensitivity-factor matrix, f ∈ R r is the factor-return vector, and u ∈ R r is the asset-specific (non-factor) returns vector. A Lagrangian relaxation of the problem is then solved by means of sub-gradient procedure [START_REF] Rockafellar | Theory of Subgradients and Its Applications to Problems of Optimization: Convex and Nonconvex Functions[END_REF] and embedded in a BB framework.

In [START_REF] Bonami | An Exact Solution Approach for Portfolio Optimization Problems under Stochastic and Integer Constraints[END_REF] Bonami and Lejeune deal with the deterministic equivalent (3.7) of the probabilistic portfolio selection problem with buy-in threshold (3.10), round lot purchasing (3.12), and diversification (3.13)-(3.14) constraints, by proposing a Nonlinear BB algorithm [START_REF] Belotti | Mixed Integer Nonlinear Optimization[END_REF] with two specific branching rules:

1. Idiosyncratic Risk Branching, consisting in selecting the fractional variable δ j or γ j , which corresponds to the asset with the highest expected return;

2. Portfolio Risk Branching: consisting in selecting the fractional variable δ j or γ j , whose integer fixing has the highest impact on the objective function (3.1a).

Buchheim et al. [START_REF] Buchheim | A Frank-Wolfe Based Branch-and-Bound Algorithm for Mixed-Integer Portfolio Optimization[END_REF] consider portfolio selection problem with objective function (3.28), constraints (3.1d) and (3.3) and integrality requirement on the decision variables, i.e.,

x ∈ Z r , (3.41)
which represents the units of assets held in the investor's portfolio. They introduce a new BB algorithm where the continuous relaxation is solved through an efficient Frank-Wolfe type method with non-monotone Armijo line-search.

Burdakov et al. [START_REF] Burdakov | Mathematical Programs with Cardinality Constraints: Reformulation by Complementarity-type Conditions and a Regularization Method[END_REF] deal with the cardinality constrained portfolio problem, by introducing a NLP reformulation, whose global minima are the same of the ones of the original problem. The NLP is solved via a sequence of regularized programs (see [START_REF] Kanzow | A New Regularization Method for Mathematical Programs with Complementarity Constraints with Strong Convergence Properties[END_REF]).

We end this section with Table 1 that summarizes the main characteristics of the papers described above. In particular, the columns report the authors, the year of publication of the paper, the objective function and constraints of the tackled problem, the proposed algorithm, the competitors employed as benchmarks, and the instances that were used for the computational experiments. Increasing Set [START_REF] Gaspero | Hybrid Metaheuristics for Constrained Portfolio Selection Problems[END_REF], [START_REF] Moral-Escudero | Selection of Optimal Investment Portfolios with Cardinality Constraints[END_REF], [START_REF] Chang | Heuristics for cardinality constrained portfolio optimization[END_REF] [216], [START_REF] Schaerf | Local Search Techniques for Constrained Portfolio Selection Problems[END_REF] 5 additional data sets 4

Table 1:

Deterministic exact approaches to mean-variance portfolio selection problem (see also [START_REF] Mansini | Twenty Years of Linear Programming Based Portfolio Optimization[END_REF] and references within). References are sorted in chronological order; papers published in the same year are sorted according to alphabetic order of the last name of the corresponding first author.

.

The portfolio problem we are interested in can be formulated as follows. Let x ∈ R r + be the decision variables representing the fraction of asset j ∈ {1, . . . , r} held in the portfolio. Each asset j is characterized by a mean expected return µ j (j ∈ {1, . . . , r}) and a non negative, possibly non-convex non-concave, cost function C j (x j ) : R + → R + (j ∈ {1, . . . , r}). The goal consists in maximizing the overall return of the portfolio in presence of cardinality constraint, risk constraint, and no-short-selling constraint. The risk of the portfolio is represented by its total variance x Σx, where Σ ∈ R r×r is the covariance matrix of the return distribution of the assets. The portfolio problem reads:

max µ x -C(x) (3.42a) s.t. 1 r x = 1 (3.42b)
x Σx σ (3.42c)

x 0 K (3.42d) x ∈ [x, x] ∪ {0}, (3.42e)
where C(x) := r j=1 C j (x j ), and σ and K are the maximum levels of risk and sparsity we aim for our portfolio, respectively. Remark 3.8.1. We assume both the mean vector µ and the covariance matrix Σ are perfectly known or coherently estimated from historical data measurements. Problem (3.42) generalizes the usual portfolio problem with concave transaction costs already considered in the literature (see, e.g., Konno and Wijayanayake [START_REF] Konno | Portfolio Optimization Problem under Concave Transaction Costs and Minimal Transaction Unit Constraints[END_REF] and Xue et al. [START_REF] Xue | Mean-Variance Portfolio Optimal Problem under Concave Transaction Cost[END_REF]). To the best of our knowledge, this is the first time this kind of portfolio problems is addressed in its entirely.

The previous MVPS problem is NP-hard (see Section 3.4.4) and can be exactly reformulated, in the sense explained in [START_REF] Liberti | Reformulations in Mathematical Programming: Definitions and Systematics[END_REF], by introducing extra binary variables δ ∈ {0, 1} r as follows:

max µ x -C(x) (3.43a) s.t. 1 r x = 1 (3.43b)
x Σx σ (3.43c)

1 r δ K (3.43d) δ x x δ x (3.43e) δ ∈ {0, 1} r . (3.43f)

. . Pointwise Reformulation

A pointwise reformulation for problem (3.43) can be obtained by replacing each function C j (x j ) with the term (1 + x j ) θ j parametrized by θ j . The only part of the problem, where the substituting terms appear, is the objective function:

µ x - r j=1 (1 + x j ) θ j = (µ -θ) x -1 r θ, (3.44)
which becomes an affine function of the continuous decision variables x j (j ∈ {1, . . . , r}). Thus, the pointwise reformulation is:

-1 r θ + max (µ -θ) x (3.45a) s.t. 1 r x = 1 (3.45b)
x Σx σ (3.45c)

1 r δ K (3.45d) δ x x δ x. (3.45e)
Remark 3.8.2. The pointwise reformulation (3.45) is spanning: the replacement terms θ j perfectly matches the replaced terms

C j (x j )
x j +1 at each feasible point x j (j ∈ {1, . . . , r}). By Lemma 2.2.6, there exist values of θ which make (3.45) a bounding reformulation for the original problem (3.42).

Example 3.8.3. The pointwise reformulation (3.45) is not exact. For instance, we consider a portfolio problem with r = 2 stocks, Σ = I 2 , K = 2, x = 0, x = 0.55 and σ = r, i.e., an essentially unconstrained portfolio problem, apart from the budget constrained

x 1 + x 2 = 1. Let µ = (1, 0.40) and C(x) = (C 1 (x 1 ), C 2 (x 2 )) = (x 1 , 0)
. For all j ∈ {1, 2}, θ j = C j 1+x j and, hence, Θ j = [C j (x)/(1 + x), C j (x)/(1 + x)]. For the example, we have Θ 1 = [0, 0.55] and Θ 2 = {0}. The objective function of the original formulation (3.42) is simply max(0.40 x 2 ) and its global optimum is x * = (0.45, 0.55) . The objective function of the pointwise reformulation is max((1 -θ 1 ) x 1 + 0.50 x 2 -θ 1 ). Now, x * = (0.45, 0.55) is a global optimum of the pointwise reformulation if and only if (1 -θ 1 ) 0.40, i.e., if and only if θ 1 > 0.60: however, these values of θ 1 do not belong to the set Θ 1 . Remark 3.8.4. The pointwise reformulation (3.45) is not efficient (see Section 3.4.4). However, since convex MIQPs can be solved nowadays with reasonable efficiency, certainly more efficiently than non-convex non-concave MINLPs, the proposed reformulation is good.

The previous approach can be easily extended to cardinality constrained portfolio selection problem with cost functions and fixed transaction costs, i.e., when a term r j=1 c j δ j is added to the objective function: in these problems the investor has to pay a fixed amount c j (j ∈ {1, . . . , r}) of money if he/she decided to buy a certain asset j (j ∈ {1, . . . , r}). Moreover, since the previous strategy consists essentially in substituting each single term in the separable function with an affine function depending on a given parameter, it can be extended to general MINLPs with separable non-convexities and non-concavities.

. . Computing MWU Costs/Gains

Since the substituting terms appear only in the objective function, we do not need to address feasibility issues in computing MWU costs/gains. In particular, we define, at each iteration t T , a r-dimensional vector α t ∈ [-1, 1] r whose components are:

α j,t =
C j (x j,t ) -(x j,t + 1) θ j,t max(|C j (x j,t )|, |(x j,t + 1) θ j,t |) (j ∈ {1, . . . , r} , t T ).

(3.46)

Remark 3.8.5. Each component α j,t represents the scaled cost/gain determined by each asset j (j ∈ {1, . . . , r}). In other words, each α j,t takes into account the contribution of the fraction of each asset held in the portfolio to the overall cost function.

Remark 3.8.6. We define α t to be a vector, instead of a scalar as described in Section 2.4.3, because this definition is more effective in presence of separable nonconvexities and non-concavities: it allows us to better follow the numerical behavior of each single cost function C j (x j ) for all j ∈ {1, . . . , r}

In order to define the costs/gains, we simply set:

ψ j,t := α j,t (j ∈ {1, . . . , r} , t T ). (3.47)
. .

Computational Experiments

The test-bed set with respect to which we analyze our algorithm and the benchmarks consists in the 20 real-world instances described by Chang et al. [START_REF] Chang | Heuristics for cardinality constrained portfolio optimization[END_REF] publicly available through the OR-Library (see Beasley [START_REF] Beasley | Obtaining Test Problems via Internet[END_REF] and Beasley [START_REF] Beasley | OR-Library: Distributing Test Problems by Electronic Mail[END_REF]) on the web site http://www.brunel.ac.uk/~mastjjb/jeb/info. html. Each instance is characterized by the number n of assets and the value of the risk level σ. We impose, as in Chang et al. [START_REF] Chang | Heuristics for cardinality constrained portfolio optimization[END_REF], x = 0.01, x = 1 and K = 10.

Transaction Cost Functions

One of the aims of the computational experiments, in addition to analyze the behavior of the MWU algorithm against other methods to solve MVPS problems, consists in empirically evaluating the performances of the MWU algorithm depending from the nonlinearity of the replaced transaction cost functions. In particular, we consider the following five univariate functions (see Figure 1):

(a) C j (x j ) = -µ j ln 20-0.06(1+x j ) 1+x j
for all j ∈ {1, . . . , r}: this cost function is increasing, concave and "almost linear".

(b) C j (x j ) = -µ j ln 0.2-0.01(0.00001+x j ) 0.00001+x j for all j ∈ {1, . . . , r}: this cost function is increasing, concave, and replicates the behavior of the transaction cost function described in Konno and Wijayanayake [START_REF] Konno | Portfolio Optimization Problem under Concave Transaction Costs and Minimal Transaction Unit Constraints[END_REF].

(c) C j (x j ) = µ j (4x j + 0.12 sin(40x j )) for all j ∈ {1, . . . , r}: this cost function has a sinusoidal behavior similar to a step function.

(d) C j (x j ) = µ j (4x j + 0.3 sin(40x j )) for all j ∈ {1, . . . , r}: this cost function is similar to the one in (c) but with a "stronger nonlinear behavior".

(e) C j (x j ) = µ j (0.5x j + sin(50x j )) for all j ∈ {1, . . . , r}: this is the "most nonlinear" transaction cost function among which we tested the methods.

Remark 3.8.7. We want to emphasize that the definition of nonlinearity of a given transaction cost function was not given in rigorous mathematical terms, but it is mostly qualitative: in fact, we carefully use quotes around each nonlinearity characterization. However, since our aim consists in empirically analyzing the behavior of the MWU algorithm with respect to a given cost function, we trust most readers will agree with our categorization, by inspection of Figure 1.

Computational Environment

Since for the MVPS problem we consider in this section there are no tailored exact or heuristic methods to solve it and since the MWU algorithm is essentially a MS algorithm with a special strategy to choose the more promising initial points both in terms of feasibility and optimality, we compare the MWU with the MS.

We used T = 20 iterations for both MWU and MS. We adopted Bonmin [START_REF] Bonmin | [END_REF] as local MINLP solver: since Bonmin exactly solve convex MINLPs, it is a reliable heuristic for non-convex MINLPs. In particular, we employ Bonmin's native Branch-and-Bound (B-BB) algorithm (see Bonami et al. [START_REF] Bonami | More Branch-and-Bound experiments in convex nonlinear integer programming[END_REF] and Gupta and Ravindran [START_REF] Gupta | Branch-and-Bound Experiments in Convex Nonlinear Integer Programming[END_REF]), since it is generally more stable for nonconvex MINLPs. We used Cplex [START_REF]ILOG CPLEX 12.2 User's Manual[END_REF] as the convex MIQP solver for the pointwise reformulation (3.45), with a 600 seconds time limit, using only one thread. All of the computational experiments were performed on an Intel Xeon CPU E5649, 2.53 GHz, using only one processor.

Localization of the MS Subsolver

Since we want to compare MWU against MS and MS is essentially a procedure composed by two steps: a random choice of a starting point, and a local descend method, we would like that Bonmin BB-B behaves like a local solver. From preliminary tests, however, the behavior of Bonmin BB-B was mostly similar to a global solver: Bonmin BB-B sets a cut-off value for the optimum based on the starting point and for our test-bed it found almost always the same solution point. In order to turn Bonmin B-BB into a truly local solver, we consider an adding local branching constraint (see Fischetti and Lodi [START_REF] Fischetti | Local Branching[END_REF]) for the original formulation, which basically defines an upper bound νn , where ν ∈ [0, 1], on the number of flips of binary variables δ: where δ is the starting point. The spirit of the constraints (3.48) consists in enforcing a local exploration in the combinatorial neighborhood of the starting point δ . After several computational experiments we decided to set ν = 0.96, since lower values made the instance infeasible excessively often.

i r δ j =0 δ j + j r δ j =1 (1 -δ j ) ν n , ( 3
Tables 23456report the computational results for each transaction cost function. Their columns are as follows:

• instance name;

• maximum risk level σ;

• number r of assets quoted on the financial market;

• objective value for the MWU algorithm;

• CPU time (in seconds) for the MWU algorithm;

• objective value for the MWU algorithm with the local branching constraint;

• CPU time (in seconds) for the MWU algorithm with the local branching constraint;

• objective value for the MS algorithm with the local branching constraint;

• CPU time (in seconds) for the MS algorithm with the local branching constraint;

• relative objective value improvement from MS to MWU computed as

Γ = val(MWU) -val(MS) |val(MS)| ; (3.49)
• time improvement ratio from MS to MWU:

Λ = cpu(MS) cpu(MWU) ; (3.50)
• relative objective value improvement from MS to MWU with the local branching constraint (see Equation (3.49));

• time improvement ratio Λ from MS to MWU with the local branching constraint (see Equation (3.50)).

Computational Results

The comparison metrics are summarized in the last three lines with the sum ( ), average (avg), and the standard deviation (std) across all 20 instances. For the CPU time we reported also the geometrical mean among all the instances. While with transaction costs (a) and (b) MS performs better than MWU, with transaction costs (c) and (e) MWU overcomes MS: we can reasonably imply that MWU behaves better with respect to "high nonlinear" transaction cost functions. In particular, for the geometric average values, the relative improvement we obtain is considerably high.

A secondary observation is about the number of assets, which composed the optimal portfolios. In MVPS problems, defining a portfolio with few assets could be a secondary goal, since for each asset we have several costs, such as transaction costs, monitoring costs, and brokerage fee (see Di Lorenzo et al. [START_REF] Lorenzo | A Concave Optimization-based Approach to Sparse Portfolio Selection[END_REF]). MWU algorithm produces small values of the number of assets compared to MS algorithm.

Finally, for the geometric mean for the CPU time, the MS defeats MWU only for the cost function (a), in all the other cases MWU overcomes MS.

. In this chapter we dealt with the Mean-Variance Portfolio Selection problem. At the beginning, we have illustrated the mathematical models proposed in the literature with a survey about the possible objective functions and constraints. Then, we have adapted the MWU algorithm to the realworld portfolio problem with separable transaction cost with respect to the strategy to define the pointwise reformulation and to compute the costs/gains necessary for the algorithm. Computational experiments on real-world instances allow us to observe that the MWU algorithm performs better than the MS algorithm for "heavily nonlinear" instances.

M U LT I P L E N O N L I N E A R K N A P S A C

K P R O B L E M S .
Let x be an m × n array of non negative real variables x = [x ij ] (i = 1, . . . , m, j = 1, . . . , n) and define M = {1, . . . , m} and N = {1, . . . , n}. We consider a multiple nonlinear knapsack problem in which

• the objective function and the capacity constraints are expressed by separable, continuously differentiable functions f j (x ij ) and g j (x ij ) (i ∈ M, j ∈ N);

• the values of f and g do not depend on i, i.e., f j (x ij ) = f j (x kj ) and g j (x ij ) = g j (x kj ) when x ij = x kj for j ∈ N and i, k ∈ M;

• f j (x ij ) and g j (x ij ) are nonlinear non negative non-decreasing functions for j ∈ N and i ∈ M;

• for each j ∈ N, the total value of x ij over all i ∈ M cannot exceed a given upper bound u j ;

• integrality requirements may be imposed on part of the variables.

Remark 4.1.1. Note that there is no further assumption on f j (x ij ) and g j (x ij ) which, in general, can be non-convex and non-concave.

The Multiple NonLinear (Separable) Knapsack Problem (MNLKP) is:

max i∈M j∈N f j (x ij ) (4.1a) s.t. j∈N g j (x ij ) c i i ∈ M (4.1b) i∈M x ij u j j ∈ N (4.1c) x ij 0 i ∈ M, j ∈ N (4.1d) x ij integer i ∈ M, j ∈ N ⊆ N, (4.1e)
which can be informally described as follows. We are given m knapsacks and n items. Each item j has a profit function f j (x ij ) and a weight function g j (x ij ) (i ∈ M), and each knapsack i has a capacity c i . For each item j we want to assign x ij quantities (some restricted to integer values) to the knapsacks so that • the overall assigned profit is maximized, see (4.1a);

59

• for each knapsack i the overall assigned weight does not exceed the corresponding capacity, see (4.1b);

• for each item j the overall assigned quantity does not exceed the corresponding upper bound, see (4.1c).

The MNLKP generalizes the classical 0-1 Multiple Linear Knapsack Problem (MLKP) (see, e.g., Martello and Toth [START_REF] Martello | Knapsack Problems: Algorithms and Computer Implementations[END_REF] and Kellerer et a. [START_REF] Kellerer | Knapsack Problems[END_REF]): in the MLKP x ij are binary decision variables, i.e., x ij ∈ {0, 1} for all i ∈ M and j ∈ N, and the profit and the weight functions are linear, i.e., f j (x ij ) = p j x ij and g j (x ij ) = w j x ij , and u j = 1 for all j ∈ N.

It follows that MNLKP is, at least, strongly NP-hard. Moreover, the MN-LKP generalizes also the (single) NonLinear Knapsack Problem (NLKP) [START_REF] Martello | Heuristic Algorithms for the General Nonlinear Separable Knapsack Problem[END_REF]: MNLKP reduces to NLKP when m = 1 and, consequently, objective function (4.1a) and constraints (4.1b) read j∈N f j (x ij ) and j∈N g j (x ij ) c 1 , respectively.

The nonlinear knapsack structure arises in many different real-world problems, such as portfolio selection, capacity and production planning, and resource allocation (see, e.g., Ibaraki and Katoh [START_REF] Ibaraki | Resource Allocation Problems[END_REF], Bretthauer and Shetty [START_REF] Bretthauer | The Nonlinear Knapsack Problem: Algorithms and Applications[END_REF], and Li and Sun [START_REF] Li | Nonlinear Integer Programming[END_REF]). For instance, we assume we have m different economical resources and n products and we want to subdivide a certain amount of advertising budget c i related to resource i in order to maximize the overall expected sales for all resources. Obviously, the profit is increasing with the advertising investment since more and more buyers happened to find out our advertisement. However, at some point, a saturation effect occurs when, despite we further increase the investments, no more buyers are interested in our products. In this example the profit function happens to be non-convex non-concave with a shape represented in Figure 3 where the advertisement cost could be linear, i.e., g j (x ij ) = x ij , if a constant unit cost is assumed or nonlinear if economies of scale are considered, i.e., unit costs decrease with size.

To the best of our knowledge no tailored exact methods or heuristics have been proposed for the MNLKP. Zhang and Hua [START_REF] Zhang | A Unified Method for a Class of Convex Separable Nonlinear Knapsack Problems[END_REF] proposed an exact method for the minimization version of the convex continuous NLKP, i.e., when all the profit and weight functions are convex, and all variables are continuous, i.e., N = ∅. Zhang and Chen [START_REF] Zhang | Heuristic and Exact Solution Method for Convex Nonlinear Knapsack Problem[END_REF] described exact and heuristic methods for the pure integer version of same problem, i.e., when N = N.

The rest of the chapter consists of three other sections. In Section 4.2 we apply the MWU algorithm for MNLPs, by specializing its main steps, namely the construction of the pointwise reformulation and the definition of costs/gains, to the MNLKP. Since this method does not fit very well this kind of optimization problems, we propose also other heuristic procedures. First of all we introduce the surrogate and Lagrangian relaxations for MN-LKPs in Section 4.3, while in Section 4.4 we discuss a constructive solution approach, whose main aspects are a greedy heuristic, two other heuristics based on the feasibility recovery of the solution produced by the surrogate relaxation, and a local search procedure to improve the quality of the heuristic solution. Extensive computational experiments are conducted both for the MWU heuristic and for the constructive one, with respect to challenging instances.

. . . Pointwise Reformulation

In order to define a pointwise reformulation for the MNLKP, we replace the "complicated" non-convex non-concave terms g j (x ij ) with affine terms in x ij , i.e., with terms θ ij x ij for all i ∈ M and j ∈ N:

max i∈M j∈N f j (x ij ) (4.2a) s.t. j∈N θ ij x ij c i i ∈ M (4.2b) i∈M x ij u j j ∈ N (4.2c) x ij 0 i ∈ M, j ∈ N (4.2d) x ij integer i ∈ M, j ∈ N ⊆ N, (4.2e)
Remark 4.2.1. The pointwise reformulation (4.2) is spanning, since the replacement terms θ ij correspond to the replaced terms

g j (x ij ) x ij
for all i ∈ M and j ∈ N. By Lemma 2.2.6, there exist values of θ which make the pointwise reformulation (4.2) a bounding reformulation for the original problem (4.1). Remark 4.2.2. The pointwise reformulation (4.2) is not efficient. Nevertheless, the feasible set of the reformulation (4.2) is polyhedral. We note that the problem (4.1) is always feasible since the zero solution, i.e., x ij = 0 for all i ∈ M and j ∈ N, is always a feasible point for the problem.

The approach, we have adopted in this section, is quite similar to the one we followed for the MVPS problem with univariate cost functions: essentially we substitute the nonlinear non-convex non-concave terms with affine terms in the decision variables.

. . Computing MWU Costs/Gains

Since the replaced terms appear only in the constraint (4.2b), we do not take into account optimality issues in defining MWU costs/gains. In particular, at each iteration t T , we implement the following strategy. We define a profit-to-weight ratio for each item x ij for i ∈ M and j ∈ N, as the ratio between its profit and weigh function values, representing the profitability of filling the knapsack i (i ∈ M) with x ij units of item j (j ∈ N), as follows:

r ij,t = f j (x ij,t ) g j (x ij,t ) (i ∈ M, j ∈ N, t T ). (4.3)
Then, we compute the feasibility costs/gains β i,t (i ∈ N, t T ) as the scaled difference from the left hand side and the right hand side of constrains (4.1b) calculated on the current point x t (see Section 2.4.3):

β i,t = max( j∈N g j (x ij,t ) -c i , 0) max s t (max( j∈N g j (x ij,s ) -c i ), 0) (i ∈ M, t T ). (4.4)
Finally, we need to spread the previous feasibility cost/gain to the item j (j ∈ N). We simply scale the feasibility cost/gain with respect to the profitto-weight ratio, defining in this way our MWU costs/gains for the MNLKP:

ψ ij,t = β i,t r ij,t j∈N r ij,t (i ∈ M, j ∈ N, t T ). (4.5)
We choose the profit-to-weight ratio since in heuristic methods for NLKP it is used to sort the items (see D'Ambrosio and Martello [START_REF] Martello | Heuristic Algorithms for the General Nonlinear Separable Knapsack Problem[END_REF]). The items with a greater ratio are more promising in terms of the trade-off represented by the objective function we want to maximize and the amount we have to pay in order to fill the knapsacks with those items. In the heuristic method proposed by D'Ambrosio and Martello [START_REF] Martello | Heuristic Algorithms for the General Nonlinear Separable Knapsack Problem[END_REF] for the NLKP, they select the item with the best profit-to-ratio at first and then fill as much as possible the knapsack with that item.

. . Computational Experiments

The MWU method was experimentally compared with the MS algorithm. We use Ipopt [126] with its default options: Ipopt is an exact solver for convex NLPs, hence it can be used as local solvers for non-convex NLPs. The number of iterations was set to T := 10 both for the MWU and MS. Moreover, we compare the MWU with Ipopt with one starting point and Couenne [59],

which is an open-source global solver for MINLPs.

For the profit and the weight functions, we used the ones proposed in D'Ambrosio and Martello [START_REF] Martello | Heuristic Algorithms for the General Nonlinear Separable Knapsack Problem[END_REF]. In particular, the profits were always obtained from

f j (x ij ) = c j 1 + b j e -a j (x ij +d j ) (4.6)
by uniformly randomly generating a j in [0.1, 0.2], b j and c j in [0, 100], and d j in [-100, 0]. These functions replicate the sigmoid shape of the profit function depicted in Figure 3. The upper bounds on the variables x ij were set to u j := 100 for all j ∈ N. For the weight, we adopted the concave increasing functions:

g j (x ij ) = p j x ij + q j - √ q j . (4.7)
Moreover, we define two sets of instances depending on the way we generated the capacities. In order to obtain challenging instances, we adopted the numerical methods described in Chapter 6 of Martello and Toth [START_REF] Martello | Knapsack Problems: Algorithms and Computer Implementations[END_REF]: (4.9) In both cases, the m-th capacity was set to:

Similar capacities c i uniformly random in 0.4 n j=1 g j (u j ) m , 0.6 n j=1 g j (u j ) m (i = 1, . . . , m -1), (4.8
c m = 0.5 n j=1 g j (u j ) - m-1 i=1 c i . ( 4.10)
The value of the number of items n varied in the set {10, 20, 50}, while the value of the number of knapsacks m varied in the set {2, 5, 10}. For each combination of (n, m) 20 real instances, i.e., with N = ∅ are generated. The total number of instance we tested was 360. All instances are available at http:// or.dei.unibo.it/library/multiple-nonlinear-knapsack-problem. If the pointwise reformulation was not solved within the time limit, we set its solution to zero.

All the experiments were performed on an Intel Xeon, CPU 3220, 2.4 GHz, using only one processor. The local method (Ipopt), the global algorithm (Couenne), and the MWU were run with a time limit of one CPU hour per problem.

Tables 7 and9 report the average solution values for the group of instances with similar and dissimilar capacities, respectively. The entries give:

• number of knapsacks;

• number of items;

• average values produced by MWU, Ipopt with a single starting point (Ipopt_1), MS, i.e., Ipopt with ten starting points (Ipopt_10), and Couenne;

Tables 8 and10 report the average CPU times (in seconds) for instances with similar and dissimilar capacities, respectively. Tables 7 and9 show that generally the solution produced by the MWU heuristic was better than the ones produced by Ipopt_1 and worse than the ones produced by Ipopt_10. For the bigger continuous similar instances, Table 7 shows the MWU algorithm was on average completely outperformed by Ipopt for all the instances with m = 5 and for the instances with m = 10 and n = 20. For the similar instances there is only one case in which the MWU algorithm is better than Ipopt_10: the dissimilar instances with m = 2 and n = 50.

Moreover, Tables 8 and10 clearly indicate the average CPU times for MWU was one order of magnitude larger than the ones of Ipopt_10 and around two orders of magnitude larger than the ones of Ipopt_1.

We point out all the knapsack problems were solved to optimality within the time limit.

From the previous observations, we can argue that the MWU method was in general outperformed by the heuristic methods available for the MNLKPs both in terms of quality of the solution found and of CPU times needed. The behavior of the proposed heuristic was not satisfying, and this is the reason why in the next sections we will introduce a different approach to heuristically solve MNLKPs relying on the discretization of the solution space and on the surrogate relaxation.

.

In this section we introduce some relevant relaxations of MNLKP.

Let (π 1 , . . . , π m ) be an m-dimensional vector of non negative multipliers. By multiplying the i-th constraints (4.1b) and summing up all the new ca-

f j (x ij ) (4.11a) s.t. i∈M π i j∈N g j (x ij ) i∈M π i c i (4.11b) i∈M x ij u j j ∈ N (4.11c) x ij 0 i ∈ M, j ∈ N (4.11d) x ij integer i ∈ M, j ∈ N ⊆ N. (4.11e)
Let val(S(MNLKP,π)) denote the optimal value of (4.11) under given multipliers π. The surrogate dual problem

min π 0 {val(S(MNLKP,π))} (4.12)
consists in finding the optimal vector of multipliers, i.e., the one producing the minimum optimal value for the surrogate relaxation, and hence the tighter upper bound for the MNLKP. Remark 4.3.1. While the surrogate relaxation of the 0-1 MLKP has strong duality property, i.e., the vector of multipliers which produced the minimum value for the surrogate relaxation is π i = k for all i ∈ M [START_REF] Martello | Knapsack Problems: Algorithms and Computer Implementations[END_REF], the same result does not hold for the MNLKP, as the following example shows. Let m = 2, n = 1, u 1 = 100, c 1 = 10, c 2 = 2, and, for i ∈ {1, 2}, f 1 (x i1 ) = x i1 , g 1 (x i1 ) = 80/ 1 + 50e -1 10 (x i1 -10) . For π 1 = π 2 = 1 the optimal solution to S(MNLKP,π) is x 11 = x 21 24 (with g 1 (x 11 ) = g 1 (x 21 )

6) and has value 48. For π 1 = 1 and π 2 = 2 such solution violates (4.11b) and the optimal solution is x 11 26 and x 21 18 (with g 1 (x 11 ) 7.2 and g 1 (x 21 ) 3.4), of value 44.

Even though the optimal surrogate multipliers are not known a priori, good multipliers can heuristically found (see Section 4.4.5) and, from the surrogate solution, a feasible solution can be easily defined (see Section 4.4.2).

Let (λ 1 , . . . , λ m ) be an m-dimensional vector of non negative multipliers, a possible Lagrangian relaxation, L(MNLKP,λ), of the MNLKP can be obtained by relaxing (4.1b): .

i∈M λ i c i + max i∈M j∈N f j (x ij ) -λ i g j (x ij ) (4.13a) s.t. i∈M x ij u j j ∈ N (4.13b) x ij 0 i ∈ M, j ∈ N (4.13c) x ij integer i ∈ M, j ∈ N ⊆ N. ( 4
In this section we describe two different kinds of heuristic algorithms for the MNLKP: the first type of heuristic is represented by a constructive procedure based on the discretization of the solution space; the second one, instead, is based on feasibility recovery strategies to restore the feasibility of the surrogate solutions with respect to the relaxed constraints.

. . Discretization Heuristic

The constructive procedure extends the heuristic method proposed by D'Ambrosio and Martello [START_REF] Martello | Heuristic Algorithms for the General Nonlinear Separable Knapsack Problem[END_REF] for the NLKP to the MNLKP.

We assume without loss of generality that the knapsacks are preliminary sorted in non-increasing order according to their capacities, i.e., c

1 c 2 • • • c m .
The algorithm is based on the discretization of the solution space. Let s be the number of sampling and δ j = u j /s for j ∈ N (or δ j = max(1, u j /s ) if j ∈ N) be the sampling step. We consider the profit-toweight ratio, meaning the ratio between the profit functions and the weight functions evaluated over the sampling points:

r jk = f j (kδ j ) g j (kδ j ) (j ∈ N, k = 1, . . . , s). (4.14)
Moreover, we assume without loss of generality that the items are sorted in non-increasing order according to their maximum profit-to-weight ratios µ j = arg max k=1,...,s {r jk } , i.e., r 1µ 1 r 2µ 2 • • • r nµ n .

We apply the same strategy as in [START_REF] Martello | Heuristic Algorithms for the General Nonlinear Separable Knapsack Problem[END_REF] by considering one single knapsack at each iteration. We take the first two items and we try to fill the knapsack as much as possible with the first one (see Procedure Construct(i)). We calculate the higher sampling point μ1 such that the ratio of the first item is greater than the ratio of the second one (Step 5), we fill the current knapsack with μ1 δ 1 units of the first item (Steps 6) and we update the remaining upper bound and the residual capacity (Step 7). Assume by the moment that the sampling points corresponding to µ 2 and µ 3 remain feasible. At the second iteration, the second item is taken into account and the knapsack is filled with μ2 δ 2 units, where μ2 is the analogue of μ1 for the second item and, again, the upper bound and the capacity are updated, and so on. If instead (Step 8) for at least one of the next two items, say 2 and 3, the sampling point corresponding to µ 2 or µ 3 is infeasible, an update of the µ values is performed on items 2, 3, . . . (and, consequently, the item order might change). Remark 4.4.1. Note that, whenever a partial solution x ij is determined, the ratios of all the unscanned items might be re-calculated, but it is not necessary to re-sort all the items: in fact, the algorithm only needs the items with the first and the second best ratio, which can be found in linear time.

The algorithm stops when only one element remains unscanned and, in this case, tries to fill the current knapsack with the last item as much as possible (Step 12). Algorithm 5 Procedure Construct(i).

1: ci := c i ; 2: ūj := u j ; 3: j := 1; 4: while j < n and ci > 0 do 5: μj := max{k : r jk r (j+1)µ j+1 , µ j k s}; 6:

x ij := μj δ j ; 7:

ūj := ūj -x ij , ci = ci -g j (x ij ); 8:
if (g j+1 (µ j+1 δ j+1 ) > ci or g j+2 (µ j+2 δ j+2 ) > ci ) then update µ for items j + 1, j + 2, . . . (and possibly update their order); x in := min(g -1 n ( ci ), ūn );

13: ūn = ūn -x in , ci = ci -g n (x in );

14: end if

We simply assume that the weight functions g j are strictly increasing and continuous, so that the inverse g -1 j exists (Step 12). If it is not the case, we consider the pseudo-inverse of g j with the largest value of the pre-image. In order to compute this value, we only need to calculate the zeros of the function g j (x ij )ci , which can be evaluated in a CPU time bounded by a constant independent from the instance size.

The algorithm can be improved through a refined search for μj . Once it has been obtained (at Step 5), the interval [ μj , μj+1 ] can be searched with a smaller sampling step and new, more precise, profit-to-weight ratios for item j can be computed. In this way a more precise point μj is obtained, and the process can be iterated by further decreasing the sampling step.

Steps 5-9 are iterated at most n times. At each time we have to determine the items with the first and the second best ratio, which can be calculated in O(n), and the (pseudo-)inverse of the function g j , which can be effectuated, as said before, in a computational time bounded by a constant. The main loop is executed up to n times (Step 4). If the refinement parameters, i.e., the number of sampling step and the number of the refinements are bounded by a constant (as usual in practice), the time complexity of Construct(i) is O(n 2 ).

Procedure Construct(i) is therefore executed for each knapsack i by considering, at each iteration, only those items whose quantities are still smaller than the upper bound (see Procedure Constructive). At the end, a greedy heuristic is applied to fill the knapsack as much as possible with the current item (Step 5). The overall time complexity of Procedure Constructive(i) is O(mn 2 ).

Algorithm 6 Procedure Constructive.

1: for i := 1 to m do Construct(i) {comment: optionally include the refined search};

2: for i := 1 to m do 3:

if ci > 0 then 4:

for j := 1 to n do {comment: increase x ij as much as possible} 5:

x ij := x ij + min(g -1 j ( ci ), ūj );

6: if j ∈ N then x ij := x ij ; 7:
ūj := u j -k∈M x kj , ci := c i -k∈N g j (x ik );

Surrogate Heuristics

In this section we introduce two different heuristic procedures based on the feasibility recovery of the solution of the surrogate relaxation (4.11). Let us first consider the problem of determining good surrogate multipliers π. A series of preliminary experiments was performed on the benchmark instances adopted for the computational experiments of Section 4.4.5, with (i) π i uniformly random in [0.0, 3.0] for all i ∈ M;

(ii) π i uniformly random in [0.8, 1.2] for all i ∈ M;

(iii) π i uniformly random in [0.9, 1.1] for all i ∈ M, and (iv) π i = 1 for all i ∈ M.

It turned out that the surrogate solutions produced by (i) were dominated by the other generations, those produced by (ii) and (iii) had about the same quality, and those produced by (iv) were, on average, clearly the best ones. Additional tests were performed using (easier) convex and concave objective functions, globally obtaining the same results. It was thus decided to always adopt option (iv). In Section 4.3 we have shown that identical multipliers (optimal solution of the surrogate dual for the linear case) are not necessary optimal for the nonlinear case. It is worth observing that they appear to be a good choice for such case too, at least for the objective functions we considered.

Let x * ij be the surrogate solution. We assume, as in Section 4.4.1, that the knapsacks are sorted in non-increasing order according to their capacities and the items are sorted in non-increasing order according to the profit-toweight ratios evaluated over the surrogate solution, i.e.,

r j = i∈M f j (x * ij ) i∈M g j (x * ij ) (j ∈ N). (4.15)
The first heuristic (see Procedure Surrogate-feas-1(x * )) starts with a zero solution -note that this solution is always feasible for the MNLKP -and is divided into two main phases.

Algorithm 7 Procedure Surrogate-feas-1(x * ). for j := 1 to n do 14:

rj :=            0 if ūj = 0
f j (min(x īj + ūj ,max(0, g -1 j ( cī +g j (x īj )) ))) g j (min(x īj + ūj ,max(0, g -1 j ( cī +g j (x īj )) )))

if ūj > 0 and i ∈ N f j (min(x īj + ūj ,max(0,g -1 j ( cī +g j (x īj ))))) g j (min(x īj + ūj ,max(0,g -1 j ( cī +g j (x īj )))))

if ūj > 0 and i / ∈ N ;

15:

end for 16:

 = arg max j∈N {r j };

17:

x ī := min(x ī + ū , max(0, g -1  ( cī + g (x ī ))));

18:

if  ∈ N then x ī := x ī ;

19:

ū := u  -i∈M x i , cī := c ī -j∈N g j (x īj ); In the first phase (Steps 6-10) the current knapsack is filled with the surrogate solution if no relaxed constraint is violated (see Step 8). Since the items are sorting according to their profit-to-weight ratios, the current knapsack is filled at first with the more promising items.

In the second phase we identify the knapsack ī with the largest residual capacity (Step 11), we update the profit-to-weight ratios (Steps 13-15), and we determine the item  with the best residual profit-to-weight ratio (Step 16). Knapsack ī is then filled as much as possible with item  (Steps 17-18). We iteratively update the upper bound and the residual capacity (Step 19), and we repeat the previous steps, as long as the largest residual capacity or the largest residual upper bound are strictly positive.

The two for loops 2-5 and 6-10 are executed mn times. The while loop is executed at most max(m, n) times. Moreover, if we assume, as in Section 4.4.1, that computing the (pseudo-)inverse of the weight function g j takes a time bounded by a constant independent from the input size, each for loop 13-15 takes n times. Hence, the overall time complexity of Procedure Surrogate-feas-1(x * ) is O(mn + n 2 ).

In short, the Procedure Surrogate-feas-1(x * ) starts with an empty solution and tries to construct a feasible solution that replicates as much as possible the (infeasible) surrogate solution. Procedure Surrogate-feas-2(x * ) starts instead with the surrogate solution and considers the items in reverse order, i.e., according to non-decreasing r j values (see (4.15)). The idea is to reduce the quantity x * ij of each item j in the current knapsack i, until the capacity constraints (4.1b) are satisfied . This is obtained by iteratively reducing the quantity of the item with the worst profit-to-weight ratio, so as to undermine as little as possible the quality of the surrogate solution.

Remark 4.4.2. The surrogate solution already satisfies the upper bound constraints (4.1c), and therefore, if we reduce the units x * ij , the new solution still meets the constraints (4.1c).

The for loop 5-8 is executed mn times. If, we assume, as above, that time to compute the (pseudo-)inverse of function g j can be bounded by constant independent from the input size, the overall time complexity of Procedure Surrogate-feas-2(x * ) is O(mn).

. . Local Search

In this section we introduce a post-processing local search procedure in order to improve the quality of the solution found by the heuristics. The local search implements pairwise exchanges of the amounts of items assigned to the same knapsack i.

For a given knapsack i, the local search considers a pair of items j and k and applies two small variations to them. Let ε be a new incremental step smaller than δ j and δ k (see Section 4.4.1). The local search simultaneously Algorithm 8 Procedure Surrogate-feas-2(x * ).

1: for i := 1 to m do for j := 1 to n do x ij := x * ij ; 2: for i := 1 to m do 3:

ci := max(0, j∈N g j (x ij ) -c i ); x ij := max(0, g -1 j (g j (x ij )ci ));

6:

if j ∈ N then x ij := x ij ;

7:

ci := max(0, j∈N g j (x ij ) -c i );

8:

if ci 0 then break;

9:
end for 10: end for increases (resp. decreases) x ij and decreases (resp. increases) x ik by ε units, respectively. Then, it computes the following variations for the objective function:

1. ∆ 1 := (f j (x ij + ε) -f j (x ij )) + (f k (x ik -ε) -f k (x ik )); 2. ∆ 2 := (f j (x ij -ε) -f j (x ij )) + (f k (x ik + ε) -f k (x ik )).
Further impose that a ∆ ( = 1, 2) takes the value 0 if the corresponding variation is infeasible, i.e., if either a right-hand side (u j , u k , or c i ) of inequalities (4.1b)-(4.1c) is exceeded or one of the two variables takes a negative value. Let ∆ = max(∆ 1 , ∆ 2 ):

• if ∆ > 0 the procedure, shown in Algorithm 9, (i) performs the corresponding variation, producing a new solution with objective function value increased by ∆;

(ii) iterates the process, for the same couple of items and ε, obviously by only computing the ∆ ( = 1 or 2) that produced ∆;

• if instead ∆ 0, i.e., both variations either worsen the solution value or are infeasible, the next couple of items is tested, or the next knapsack is considered (when j = n -1 and k = n).

Preliminary computational experiments are conducted by varying the selection procedure for the knapsack i and for the items j and k: for instance, we select the knapsack and the couple of items randomly. However, the deterministic version produces the most satisfactory results. define an appropriate value ε < min(δ j , δ k ); if ∆ > 0 then apply the variation corresponding to ∆; The inner repeat-until loop is executed mn 2 times. This loop can theoretically take a pseudo-polynomial time but, in practice, it is executed for a number of times bounded by a constant independent from the input size. Over the 3,360 instances we tested the algorithm (see Section 4.4.5), the number of iterations was normally between 1 and 2, and it never attained 10. The overall computational time of the local search is O(mn 2 ) time.

. . Overall Algorithm

The overall heuristic algorithm for the MNLKP can be stated as follows:

Algorithm 10 DMM.

. . Computational Experiments

In order to evaluate the computational behavior of the algorithms of Section 4.4.4, we compared it with open-source local solvers for nonlinear programming (Ipopt [126] for real instances and Bonmin [START_REF] Bonmin | [END_REF], with option bonmin. algorithm B-BB, for integer instances) and with a global solver (Couenne [59], both for real and integer instances). Ipopt and Bonmin are exact solvers for convex NLPs and convex MINLPs, respectively: hence, they can be used to min(δ j , δ k )/2. All computations of the zero of a weight function needed by DMM were performed through a binary search over the definition range. (The impact on the overall CPU time was however negligible.) report the results for nonlinear weights, with real and integer variables x ij for similar capacities. The entries are, for real variables: the average values produced by DMM, Ipopt with a single starting point (Ipopt_1), Ipopt with 10 starting points (Ipopt_10), and Couenne, and for integer variables: the average values produced by DMM, Bonmin with a single starting point (Bonmin_1), Bonmin with 10 starting points at the root node (Bonmin_-10), and Couenne. For the group of instances, for which the solvers do not succeed in finding a feasible solution within the time limit, the tables report also in brackets the number of non-solved instances.

Tables 13-14 report the same statistic for dissimilar capacities. The tables with odd numbering present the average solution values produced over the 20 generated instances, while those with even numbering report the corresponding average CPU times (in seconds).

Tables 15-18 report the same information for the case of linear weights for similar and dissimilar capacities.

Tables 11 and13 clearly show that, on the nonlinear instances, a part for few smaller instances with n = 10 (and a single case for n = 20), the proposed algorithm almost always outperforms both the exact and the heuristic solvers. Tables 12 and14 show that on the integer instances DMM is always the fastest method with regard to the elapsed CPU time. For the real instances Ipopt_1 is generally faster, but the solution values it produces are definitely worse (by over 10% on average). Overall, DMM seems to be a reasonable algorithm with respect to the trade-off between the quality of the produced solution value and the computational time.

For the instances with linear weight functions, Tables 15 and17 show that Ipopt_10 and Bonmin_10 often provide better solutions for smaller instances, while, instead, DMM always performs better for n 200. Concerning the average CPU times (Tables 16 and18), DMM is again the clear winner on the integer instances. For the real instances, Ipopt_1 is faster for m = 2 (but it produces on the other hand worse solution values), while DMM always outperforms the open-source solvers for m 5 (the difference is particularly high on the larger instances).

It turns out that the instances with linear weights are more difficult to solve to optimality than those with nonlinear weights. Although this can appear surprising, there is no theoretical result implying that one case must be easier than the other. Couenne, for example, transforms the objective function so as it becomes linear, while its nonlinear terms become additional For all the real instances, all the algorithms (DMM, Ipopt_1, and Ipopt_-10) are able to provide a feasible solution within the time limit (both for the case of linear and nonlinear weights); while, for the larger integer instances, DMM is the only heuristic with can determine a feasible solution within the time limit.

Tables 19 and20 compare the solution values produced by DMM and Couenne, for those instances, with respect to which Couenne is able to find the global maximum. Table 19 reports the results for nonlinear weights and both similar and dissimilar capacities, while Table 20 is the analogue for linear weights. DMM computes solutions with an average error of 3.65%, with a minimum of 0% and a maximum of 18% with regards to the global optimum. We observe that Couenne optimally solved more instances with nonlinear weights (5.6%) than with linear weights (1.8%).

Finally, for 10 of the instances globally solved by Couenne, it turned out that the software behaved incorrectly, finding a solution value lower that the one produced by DMM. Hence, at the suggestion of Couenne developer, we re-execute the solver by disabling several default options, namely: aggressive _ fbbt, optimality _ bt, and redcost _ bt. For instance, the solution value of the instance 6 with m = 5, n = 10, dissimilar capacities, integer variables and nonlinear weights, produced by Couenne with the default options was 286.86, while DMM provide a feasible solution of value 294.76. We re-run Couenne without the suggested options and the solution value was 305.67.

. In this chapter we considered the MNLKP and we discussed the possible relaxations for this problem. We implemented the MWU algorithm, adapting its main steps to this class of problems. However, computational experiments illustrated that this methodology was not successful in this situation; hence, we proposed different constructive heuristic strategies followed by a local search post-processing. Extensive computational tests clearly showed that this approach outperformed the other heuristic and exact algorithms available for the MNLKP both in terms of solution quality and of CPU time elapsed.

Conclusions

C O N C L U S I O N S

In this thesis we have described a new heuristic algorithm for Mixed Integer NonLinear Problems, based on the general Multiplicative Weights Update algorithm. The contribution of the thesis is twofold: from one side we have theoretically introduced the algorithm framework with a new class of mathematical reformulation, namely the pointwise reformulation, whose theoretical properties are inspected in detail; from the other side we have illustrated two real-world optimization problems, namely the Mean-Variance Portfolio Selection and the Multiple NonLinear Knapsack Problems, with regard to which we practically applied and implemented the algorithm.

In the theoretical part, we analyzed in depth the steps of the algorithm and the general properties of the pointwise reformulation with a particular emphasis on exactness and efficiency. A pointwise reformulation is exact if there exists a parameter such that the global optimality properties of the original formulation are preserved; while a reformulation is efficient if it can be solved in polynomial time. In case the reformulation is exact and efficient, it is sufficient to solve the (easier) reformulation in order to obtain a global solution for the original problem. The Multiplicative Weights Update for Mixed Integer NonLinear Problems is, in fact, a MultiStart type algorithm with a specific choice of the starting points, conducted according to the Multiplicative Weights Update algorithm. Moreover, we described several general building strategies for the pointwise reformulation with respect to given classes of mathematical optimization problems, such as Polynomial, Bilinear, and Quadratic Mixed Integer NonLinear Problems.

In the second part, we examined two different real-world problems and we adapted the general framework to these programs. First, we gave a sufficiently complete survey of the Mean-Variance Portfolio Selection Problems about the ways to model the behaviors of investors and the restrictions of the financial markets, robust and stochastic approaches, exact reformulations and inner approximations, and exact algorithms to solve these challenging problems. Then, we introduced a Multiplicative Weights Update for a general class of portfolio optimization problems, specifying how to generate the pointwise reformulation and how to choose the promising starting points for the original formulation. The second problem we dealt with is the Multiple NonLinear (Separable) Knapsack Problem, characterized by a set of items and a set of knapsacks and whose target consists in filling the knapsacks with units of items in order to maximize a given profit function and meeting knapsack and upper bound constraints. In the general scheme, several variables can be required to be integer. We apply the Multiple Weights Update 91 to this kind of problems, but the computational behavior of this algorithm was quite poor in terms of performances: in facts, it was almost always overcome by the MultiStart procedure.

Therefore, we considered a constructive heuristic based on three elements: (i) a greedy type heuristic, already proposed for the single knapsack problem, and here extended to the multiple case; (ii) a feasibility recovery strategy which, starting from the empty solution, tries to fill the knapsacks as much as possible with the solution of the surrogate relaxation; (iii) another feasibility recovery procedure which, starting from the surrogate solution, progressively removes the exceeding capacities. Then, we took the best solution produced by procedures (i)-(iii) and ran a local search postprocessing. The computational tests indicate the resulting algorithm outperformed heuristic and exact solution methods available for these knapsack problems.

We hope that the Multiplicative Weights Update algorithmic methodology will be applied with success also to other Mixed Integer NonLinear Problems. However, several questions are still open. It could be interesting to consider an analysis of the goodness of a given pointwise reformulation a priori, i.e., a way to theoretically evaluate several formulations and choose between them. We felt that the general approach presented for the portfolio and the knapsack problems could be easily exploited for other problems characterized by separable non-convexities and non-concavities. Moreover, it could be an engaging target defining other automatic strategies to build the pointwise reformulation for other classes of problems and testing them over general Mixed Integer NonLinear Problems Libraries. Therefore, there could be also the possibility to employ several different reformulations at the same time, alternating from one iteration to another the set with respect to which the parameter θ is defined: this could be the case of bilinear optimization problems where the decision variables can be divided into two different disjunctive sets. Finally, we point out that the definition of the pointwise reformulation could be exploited also in more general situations than the Multiplicative Weights Update Algorithm: it generalizes the definition of relaxations to bounding reformulations since by solving them we always obtain a lower bound for the original optimization problem.

Figure

  Figure Examples of transaction cost functions. 52 Figure MVPS, CPU time vs. size of the problem n (#assets). 58 Figure Example of profit function. 61

.

  The pseudocode of the MWU algorithm for MINLPs is shown in Algorithm 4. It takes a MINLP formulation P as input and produces a local solution as output. In the next sections we analyze in more detail the steps of the algorithm.Algorithm 4 MWU Algorithm for MINLPs1: assign weights w 1 := 1, incumbent (x * , y * ) := (∞, ∞) and t := 1 2: while t T do 3:

Theorem 3 . 3 . 7 .

 337 (Bonami and Lejeune [26]) Let p ∈ [0.5, 1). If the skewness of the probability distribution of ξ x is positive, then the deterministic constraint (3.7b), equivalent to (3.5), is a SOCC.

Theorem 3 . 6 . 5 .

 365 (Filomena and Lejeune [89]) Program (3.36) is equivalent to the following NLP: min x,h h s.t. (3.36b), (3.36c), (3.36d), (3.36e).

  38): min x,h,b h s.t. (3.38b), (3.38c), (3.38d), (3.38e), (3.38f).

  1c), (3.1b), (3.1d), (3.10),(3.16)

  Hard trigonometric transaction costs.

Figure 1 :

 1 Figure 1: Examples of transaction cost functions.

  "Easy" concave transaction costs.

  "Hard" concave transaction costs.

Figure 2 :

 2 Figure 2: MVPS, CPU time vs. size of the problem n (# assets).

Figure 3 :

 3 Figure 3: Example of profit function.

  1, . . . , m -1).

11 :

 11 if ci > 0 then {comment: fill the residual capacity with item n} 12:

1 : 4 : 7 : 8 :

 1478 for j := 1 to n do ūj := u j ; 2: for i := 1 to m do 3: ci := c i ; for j := 1 to n do x ij := 0; 5: end for 6: for i := 1 to m do for j := 1 to n do if g j (x * ij ) ci then x ij := x * ij , ci := ci -g j (x * ij ), ūj := ūj -x * ij ; end for 11: ī = arg max i∈M { ci }, cmax := cī , ūmax := max j∈N { ūj }; 12: while cmax > 0 and ūmax > 0 do 13:

20 :ī

 20 = arg max i∈M { ci }, cmax := cī , ūmax := max j∈N { ūj }; 21: end while

4 :

 4 for j := n back to 1 do 5:

Algorithm 9 2 : 3 :

 923 Procedure Local Search. 1: for i := 1 to m do for j := 1 to n -1 do for k := j + 1 to n do 4:

compute ∆ 1 ,

 1 ∆ 2 , and ∆; 7:

7 :

 7 return the solution corresponding to Z l .

  follows that, in any case, the feasible region is defined by nonlinear constraints, and the evolution of the branching search is unpredictable.

Table 20 :

 20 MNLKP, linear weights. Solution values for instances globally solved by Couenne.

  Titre : L'Algorithme Multiplicative Weights Update pour la Programmation non linéaire en nombres entiers: Théorie, Applications et Limites Mots clefs : Algorithme Multiplicative Weights Update, Programmation non linéaire en nombres entiers, Sélection du portefeuille moyenne-variance, Sac à dos multiple non linéaire Résumé : L'objectif de cette thèse consiste à présenter un nouvel algorithme pour la programmation non linéaire en nombres entiers, inspirée par la méthode Multiplicative Weights Update et qui compte sur une nouvelle classe de reformulations, appelées les reformulations ponctuelles. La programmation non linéaire en nombres entiers est un sujet très difficile et fascinant dans le domaine de l'optimisation mathématique à la fois d'un point de vue théorique et computationnel. Il est possible de formuler de nombreux problèmes dans ce schéma général et, habituellement, ils posent de réels défis en termes d'efficacité et de précision de la solution obtenue quant aux procédures de résolution. La thèse est divisée en trois parties principales: une introduction composée par le Chapitre 1, une définition théorique du nouvel algorithme dans le Chapitre 2 et l'application de cette nouvelle méthodologie à deux problèmes concrets d'optimisation, tels que la sélection optimale du portefeuille avec le critère moyenne-variance dans le Chapitre 3 et le problème du sac à dos non linéaire dans le Chapitre 4. Conclusions et questions ouvertes sont présentées dans le Chapitre 5. Title : The Multiplicative Weights Update Algorithm for Mixed Integer NonLinear Programming: Theory, Applications, and Limitations Keywords : Mixed Integer NonLinear Programming, Multiplicative Weights Update Algorithm, Mean-Variance Portfolio, Multiple NonLinear Knapsack Problem Abstract : This thesis presents a new algorithm for Mixed Integer NonLinear Programming, inspired by the Multiplicative Weights Update framework and relying on a new class of reformulations, called the pointwise reformulations. Mixed Integer NonLinear Programming is a hard and fascinating topic in Mathematical Optimization both from a theoretical and a computational viewpoint. Many real-world problems can be cast this general scheme and, usually, are quite challenging in terms of efficiency and solution accuracy with respect to the solving procedures. The thesis is divided in three main parts: a foreword consisting in Chapter 1, a theoretical foundation of the new algorithm in Chapter 2, and the application of this new methodology to two real-world optimization problems, namely the Mean-Variance Portfolio Selection in Chapter 3, and the Multiple NonLinear Separable Knapsack Problem in Chapter 4. Conclusions and open questions are drawn in Chapter 5. Université Paris-Saclay Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

Table 13 MNLKP

 13

	54

47 Table MVPS, comparative results of MS and MWU for the transaction cost function (a). , nonlinear weights, dissimilar capacities. Average solution values over 20 instances (#no solution).

Table MNLKP ,

 MNLKP nonlinear weights, dissimilar capacities. Average CPU times over 20 instances (#no solution). 80 Table MNLKP, linear weights, similar capacities. Average solution values over 20 instances (#no solution).

			82
	Table	MNLKP, linear weights, similar capacities. Average
		CPU times over 20 instances (#no solution).	83
	Table	MNLKP, linear weights, dissimilar capacities. Aver-
		age solution values over 20 instances (#no solution).	84
	Table	MNLKP, linear weights, dissimilar capacities. Aver-
		age CPU times over 20 instances (#no solution).	85
	Table	MNLKP, nonlinear weights. Solution values for in-
		stances globally solved by Couenne.	86
	Table	MNLKP, linear weights. Solution values for instances
		globally solved by Couenne.

  Proposition 1.2.4. Let Assumption 1.2.2 hold and suppose the problem is feasible. Let {(x * t , y * t )} be a sequence of points generated by the Algorithm 1. Then, for any ε > 0, we have

  3 we illustrate several examples of automatic procedures to build the pointwise reformulation for specific classes of 17 MINLPs. Finally, in Section 2.4 a complete description of the MWU algorithm for MINLPs is specified.

	Definition 2.2.1. Given a MINLP P as in (1.1), a pointwise reformulation R θ =
	ptw	(P) is a family of MINLP formulations, depending on a parameter θ =
	t←t'(θ)
	(θ s | s r), which are obtained by replacing terms t 1 , . . . , t r in P by corresponding
	parametrized terms t s (θ s ) (for s r).

.

Definition 2.2.4. Given

  a MINLP P and its pointwise reformulation R θ : (a) R θ is spanning if, for any x ∈ R n and y ∈ Z p , there are values of θ such that evaluating the functions of P and of R θ at (x, y) determines the same value, i.e., such that

			.1d)
	Remark 2.2.3. The number of the constraints in the pointwise reformulation (2.1)
	is the same as in the original formulation (1.1); however, since in the pointwise
	reformulation (2.1) several terms are substituted by means of parameters θ, the
	number of variables could be different from the one of the original formulation (1.1).
	We introduce the following classification for the pointwise reformulations.
	∀s r	D s ⊂	θ s ∈Θ s

  1 , . . . , h t ) and d = (d 1 , . . . , d t ), let

		µ d h =	d j x h j
		j t
	be a monomial occurring in P. Let us see the MILP pointwise reformulation
	first:
	1. introduce new parameters θ 1 , . . . , θ t ;
	2. for each monomial µ d h in P:
	a) choose	t as = argmin j t {d j }
	b) replace µ by the linear term a µ x h , where

Remark 2.3.4. Let

  λ min be the minimum eigenvalue of Q , and 1 n×n indicate the identity matrix of order n. B = λ min 1 n×n is a good solution for both (2.8) and (2.9) (see Section 4 in[START_REF] Frangioni | Perspective Cuts for a Class of Convex 0-1 Mixed Integer Programs[END_REF]).

	Remark 2.3.5. By construction, the pointwise convex MIQP reformulation previ-
	ously introduced is spanning for all θ = x.
			.8c)
	A convex relaxation of (2.8) can be obtained by replacing, in the objective
	function, the 0 -norm with 1 -norm, namely
	min B 1	(2.9a)
	s.t. A	0	(2.9b)
	A + B = Q .	(2.9c)

Table 2 :

 2 MVPS, comparative results of MS and MWU for the transaction cost function (a).

Table 3 :

 3 MVPS, comparative results of MS and MWU for the transaction cost function (b).

Table 4 :

 4 MVPS, comparative results of MS and MWU for the transaction cost function (c).

	3.231
	-35.917%
	3.034

Table 5 :

 5 MVPS, comparative results of MS and MWU for the transaction cost function (d).

	4.306
	-77.499%
	4.887

Table 6 :

 6 MVPS, comparative results of MS and MWU for the transaction cost function (e).

	geo (CPU)	std	avg		port-orlib5 _ 150	port-orlib5 _ 100	port-orlib5 _ 050	port-orlib5 _ 000	port-orlib4 _ 150	port-orlib4 _ 100	port-orlib4 _ 050	port-orlib4 _ 000	port-orlib3 _ 150	port-orlib3 _ 100	port-orlib3 _ 050	port-orlib3 _ 000	port-orlib2 _ 150	port-orlib2 _ 100	port-orlib2 _ 050	port-orlib2 _ 000	port-orlib1 _ 150	port-orlib1 _ 100
					0.0003272876	0.0003918260	0.0005150277	0.0016485224	0.0001613979	0.0003059553	0.0006828450	0.0029387241	0.0002239117	0.0003215941	0.0005849758	0.0015166351	0.0001663200	0.0002704062	0.0004953237	0.0028352430	0.0007158421	0.0010585969
					225	225	225	225	98	98	98	98	89	89	89	89	85	85	85	85	31	31
		(0.0145719)	0.0269742	0.5394835	0.0143170	0.0078564	0.0031309	0.0326302	0.0236738	0.0135301	0.0583659	0.0130890	0.0284570	0.0342867	0.0347923	0.0116564	0.0315010	0.0283756	0.0383876	0.0136601	0.0399172	0.0432059
	135.816	(4819.775)	2489.849	49796.978	968.684	5516.043	302.041	1448.934	15929.002	1082.307	32.832	4.641	15650.882	903.586	19.026	9.548	3642.877	4226.211	39.267	10.510	2.588	2.984
		(2.685)	7.950	159	10	7	9	6	10	10	10	2	10	9	8	2	10	10	7	4	9	10
		(0.0155312)	0.0264158	0.5283151	0.0143169	0.0035668	0.0073133	0.032836	0.0260514	0.0078843	0.0541541	0.0123611	0.0269692	0.0392775	0.0348732	0.0029338	0.0307374	0.0166890	0.0387789	0.0133721	0.0376060	0.0413710
	247.786	(5846.227)	2705.947	54118.943	730.436	4774.516	370.425	788.687	24604.167	696.938	91.192	66.034	11231.379	2142.790	102.802	22.035	3989.962	4382.785	33.822	8.726	6.172	4.962
		(1.917)	8.900	178	10	8	10	6	10	10	10	10	10	10	8	8	10	10	6	3	9	10
		(0.0135588)	0.0244296	0.4885919	0.0143170	0.0103946	0.0018351	0.0017349	0.0252765	0.0248876	0.0239543	0.0233697	0.0255798	0.0188460	0.0239912	0.0234077	0.0268662	0.0207083	0.0186861	0.0184047	0.0388779	0.0498939
	426.464	(1891.435)	1256.280	25125.592	913.920	1525.448	350.660	5513.344	5982.823	513.718	272.405	432.908	5050.864	719.657	441.404	240.592	1599.415	524.517	522.376	375.738	6.757	58.515
		(0.308)	9.900	198	10	10	9	10	10	10	10	10	10	10	10	10	10	10	10	10	9	10
		(56.991%)	-5.209%	-104.177%	0.000%	-32.308%	41.386%	94.683%	-6.770%	-83.943%	58.958%	-78.545%	10.111%	45.034%	31.044%	-100.814%	14.713%	27.021%	51.323%	-34.734%	2.604%	-15.479%
		(21.729)	13.244	264.875	0.943	0.277	1.161	3.805	0.376	0.475	8.297	93.279	0.323	0.796	23.200	25.198	0.439	0.124	13.303	35.751	2.611	19.610	11.903
		(172.079%)	-46.210%	-924.200%	0.000%	-191.430%	74.907%	94.717%	2.974%	-215.660%	55.766%	-89.058%	5.152%	52.018%	31.204%	-697.857%	12.594%	-24.083%	51.814%	-37.635%	-3.382%	-20.601%	-2.415%
		9.900	5.575	111.494	1.251	0.319	0.947	6.991	0.243	0.737	2.987	6.556	0.450	0.336	4.294	10.919	0.401	0.120	15.445	43.060	1.095	11.793	0.665

Table 7 :

 7 MNLKP, nonlinear weights, similar capacities. Average solution values over 20 instances.

				Real Variables
	m	n	MWU Ipopt_1 Ipopt_10	Couenne
	2	10	326.57	313.50	351.70	362.63
	2	20	622.24	569.20	635.35 593.43(12)
	2	50 1,777.83 1,714.74	1,799.78	n/a(20)
	2 total 2,726.64 2,597.44	2,786.83	-
	5	10	266.09	271.61	312.63	299.06(5)
	5	20	677.78	687.95	733.25	n/a(20)
	5	50 1,705.66 1,705.86	1,782.99	n/a(20)
	5 total 2,649.53 2,655.42	2,828.87	-
	10	10	195.08	187.40	218.59	212.07(2)
	10	20	658.91	664.58	702.21	n/a(20)
	10	50 1,873.94 1,864.39	1,927.13	n/a(20)
	10 total 2,727.93 2,716.37	2,847.93	-
	total total 8,104.10 7,979.23	8,463.63	-

Table 8 :

 8 MNLKP, nonlinear weights, similar capacities. Average CPU times over 20 instances.

				Real Variables
	m	n MWU Ipopt_1 Ipopt_10	Couenne
	2	10	1.48	0.05	0.51	1,113.88
	2	20	8.08	0.10	1.18 3,601.81(12)
	2	50	33.01	0.34	3.45	n/a(20)
	2 total	42.57	0.49	5.14	-
	5	10	8.75	0.16	1.72	3,600.61(5)
	5	20	22.82	0.32	3.57	n/a(20)
	5	50	68.97	1.17	11.76	n/a(20)
	5 total 100.54	1.65	17.05	-
	10	10	15.88	0.39	3.78	3,600.54(2)
	10	20	43.02	0.93	8.57	n/a(20)
	10	50 122.74	3.35	30.70	n/a(20)
	10 total 181.64	4.67	43.05	-
	total total 324.75	6.81	65.24	-

Table 9 :

 9 MNLKP, nonlinear weights, dissimilar capacities. Average solution values over 20 instances.

				Real Variables
	m	n	MWU Ipopt_1 Ipopt_10	Couenne
	2	10	311.18	298.24	339.50	355.19
	2	20	600.98	547.68	610.20	564.65(9)
	2	50 1,766.38 1,619.80	1,737.00	n/a(20)
	2 total 2,678.54 2,465.72	2,686.70	-
	5	10	288.27	285.28	318.68	329.51(1)
	5	20	688.35	670.57	717.91 659.84(17)
	5	50 1,654.93 1,614.39	1,723.24	n/a(20)
	5 total 2,631.55 2,570.24	2,759.83	-
	10	10	316.10	314.88	345.18	354.52
	10	20	687.13	678.59	741.04 699.31(11)
	10	50 1,825.00 1,807.00	1,909.65	n/a(20)
	10 total 2,828.23 2,800.47	2,995.87	-
	total total 8,138.32 7,836.43	8,442.40	-

Table 10 :

 10 MNLKP, nonlinear weights, dissimilar capacities. Average CPU times over 20 instances.

				Real Variables
	m	n MWU Ipopt_1 Ipopt_10	Couenne
	2	10	1.84	0.05	0.53	794.39
	2	20	8.59	0.11	1.14	3,307.95(9)
	2	50	33.88	0.34	3.36	n/a(20)
	2 total	44.31	0.50	5.03	-
	5	10	9.27	0.14	1.53	2,790.28(1)
	5	20	24.62	0.34	3.23 3,600.18(17)
	5	50	68.75	1.16	10.25	n/a(20)
	5 total 102.64	1.64	15.01	-
	10	10	16.06	0.31	2.62	2,266.36
	10	20	42.70	0.79	6.10 3,599.61(11)
	10	50 126.30	2.89	23.01	n/a(20)
	10 total 185.06	3.99	31.73	-
	total total 332.01	6.13	51.77	-

  OtherLagrangian relaxations can be obtained by multiplying the upper bound constraints (4.1c) by a vector of non negative multipliers or relaxing both the constraints (4.1b) and (4.1c). Preliminary computational experiments have shown that the solutions produced by the Lagrangian relaxations are generally worse than the ones produced by the surrogate relaxation.

.13d) Remark 4.3.2. The Lagrangian relaxation (4.13) can be decomposed into n independent subproblems, one for each item j, with nonlinear objective functions.

  : execute procedures Surrogate-feas-1(x * ) and Surrogate-feas-2(x * ) of Section 4.4.2 and let Z s be the best solution value; 5: if Z s = U then terminate else Z := max(Z h , Z s ); 6: execute Local Search of Section 4.4.3 on the solution corresponding to Z and let Z l be the solution value;

	1: solve the surrogate relaxation (4.11) and let U and x * be respectively the
	resulting upper bound and the corresponding solution;
	2: execute procedure Constructive of Section 4.4.1 and let Z h be the solu-
	tion value;

3: if Z h = U then terminate; 4

Table 11 :

 11 MNLKP, nonlinear weights, similar capacities. Average solution values over 20 instances (# no solution).

	Integer Variables	DMM Bonmin_1 Bonmin_10 Couenne	350.47 327.21 313.87 364.45	645.63 591.26 589.88 561.83	1,902.65 1,769.76 1,800.63 1,481.36	3,742.14 3,413.98 3,454.54 -	7,168.07 6,429.80(5) 6,419.05 -	18,376.30 16,286.80(13) 16,404.90(5) -	37,058.30 n/a(20) n/a(20) -	69,243.56 ---	321.18 288.19 281.05 299.35	727.24 706.88 716.23 620.54(2)	1,821.57 1,753.76(2) 1,765.72 1,245.52(1)	3,868.62 3,719.86(1) 3,825.23 -	7,850.18 n/a(20) 7,677.54(16) -	19,273.73 n/a(20) n/a(20) -	38,541.14 n/a(20) n/a(20) -	72,403.66 ---	232.58 183.37(2) 201.60(1) 222.29	727.85 719.17(3) 681.76(2) 489.35(5)	1,983.72 1,834.39(16) 1,937.89 1,275.43(16)	3,957.84 n/a(20) 3,852.28(19) -	7,655.04 n/a(20) n/a(20) -	19,645.34 n/a(20) n/a(20) -	39,721.93 n/a(20) n/a(20) -	73,924.30 ---	215,571.52 ---
	Real Variables	m n DMM Ipopt_1 Ipopt_10 Couenne	2 10 350.49 313.50 351.70 362.63	2 20 641.73 569.20 635.35 593.43(12)	2 50 1,900.41 1,714.74 1,799.78 n/a(20)	2 100 3,741.92 3,341.04 3,548.40 -	2 200 7,167.80 6,302.04 6,683.39 -	2 500 18,375.35 16,153.10 16,942.30 -	2 1000 37,051.16 32,117.80 33,899.30 -	2 total 69,228.86 60,511.42 63,860.22 -	5 10 322.76 271.61 312.63 299.06(5)	5 20 729.46 687.95 733.25 n/a(20)	5 50 1,822.58 1,705.86 1,782.99 n/a(20)	5 100 3,865.64 3,701.10 3,818.14 -	5 200 7,846.03 7,372.56 7,681.43 -	5 500 19,272.60 18,180.80 18,825.70 -	5 1000 38,540.16 36,391.00 37,485.80 -	5 total 72,399.23 68,310.88 70,639.94 -	10 10 216.19 187.40 218.59 212.07(2)	10 20 734.75 664.58 702.21 n/a(20)	10 50 1,983.10 1,864.39 1,927.13 n/a(20)	10 100 3,952.78 3,732.51 3,843.88 -	10 200 7,652.05 7,311.37 7,431.36 -	10 500 19,640.75 18,788.40 19,097.90 -	10 1000 39,717.69 36,316.50 38,266.70 -	10 total 73,897.31 68,865.15 71,487.77 -	total total 215,525.40 197,687.45 205,987.93 -

Table 12 :

 12 MNLKP, nonlinear weights, similar capacities. Average CPU times over 20 instances (# no solution).

	2	m
	10 1.08 0.05 0.51	n DMM Ipopt_1 Ipopt_10	Real Variables
	1,113.88	Couenne
	1.09	DMM
	3.71 1.47	Bonmin_1 Bonmin_10	Integer Variables
	847.17	Couenne

Table 13 :

 13 MNLKP, nonlinear weights, dissimilar capacities. Average solution values over 20 instances (# no solution).

	Integer Variables	DMM Bonmin_1 Bonmin_10 Couenne	339.89 299.21 308.60 354.62	636.54 569.67 573.25 555.90	1,871.45 1,688.48 1,653.26 1,362.62	3,722.43 3,340.35(1) 3,294.41 -	7,072.53 6,188.82(3) 5,905.41 -	18,287.28 15,883.50(14) 15,862.40(5) -	37,086.84 n/a(20) n/a(20) -	69,016.96 ---	313.70 305.05(1) 299.28(1) 319.22(1)	720.97 692.01(2) 664.84 598.87(5)	1,772.33 1,600.30(4) 1,627.85 1,228.18(4)	3,805.91 3,406.69(2) 3,400.23 -	7,595.80 n/a(20) n/a(20) -	19,100.00 n/a(20) n/a(20) -	38,472.11 n/a(20) n/a(20) -	71,780.82 ---	321.47 315.20(1) 307.11 340.21(1)	737.34 703.57(3) 681.37(1) 633.95(8)	1,883.62 1,817.64(9) 1,810.51(1) 1,499.67(7)	3,615.12 n/a(20) 3,804.82(19) -	7,395.62 n/a(20) n/a(20) -	19,242.96 n/a(20) n/a(20) -	39,403.01 n/a(20) n/a(20) -	72,599.14 ---	213,396.92 ---
	Real Variables	m n DMM Ipopt_1 Ipopt_10 Couenne	2 10 340.62 298.24 339.50 355.19	2 20 634.68 547.68 610.20 564.65(9)	2 50 1,870.83 1,619.80 1,737.00 n/a(20)	2 100 3,722.84 3,236.18 3,437.91 -	2 200 7,074.48 5,908.90 6,313.01 -	2 500 18,280.25 15,286.70 16,081.20 -	2 1000 37,089.65 31,168.30 32,825.20 -	2 total 69,013.35 58,065.80 61,344.02 -	5 10 321.28 285.28 318.68 329.51(1)	5 20 728.42 670.57 717.91 659.84(17)	5 50 1,792.76 1,614.39 1,723.24 n/a(20)	5 100 3,810.80 3,392.78 3,625.14 -	5 200 7,647.44 6,778.96 7,180.81 -	5 500 19,108.46 17,131.00 17,967.30 -	5 1000 38,447.52 34,559.20 36,099.60 -	5 total 71,856.68 64,432.18 67,632.68 -	10 10 343.95 314.88 345.18 354.52	10 20 757.28 678.59 741.04 699.31(11)	10 50 1,972.11 1,807.00 1,909.65 n/a(20)	10 100 3,873.07 3,441.87 3,684.99 -	10 200 7,543.03 6,708.00 7,094.40 -	10 500 19,271.73 17,027.60 17,950.80 -	10 1000 39,486.36 35,248.30 37,237.50 -	10 total 73,247.53 65,226.24 68,963.56 -	total total 214,117.56 187,724.22 197,940.26 -

Table 14 :

 14 MNLKP, nonlinear weights, dissimilar capacities. Average CPU times over 20 instances (# no solution).

	5 50 5.16 1.16 10.25 n/a(20)	5 20 2.03 0.34 3.23 3,600.18(17)	5 10 1.20 0.14 1.53 2,790.28(1)	2 total 210.01 48.14 465.36 -	2 1000 113.39 33.71 318.10 -	2 500 57.34 10.64 105.76 -	2 200 20.87 2.40 26.51 -	2 100 10.23 0.89 9.96 -	2 50 5.08 0.34 3.36 n/a(20)	2 20 2.02 0.11 1.14 3,307.95(9)	2 10 1.08 0.05 0.53 794.39	m n DMM Ipopt_1 Ipopt_10 Couenne	Real Variables
	5.17 2,716.96(4) 773.89 3,604.77(4)	2.05 250.53(2) 588.38 3,600.92(5)	1.23 29.11(1) 196.66(1) 2,013.34(1)	209.39 ---	115.93 n/a(20) n/a(20) -	54.07 3,265.7(14) 3,456.18(5) -	20.71 3,254.23(3) 3,099.98 -	10.26 1,731.33(1) 813.99 -	5.30 277.57 69.66 3,602.78	2.03 20.41 8.18 3,430.29	1.09 2.91 1.11 364.44	DMM Bonmin_1 Bonmin_10 Couenne	Integer Variables

Table 15 :

 15 MNLKP, linear weights, similar capacities. Average solution values over 20 instances (# no solution).

	total total 232,056.46 215,598.37 222,279.86 -	10 total 79,618.76 71,913.64 75,227.95 -	10 1000 42,766.36 36,853.00 39,276.10 -	10 500 21,309.28 20,091.20 20,535.40 -	10 200 8,376.55 8,063.39 8,224.59 -	10 100 4,040.29 3,860.72 4,008.91 -	10 50 2,095.67 2,042.56 2,122.72 1,572.96(18)	10 20 790.00 778.73 816.61 n/a(20)	10 10 240.61 224.04 243.62 211.05(15)	5 total 77,856.85 74,359.53 75,950.48 -	5 1000 41,748.25 39,406.90 40,231.00 -	5 500 20,883.19 20,123.60 20,424.60 -	5 200 8,135.86 7,927.23 8,125.15 -	5 100 3,942.02 3,808.03 3,955.77 -	5 50 2,010.44 1,990.01 2,038.93 n/a(20)	5 20 762.51 739.22 782.22 n/a(20)	5 10 374.58 364.54 392.81 243.95(1)	2 total 74,580.85 69,325.20 71,101.43 -	2 1000 39,562.06 36,752.90 37,448.70 -	2 500 20,111.10 18,685.00 19,222.20 -	2 200 7,966.74 7,401.99 7,645.14 -	2 100 3,873.34 3,585.02 3,727.55 -	2 50 1,943.44 1,841.93 1,921.17 n/a(20)	2 20 780.54 739.70 786.05 n/a(20)	2 10 343.63 318.66 350.62 343.65(6)	m n DMM Ipopt_1 Ipopt_10 Couenne	Real Variables
	232,115.98 ---	79,650.43 ---	42,782.19 n/a(20) n/a(20) -	21,321.39 n/a(20) n/a(20) -	8,382.21 n/a(20) n/a(20) -	4,043.65 n/a(20) n/a(20) -	2,097.90 2,249.09(19) 2,249.09(19) 486.79(6)	787.61 804.24(1) 807.22 329.71(4)	235.48 229.52 238.99 170.29(4)	77,890.48 ---	41,755.37 n/a(20) n/a(20) -	20,900.04 n/a(20) n/a(20) -	8,141.77 n/a(20) n/a(20) -	3,942.90 3,978.83(15) 3,900.36(14) -	2,013.13 1,991.96 1,997.33 1,366.98(11)	761.99 750.92 768.77 485.02(4)	375.28 366.75 386.39 286.99(2)	74,575.07 ---	39,563.97 n/a(20) n/a(20) -	20,107.13 18,438.10(18) 18,438.10(18) -	7,966.32 7,510.38(4) 7,574.34(1) -	3,872.29 3,681.92 3,740.08 -	1,941.38 1,896.66 1,939.18 1,494.77(10)	780.03 752.42 784.46 676.70	343.95 332.91 351.35 335.53	DMM Bonmin_1 Bonmin_10 Couenne	Integer Variables

Table 16 :

 16 MNLKP, linear weights, similar capacities. Average CPU times over 20 instances (# no solution).

	Integer Variables	DMM Bonmin_1 Bonmin_10 Couenne	1.07 3.67 2.91 3,592.36	2.02 16.76 12.94 3,602.05	5.08 145.37 195.66 3,602.23(10)	10.22 1,499.39 1,312.30 -	20.81 2,978.50(4) 2,666.15(1) -	58.50 1,4671.80(18) 11,718.9(18) -	118.58 n/a(20) n/a(20) -	216.28 ---	1.13 350.92 330.21 3,601.28(2)	2.04 550.75 675.76 3,601.41(4)	5.14 3,110.13 3,182.53 3,608.20(11)	10.32 3,585.95(15) 3,292.52(14) -	21.18 n/a(20) n/a(20) -	57.23 n/a(20) n/a(20) -	128.10 n/a(20) n/a(20) -	225.14 ---	1.27 1,088.11 934.47 3,601.56(4)	2.06 3,130.17(1) 3,097.31 3,602.08(4)	5.18 3,583.36(19) 3,417.98(19) 3,610.75(6)	10.57 n/a(20) n/a(20) -	22.18 n/a(20) n/a(20) -	63.64 n/a(20) n/a(20) -	163.09 n/a(20) n/a(20) -	267.99 ---	709.41 ---
	Real Variables	m n DMM Ipopt_1 Ipopt_10 Couenne	2 10 1.06 0.04 0.56 2,250.10(6)	2 20 2.02 0.13 1.40 n/a(20)	2 50 5.06 0.50 6.02 n/a(20)	2 100 10.21 1.69 17.64 -	2 200 20.71 5.90 58.31 -	2 500 56.59 26.18 273.87 -	2 1000 115.69 82.53 857.74 -	2 total 211.34 116.97 1,215.54 -	5 10 1.11 0.13 1.46 3,600.77(1)	5 20 2.03 0.36 3.78 n/a(20)	5 50 5.10 1.41 14.74 n/a(20)	5 100 10.34 5.72 47.38 -	5 200 20.85 15.15 145.10 -	5 500 55.02 82.17 834.01 -	5 1000 119.42 229.58 2,357.48 -	5 total 213.87 334.52 3,403.95 -	10 10 1.21 0.34 3.36 3,602.44(15)	10 20 2.04 0.81 8.00 n/a(20)	10 50 5.11 3.46 35.50 3,596.06(18)	10 100 10.31 13.30 148.41 -	10 200 21.09 46.59 497.50 -	10 500 56.17 195.58 1,999.56 -	10 1000 170.69 398.56 4,219.18 -	10 total 266.62 658.64 6,911.51 -	total total 691.83 1,110.13 11,531.00 -

Table 17 :

 17 MNLKP, linear weights, dissimilar capacities. Average solution values over 20 instances (# no solution).

	10 20	10 10	5 total	5 1000	5 500	5 200	5 100	5 50	5 20	5 10	2 total	2 1000	2 500	2 200	2 100	2 50	2 20	2 10	m n
	816.72	308.99	76,737.59	41,112.47	20,598.05	7,991.03	3,906.82	1,985.20	749.42	394.60	73,725.86	39,193.30	19,806.13	7,856.57	3,859.95	1,916.44	760.08	333.39	DMM
	814.43 854.19 n/a(20)	279.20 310.72 313.65(11)	70,679.55 72,806.91 -	37,430.50 38,623.60 -	19,076.70 19,540.10 -	7,535.88 7,739.56 -	3,657.64 3,776.51 -	1,907.99 1,965.24 n/a(20)	699.44 757.70 600.46(18)	371.40 404.20 434.82(10)	66,316.08 67,985.77 -	35,220.40 35,900.90 -	17,817.50 18,240.70 -	7,044.65 7,263.35 -	3,451.60 3,625.93 -	1,789.24 1,854.16 n/a(20)	696.13 758.70 674.85(16)	296.56 342.03 342.20(3)	Ipopt_1 Ipopt_10 Couenne	Real Variables
	817.87 824.66(1) 844.13 647.17	307.06 283.12 304.86 273.38(1)	76,769.80 ---	41,144.61 n/a(20) n/a(20) -	20,597.93 n/a(20) n/a(20) -	8,000.10 6,768.74(19) 7,021.51(18) -	3,899.92 3,657.99(8) 3,733.00(9) -	1,986.83 1,953.95 1,978.85 1,030.52(8)	746.29 708.66 739.91 481.27(6)	394.12 379.34 401.04 378.30(3)	73,781.80 ---	39,183.30 31,516.00(19) n/a(20) -	19,870.13 19,164.10(18) 18,296.70(14) -	7,854.79 7,209.46(2) 7,348.76(1) -	3,864.18 3,614.56 3,696.08 -	1,920.41 1,866.97 1,887.10 1,407.08(5)	757.24 705.37 763.36 652.22(1)	331.75 312.05 339.32 336.78	DMM Bonmin_1 Bonmin_10 Couenne	Integer Variables

Table 18 :

 18 MNLKP, linear weights, dissimilar capacities. Average CPU times over 20 instances (# no solution).

	Integer Variables	DMM Bonmin_1 Bonmin_10 Couenne	1.07 2.24 1.85 2,017.11	2.02 9.99 10.64 3,601.04(1)	5.08 530.30 137.54 3,602.40(5)	10.21 1,699.74 1,571.01 -	20.77 2,655.79(2) 2,601.08(1) -	58.30 1,493.78(18) 4,318.67(14) -	118.16 50,892.30(19) n/a(20) -
	Real Variables	m n Couenne DMM Ipopt_1 Ipopt_10	0.04 0.58 1613.09(3) 2 10 1.06	0.12 1.46 3,602.01(16) 2 20 2.02	n/a(20) 0.48 5.78 2 50 5.06	-1.55 16.23 2 100 10.21	-5.13 52.37 2 200 20.69	-22.90 237.33 2 500 56.61	-70.27 741.25 2 1000 115.26	2 total 210.91

Table 19 :

 19 MNLKP, nonlinear weights. Solution values for instances globally solved by Couenne.

	19 218.75 227.50	18 358.20 362.85	17 388.22 401.22	16 244.25 244.25	15 --	14 374.85 390.88	13 328.83 353.59	12 291.93 305.05	11 355.74 396.79	10 314.74 330.25	9 --	8 --	7 286.71 304.76	6 --	5 553.69 584.99	4 309.95 317.61	3 373.14 381.21	2 338.81 349.23	1 462.02 470.98	2 10 0 223.20 230.69	m n instance DMM Couenne	Real
																							Similar
	217.83 227.50	357.65 362.05	391.21 401.16	244.25 244.25	440.29 440.79	375.30 389.56	328.73 352.30	294.18 302.82	353.70 395.00	314.74 326.33	440.18 454.90	443.08 455.53	285.50 305.06	286.42 303.57	553.92 583.28	--	368.56 379.98	--	457.49 471.87	205.37 228.65	DMM Couenne	Integer
	227.50 227.50 227.50 227.50	332.11 359.07 330.58 357.38	350.11 352.39 351.32 352.11	244.25 244.25 244.25 244.25	409.60 428.26 411.79 422.74	361.39 390.92 362.07 389.32	329.20 349.62 329.08 349.50	298.81 299.82 292.96 299.39	374.36 395.14 369.83 393.58	327.94 330.40 322.26 326.86	--407.35 455.29	--438.36 455.49	217.68 221.78 217.83 221.64	250.55 264.43 250.55 263.71	553.66 587.48 554.02 585.63	309.04 316.53 312.20 316.17	367.99 382.72 367.83 382.82	340.36 349.29 339.04 349.21	--457.04 471.50	211.03 230.62 211.98 228.33	DMM Couenne DMM Couenne	Real Integer	Dissimilar

Translation: "It is impossible to separate a cube into two cubes, or a fourth power into two fourth powers, or in general, any power higher than the second, into two like powers. I have discovered a truly marvelous proof of this, which this margin is too narrow to contain". (see[114, p. 144-145]).

Available at URL http://miplib.zib.de/ Available at URL http://people.brunel.ac.uk/~mastjjb/jeb/orlib/portinfo.html Available at URL http://www.di.unipi.it/optimize/Data/MV.html Available at URL http://w3.uniroma1.it/tardella/datasets.html

A C K N O W L E D G M E N T S

The last four years have been an unforgettable experience for me and actually this is due to many persons I had the chance to meet. First of all I want to thank my supervisors for their incredible help and support also during the moments of greater difficulty.

as heuristics for non-convex problems like MNLKPs.

The first test-bed with respect to which we computationally test the heuristic algorithm is the same as in Section 4.2.3. Furthermore, we consider a second test-bed, characterized by linear weight functions, namely

with w j uniformly random in [START_REF] Alizadeh | Second-Order Cone Programming[END_REF][START_REF] Gablonsky | Modifications of the DIRECT Algorithm[END_REF].

The number of items n took the value in {10, 20, 50, 100, 200, 500, 1000}, while, as in the Section 4.2.3, the number of knapsacks m took the value in {2, 5, 10}. For each test-bed, we generate two groups of instances Similar capacities and Dissimilar capacities (see Section 4.2.3). For each combination of (n, m), 20 real instances, i.e., with N = ∅, and 20 integer instances, i.e., with N = N, are generated. The total number of tested instances was thus 3,360.

All the experiments were performed on an Intel Core 2, CPU 6600, 2.4 GHz, 1.94 GB ram, using only one processor.

We ran all the algorithms (Ipopt, Bonmin, Couenne, and DMM) with a time limit of one CPU hour per instance. Couenne was executed with its default values (with an exception mentioned at the end of the present section) as the use of other options strongly increases its computing times. Couenne was executed only for n 50, since already in those smaller instances its performances are quite poor. For DMM, we compute the solution of the surrogate relaxation (4.11) by means of Couenne with a time limit of n/10 seconds. If Couenne did not find a feasible solution within the time limit, we did not execute the surrogate heuristics. Moreover, the time limit of each local search (Step 6) is 5 CPU seconds.

We compare the heuristic against Bonmin and Ipopt with one and ten random starting points. We also tested Bonmin with ten random starting points at each branch-and-bound decision node, but this only resulted in few improvements for the small instances, and, however, in high computational time, so we do not report the corresponding computational results, which, however, can be found in the technical report [START_REF] Mencarelli | Relaxations and Heuristics for the General Multiple NonLinear Knapsack[END_REF]. Preliminary computational tests were conducted by using Scip [223], as exact solver for non-convex MINLPs, resulting in significantly worse performances than the ones obtained by Couenne.

Procedure Construct(i) (within Constructive) was executed for s ∈ {1, 10, 50, 100}, and the best solution was selected. The refined search for μj (see Section 4.4.1) was obtained: (i) by trying up to 5 consecutive refinement rounds, each time dividing the current sampling step by 2; (ii) by trying a single refinement round twice (dividing the initial sampling step by 5 and 10, respectively), and (iii) taking the best solution. We set the value of ε mean return vector of the assets ξ ∈ R r random vector of expected returns µ ∈ R r mean of the r-variate distribution of ξ ψ ∈ R r normalized portfolio return x, x ∈ R r lower and upper bound for the fraction of the portfolio value S ∈ R r size of the batches of the assets q ∈ R r + market value of the quoted assets x (0) ∈ R r fraction of the portfolio value already invested

asset-specific (non-factor) returns vector Ξ ∈ R r×r sensitivity-factor matrix Σ ∈ R r×r covariance matrix of the r-variate distribution of ξ Σ ∈ R r×r covariance return matrix of the assets

Sets:

L l ⊆ {1, . . . , r} set of assets for economic sector l (l = 1, . . . , ) C k ⊆ {1, . . . , r} set of indexes of all risky assets connected with the category k (k = 1, . . . , n)

Decision Variables:

x ∈ R r (fraction of the) portfolio value invested per asset η ∈ Z r + integer multiple of the lot-size S δ ∈ {0, 1} r additional binary variables such that δ j = 1 (j ∈ {1, . . . , r}) iff x j > 0 (j ∈ {1, . . . , r}) γ ∈ Z r + additional vector of general integer variables ζ ∈ {0, 1} n additional binary variables such that

additional binary variables such that y l = 1 (l ∈ {1, . . . , }) iff j∈L l δ j = 1 z ∈ {-1, 0, 1} r additional ternary variables such that z j = 1 (j ∈ {1, . . . , r}) iff x j > 0 (j ∈ {1, . . . , r}) and z j = -1 (j ∈ {1, . . . , r}) iff x j < 0 (j ∈ {1, . . . , r}) s ∈ R r additional continuous variables such that s j = √ x j (j ∈ {1, . . . , r}) b ∈ R T additional continuous variables such that b t = r j=1 (ν jtµ j )x j (t ∈ {1, . . . , T }) Functions: F (x) : R → R cumulative distribution of the normalized portfolio return φ : R → R penalty function for the minimum buy-in threshold constraint (3.9) θ : R → R penalty function for the round lot purchasing constraints (3.11)