Operations Research Letters 38 (2010) 39–46 Contents lists available at ScienceDirect Operations Research Letters journal homepage: www.elsevier.com/locate/orl Piecewise linear approximation of functions of two variables in MILP models Claudia D’Ambrosio, Andrea Lodi , Silvano Martello DEIS, University of Bologna, Viale Risorgimento 2, 40136 Bologna, Italy article info Article history: Received 1 February 2009 Accepted 24 July 2009 Available online 6 October 2009 Keywords: Mixed integer linear programming Non-linear functions of two variables Piecewise linear approximation abstract We consider three easy-to-implement methods for the piecewise linear approximation of functions of two variables. We experimentally evaluate their approximation quality, and give a detailed description of how the methods can be embedded in a MILP model. The advantages and drawbacks of the three methods are discussed on numerical examples. © 2009 Elsevier B.V. All rights reserved. 1. Introduction In recent years, the increased efficiency of mixed integer linear programming (MILP) software tools has encouraged their use also in the solution of non-linear problems, bringing to the need for ef- ficient techniques to linearize non-linear functions of one or more variables. The standard methodologies consist in the piecewise lin- ear approximation of such functions. For functions of a single variable, say, f (x), the piecewise linear approximation is obtained by introducing a number n of sampling coordinates x 1 ,..., x n on the x axis (breakpoints) on which the function is evaluated, with x 1 and x n coinciding with the left and right extremes of the x domain (see Fig. 1). The function is then approximated by the linear segments [(x i , f (x i )), (x i+1 , f (x i+1 ))] (i = 1,..., n 1). More precisely, for any given x value, say, ¯ x, with x i ≤¯ x x i+1 , the function value is approximated by convex combination of f (x i ) and f (x i+1 ). Let λ be the (unique) value in [0, 1] such that ¯ x = λx i + (1 λ)x i+1 . (1) Then the approximated value is f a ( ¯ x) = λf (x i ) + (1 λ)f (x i+1 ). (2) This methodology can alternatively be described through the slope (f (x i+1 ) f (x i ))/(x i+1 x i ) of the interpolating function, namely f a ( ¯ x) = f (x i ) + ( ¯ x x i ) f (x i+1 ) f (x i ) x i+1 x i . (3) From (3), one has λ = (x i+1 −¯ x)/(x i+1 x i ). Corresponding author. E-mail addresses: c.dambrosio@unibo.it (C. D’Ambrosio), andrea.lodi@unibo.it (A. Lodi), silvano.martello@unibo.it (S. Martello). In order to use the above technique in a MILP solver, it is nec- essary to include in the model variables and constraints that force any x value to be associated with the proper pair of consecutive breakpoints (or with a single one, in case x ∈{x 1 ,..., x n }). Let us introduce a continuous variable α i for each breakpoint i, such that α i ∈[0, 1] (i = 1,..., n). Let h i be a binary variable associated with the ith interval [x i , x i+1 ] (i = 1,..., n 1), with dummy val- ues h 0 = h n = 0 at the extremes. The approximate value f a can then be obtained by imposing the following constraints: n1 i=1 h i = 1 (4) α i h i1 + h i (i = 1,..., n) (5) n i=1 α i = 1 (6) x = n i=1 α i x i (7) f a = n i=1 α i f (x i ). (8) Constraint (4) imposes that only one h i , say, h ¯ ı , takes the value 1. Hence, constraints (5) impose that the only α i values different from 0 can be α ¯ ı and α ¯ ı+1 . It follows from (6) and (7) that α ¯ ı = λ and α ¯ ı+1 = 1 λ (see (1)). Constraint (8) ensures then the correct computation of the approximate value according to (2). In contexts of this type, the MILP constraints can be simplified by the so-called special ordered sets, introduced by Beale and Tomlin [1], and extensively studied in [2–4]. By defining a set of variables to be a Special Ordered Set of type k (SOSk), one imposes that at most k such variables can take a non-zero value, and that they must be adjacent. Most modern MILP solvers are capable of 0167-6377/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.orl.2009.09.005