Sandra Ulrich Ngueveu
email: ngueveu@laas.fr

Piecewise linear bounding of univariate nonlinear functions and resulting mixed integer linear programming-based solution methods

Keywords: OR in energy, Nonlinear programming, Combinatorial optimization, Piecewise linear bounding 2010 MSC: 90C30, 90C59, 90C90, 68W25

Piecewise linear bounding of univariate nonlinear

functions and resulting mixed integer linear programming-based solution methods.

Introduction

Various optimization problems result from the introduction of nonlinear terms into combinatorial optimization problems and can therefore be modeled as mixed integer nonlinear problems (MINLP). In the context of energy optimization, such nonlinear terms model energy conversion or demand/cost functions. Let us consider, for example, energy optimization in hybrid electric vehicles (HEV). In such vehicles the electrical powertrain system has multiple energy sources that it can gather power from to satisfy the propulsion power requested by the vehicle at each instant. The problem usually consists in finding at each instant the optimal power split between the multiple energy sources to satisfy the power demand of a driver on a predefined road section. The variant of the problem called offline assumes that the power demand profile is known a priori. The objective is to minimize the total fuel consumption of the vehicle performing a predefined mission, taking into account the characteristics and the limitations of each energy source, such as the energy losses happening during any energy transfer. Let (P) denote such problem for a HEV operating with two energy sources: (1) a Fuel Cell (FC) stack able to produce power from P 1 min up to P 1 max at each instant i ∈ {1...I}, (2) a Storage Element (SE) able to retrieve power up to P 2 min and to produce power up to P 2 max at each instant i ∈ {1...I}. The amount of energy stored in the SE is also called State Of Charge (SOC). To avoid a premature aging of SE, its state of charge is only allowed to vary between E min and E max , typically 25% and 100% of its energy capacity. Problem (P) can be modeled with equations ( 1)-( 5) where x 1 i , x 2 i and x 3 i are the amount of energy produced by the FC, produced by the SE, and retrieved by the SE at instant i ∈ {1...I}. Problem (P) is a (MI)NLP because of nonlinear energy conversion functions f 1 , f 2 , and f 3 continuous on intervals [P 1 min , P 1 max ], ]0, P 2 max ], and ]0, -P 2 min ] respectively, and verifying f 1 (0) = f 2 (0) = f 3 (0) = 0, often with a discontinuity at P 1 min or 0 that requires the introduction of binary variables to be modeled. The mathematical model can be expressed as follows:

(P) min

I i=1 f 1 (x 1 i ) //minimize total cost on FC (1)
subject to (s.t.)

x 1 i + x 2 i -x 3 i ≥ P i , ∀i ∈ {1...I} //power demand satisfaction (2)

I i=1 (f 2 (x 2 i ) -f 3 (x 3 i )) ≤ 0 //final SOC ≥ initial SOC (3) E 0 -E max ≤ i k=1 (f 2 (x 2 k ) -f 3 (x 3 k )
) ≤ E 0 -E min , ∀i ∈ {1...I} //SOC limits (4)

x 1 i ∈ {0} ∪ [P 1 min , P 1 max ], x 2 i ∈ [0, P 2 max ], x 3 i ∈ [0, -P 2 min ] ∀i ∈ {1...I} //domain definition. (5)
General MINLPs are NP-hard, but some subclasses such as convex MINLP may be easier to solve although still NP-hard. Convex MINLPs have a convex objective function to minimize, and their constraint functions are all convex and upper bounded. In these instances, when the integrity constraints are relaxed the feasible set is convex, and there exist efficient algorithms to solve the resulting convex NLP. However, the (MI)NLP modeling an energy optimization problem may be neither convex (nor concave) even when all energy conversion functions are convex or concave. Therefore, only small instances may be tractable using standard MINLP solvers. Several real-world applications have been addressed using piecewise linear approximations of the nonlinear functions of the MINLP to obtain a MILP which is easier to solve (see for example [START_REF] Borghetti | An MILP approach for short-term hydro scheduling and unit commitment with head-dependent reservoir[END_REF], [START_REF] Boukouvala | Global optimization advances in mixed-integer nonlinear programming, MINLP, and constrained derivative-free optimization, CDFO[END_REF], [START_REF] Camponogara | Compressor scheduling in oil fields: A piecewise-linear formulation[END_REF], [START_REF] D'ambrosio | Piecewise linear approximation of functions of two variables in MILP models[END_REF]). Such approach presents the main advantage of producing solutions faster than purely MINLP-based approaches if not too many additional

(v3)

Technical report LAAS-CNRS n • 16358 2/ 29 binary variables or combinatorial conditions have been introduced in the process, meaning that the number of pieces of the piecewise linear (pwl) functions used should be limited. [START_REF] Geißler | Using Piecewise Linear Functions for Solving MINLPs[END_REF] explain that the approach suffers from a few drawbacks because the nature of the nonlinearities inherent to the problem may be lost and the solutions obtained from the MILP may have no meaning for the MINLP. If the solution obtained is not satisfactory in terms of feasibility or optimality, then a new pwl approximation may be performed using a higher number of pieces, to obtain a new MILP to solve. This yields an iterative solution procedure with a number of iterations unknown a priori which translates into high computing times, either because several iterations needed to be performed, or because an unnecessary large number of pieces were introduced upfront, resulting into an unnecessarily large MILP that required a high solution time.

As an alternative we propose a straightforward two-phase solution method. The first phase consists in bounding each nonlinear term from above and below using a pair of piecewise linear functions satisfying conditions that will be specified in the core of the paper. Contrary to most publications on piecewise linear approximation which focus on the minimization of the approximation error for a given number of pieces or breakpoints that may or may not be equidistant, we aim at minimizing the number of pieces for a given error bound. The second phase of the solution method consists in solving two MILPs obtained from the replacement of the nonlinear functions with one of the two piecewise linear functions. This paper addresses the first phase of the iterative procedure, assuming that the resulting MILPs can be solved efficiently with a MILP solver. If it is not the case, then a specific solution method may need to be designed for them. An example of such case is presented in [START_REF] Ngueveu | Scheduling under a non-reversible energy source: An application of piecewise linear bounding of non-linear demand/cost functions[END_REF], which focused on the solution of the resulting MILP without specifying how to obtain the piecewise linear functions.

State of the art

Several publications exists on the application of piecewise linear (pwl) approximation on nonlinear univariate functions to solve MINLP problems, but the issue addressed in the large majority of them is to minimize the approximation error given a predefined number of breakpoints or pieces. To the best of our knowledge, only three papers focused on the specific problem of minimizing the number of breakpoints for a given tolerance or bounded approximation error. [START_REF] Rosen | Global minimization of large-scale constrained concave quadratic problems by separable programming[END_REF] were the first to propose the computation of breakpoints for a given error tolerance for concave quadratic functions. The pwl interpolators were built using equidistant breakpoints, and because the input functions were concave their interpolators were underestimators. They identified conditions on the number of breakpoints required to achieve a given error tolerance. [START_REF] Geißler | Using Piecewise Linear Functions for Solving MINLPs[END_REF] showed that certain cases of general MINLPs can be solved by just applying techniques purely from the mixed integer linear programming by approximating the nonlinearities with pwl functions. They proposed to compute a priori errors for pwl approximations or a priori errors for over-and under-pwl estimators. However, they did not focus on the computation of optimal (minimal) breakpoint positions. [START_REF] Rebennack | Continuous piecewise linear delta-approximations for univariate functions: Computing minimal breakpoint systems[END_REF] propose two exact approaches and two heuristics for the computation of optimal continuous pwl approximators for univariate continuous functions

2018 (v3) Technical report LAAS-CNRS n • 16358 3/ 29 0 1 2 3 4 5 0 1 2 continous δ-approximator maximizing interval length (unique) optimal continuous δ-approximator f (x) δ-tube around f (x) x f (x) f (x) =        1 if x ∈ [0, 2[ 0.75x -0.5 if x ∈ [2, 3[ -0.5x + 3.25 if x ∈ [3, 4[ 0.5x -0.75 if x ∈ [4, 5] (6)
Figure 1: Maximizing the length of the intervals successively is not optimal, in general (source: [START_REF] Rebennack | Continuous piecewise linear delta-approximations for univariate functions: Computing minimal breakpoint systems[END_REF], page 628) over a compactum D = [X -, X + ]. Their algorithms handle more general functions than the ones of [START_REF] Rosen | Global minimization of large-scale constrained concave quadratic problems by separable programming[END_REF]. In addition, their breakpoints are distributed freely and shifts from the function are allowed at breakpoints, which were shown to be important degrees of freedom contributing to a significant reduction of the number of breakpoints, up to 80% in some cases. The work of [START_REF] Rebennack | Continuous piecewise linear delta-approximations for univariate functions: Computing minimal breakpoint systems[END_REF] differs from [START_REF] Geißler | Using Piecewise Linear Functions for Solving MINLPs[END_REF] in the following aspects. [START_REF] Geißler | Using Piecewise Linear Functions for Solving MINLPs[END_REF] do not target on computing the minimal number of breakpoints and do not consider shift variables at breakpoints whereas [START_REF] Rebennack | Continuous piecewise linear delta-approximations for univariate functions: Computing minimal breakpoint systems[END_REF] do so. Since their work is close to ours in some aspects, let us focus on their contributions before highlighting the main differences and contributions of this paper. [START_REF] Rebennack | Continuous piecewise linear delta-approximations for univariate functions: Computing minimal breakpoint systems[END_REF] show that ensuring that the approximator and the original function do not deviate more than a predefined tolerance δ from each other leads to sets of constraints which have to hold over a continuum, resulting in a semi-infinite programming (SIP) problem denoted OBSC for "Optimal Breakpoint System using a Continuum approach for x". The authors show that it is NP-hard to compute a δ-approximator for an arbitrary continuous function and propose an iterative solution procedure based on the evaluation of the continuum conditions only on a discrete set of grid points, resulting into a MINLP model which is a relaxation of the SIP. The feasibility of the resulting solution with regards to OBSC is then evaluated by solving an NLP on each interval corresponding to a line-segment of the pwl approximator, to compute the true maximal deviation between the approximator and the original function. The algorithm stops if the true deviation is less or equal to δ on all line-segments. In this case the solution obtained is optimal for OBSC. Otherwise the grid is refined to obtain a new MINLP to be solved. OBSC and the MINLPs are in general too large and difficult to solve to optimality, therefore the authors proposed two heuristic methods. The heuristics methods were based on the successive computation of the breakpoints, from X -to X + , maximizing at each iteration the length of the interval on x. This meant solving at each iteration a problem denoted BSB that can be expressed as follows: given the breakpoint x i ending the i th line-segment (which corresponds by continuity to the beginning of the i + 1 th line-segment), compute the next breakpoint x i+1 (end of the i + 1 th line-segment) so as to maximize x i+1 while ensuring a deviation of at most δ between the original function and the linear approximation on interval

[x i , x i+1 ]. Rebennack 2018 (v3) Technical report LAAS-CNRS n • 16358 4/ 29
and Kallrath (2015) provided a counter-example showing that maximizing the length of the intervals do not necessarily lead to an optimal breakpoint system i.e., to a δ-approximator with the least number of breakpoints, contrary to what is stated for example in [START_REF] Frenzen | On the number of segments needed in a piecewise linear approximation[END_REF]. The counter-example is illustrated on Figure 1. It shows, for a function f (x) that the authors proposed with [X -, X + ] = [0, 5] and δ = 0.25, that the unique optimal 0.25approximator uses two line-segments whereas maximizing the interval length successively from X -to X + produces a 0.25-approximator using three line-segments. The Forward Heuristic with Moving Breakpoints (FHMB) solved each BSB problem to optimality with the iterative "grid discretization + NLP solution" approach. Its main limitation was the necessity to solve several NLPs. The α-Forward Heuristic with Backward Iterations (FHBI) solved each BSB problem heuristically by trying different decreasing values for x i+1 with a predefined step parameter α. FHBI solved less NLPs than FHBM, and therefore required less computing time, but FHBI obtained better solutions. Using any of the heuristics, it is possible to obtain breakpoint positions satisfying the required δ-tolerance, and provide an upper bound on the minimal number of breakpoints.

It is worth mentioning that there exists publications on piecewise linear approximation with a minimum number of pieces given a predefined bound on the absolute error in the fields of data reduction, pattern recognition or classification, and ECG waveform preprocessing (Tomek (1974a), [START_REF] Tomek | Two algorithms for piecewise-linear continuous approximation of functions of one variable[END_REF], [START_REF] Gritzali | A fast piecewise linear approximation algorithm[END_REF]). The main difference with our problem is that those publications consider as an input a discrete set of points. Their objective is to find the piecewise linear function with a minimum number of pieces, such that the error for each of the sample points is less than the allowed value δ. Even in cases where an analytical expression of a continuous function was available, the function was sampled and the approximation was performed on the set of sample points. The algorithms proposed in these research fields did not ensure the respect of the predefined approximation error on the entire interval [X -, X + ] and are therefore not applicable to our problem. Finally, the field of piecewise linear approximation of planar curves could be mentionned [START_REF] Dunham | Optimum uniform piecewise linear approximation of planar curves[END_REF], [START_REF] Papakonstantinou | Parallel approaches to piecewise linear approximation[END_REF]), with applications related to shape analysis or pattern classification, since a nonlinear function f (x) could be represented as a parametric curve x(t) = t, y(t) = f (t). However, the algorithms proposed in that research field are not applicable to our problems for two main reasons: (i) input curves considered are discrete or digitised, not continuous ones as ours, and (ii) the error metric is the Euclidean distance, which does not correspond to the approximation error between the original function and the piecewise linear function.

In the view of the state-of-the-art, the contributions of the paper are the following: (1) instead of using continuous pwl δ-approximators, we propose to approximate general univariate continuous functions with non necessarily continuous piecewise linear δ-approximators, adding an additional degree of freedom to obtain a breakpoint system of equal or less linesegments, (2) we prove that when discontinuity is allowed, maximizing the interval length produces an optimal pwl δ-approximator, leading to an exact solution procedure based on the iterative solution of adapted BSB problems, (3) we introduce relative -tolerance and show the benefit of using it instead of the absolute δ-tolerance, (4) models and algorithms to compute the discontinuous pwl under-and over-estimator with absolute or relative tolerance and with an additive worst case guarantee are presented, (5) for solving MINLP involving nonlinear univariate energy conversion functions, we propose a solution method based on the

(v3)

Technical report LAAS-CNRS n • 16358 5/ 29 upper and lower bounding of energy conversion expressions using non necessarily continuous piecewise linear functions with a relative -tolerance, and the solution of a pair of mixed integer linear programs, (6) we show that such approach yields a performance guarantee when the nonlinearity is restricted to the objective function, and, finally (7) computational results on energy optimization problems for hybrid electric vehicles illustrate the efficiency of the method in comparison to state-of-the-art methods, including solution procedures based on approximations with absolute δ-tolerance.

3 Non necessarily continuous pwl δ-approximation of continuous nonlinear functions

Let f : D = [X -, X + ] → R be a function on the compact interval D ⊂ R. A function g : D = [X -, X + ] → R is a pwl function with n g ∈ N line-segments if ∃a ∈ R ng , b ∈ R ng , x min ∈ [X -, X + ] ng and ∀i ∈ [1...n g ], ∃x max i ∈]x min i , X + ],
such that the following equations are verified:

g(x) = a i x + b i , ∀i ∈ [1...n g ], ∀x ∈ [x min i , x max i ] (7)
x max i = x min i+1 , ∀i ∈ [1...n g -1] (8)
x

min 1 = X - (9)
x max ng = X + (10) Such pwl function is said to be defined by

G = ng i=1 ([a i , b i ], [x min i , x max i ])
, and the two endpoints x min i and x max i of each line-segment i are called breakpoints. A pwl function g is:

• continuous iff it verifies all following equations:

a i x max i + b i = a i+1 x min i+1 + b i+1 , ∀i ∈ [1, n g -1] (11)
• or discontinuous if it does not verify all equations (11), i.e. ∃j ∈ [1, n g -1] such that

(a i x max i + b i ) = (a i+1 x min i+1 + b i+1 )
• or non necessarily continuous (nnc) if the satisfaction of equations ( 11) is neither required nor forbidden

Definition 3.1 (δ-approximator): A pwl function g : D = [X -, X + ] → R is called δ- approximator of a function f : D = [X -, X + ] → R with δ > 0, iff : max x∈D |g(x) -f (x)| ≤ δ (
2018 (v3) Technical report LAAS-CNRS n • 16358 7/ 29 Proof Let l = n l i=1 ([a i , b i ], [x min i , x max i ]
) be an optimal continuous δ-pwl approximator (n l = n * ) of a continuous function f : D = [X -, X + ] → R. Let y ∈ D be the solution value of the problem defined as follows:

(PY)        max y //maximize ending breakpoint s.t. |a 1 x + b 1 -f (x)| ≤ δ, ∀x ∈ [x min 1 , y] //δ-approximation constraint y ∈ [x min 1 , x max i ] //domain definition Let q ∈ [1...n l ] be the piece number that verifies x min q ≤ y ≤ x max q . Breakpoint x max 1 verifies |a 1 x max 1 + b 1 -f (x)| ≤ δ, therefore y ≥ x max 1 . If y = x max 1
, then l = l is a pwl continuous δapproximator of f with n = n * pieces where the projection of the first line-segment on interval D is of maximal length. Otherwise y > x max 1 , in which case a discontinuous pwl δ-approximator l with n l ≤ n * line-segments where the projection of the first line-segment on interval D is of maximal length can be defined by L =

n l i=1 ([ a i , b i ], [ x min i , x max i ]) = ([a 1 , b 1 ], [x min 1 , y]) ∪ ([a q , b q ], [y, x max q ]) ∪ n l i=q+1 ([a i , b i ], [x min i , x max i ]) .
Theorem 3.3 For any continuous function f : D = [X -, X + ] → R and any scalar δ > 0, there exists an optimal non necessarily continuous pwl δ-approximator g defined by

G = ng i=1 ([a i , b i ], [x min i , x max i ]), such that each line-segment i has a maximal length projection on the interval [x min i , X + ].
Proof There exists a continuous δ-approximator function for any continuous function f on a compactum D and any scalar δ > 0 (Duistermaat and Kol, 2004). Proposition 3.2 can be applied iteratively on an optimal continuous δ-approximator to obtain an optimal non necessarily continuous pwl δ-approximator with intervals of maximal length on D.

Figure 2 illustrates the three possible cases related to Proposition 3.2. Figure 3(a) illustrates the implications of Theorem 3.3 for the function f defined with equation ( 6): maximizing the length of intervals leads to a discontinous pwl 0.25-approximator of two line-segments, which is the optimal number of line-segments for any non necessarily continuous pwl 0.25approximator of f . In this example optimal continuous and non necessarily continuous pwl approximators have the same number of line-segments. Figure 3(b) illustrates an example where it is not the case: one has more line-segments than the other. In the general case, it can be easily proven that given a continuous function f : D = [X -, X + ] → R and given a scalar δ ∈ R + , if n nnc is the number of line-segments of an optimal non necessarily continuous pwl δ-approximator of f and if n c is the number of line-segments of an optimal continuous pwl δ-approximator of f , then nc+1 2 ≤ n nnc ≤ n c . The mathematical models and algorithms from [START_REF] Rebennack | Continuous piecewise linear delta-approximations for univariate functions: Computing minimal breakpoint systems[END_REF] for computing continuous pwl δ-approximators are modified to compute non necessarily continuous pwl δ-approximators by considering two shifts per breakpoint instead of one, namely s - b used for the line-segment ending a breakpoint b and s + b used for the line-segment starting at breakpoint b. In particular, thanks to Theorem 3.3, Algorithm 1 resulting from the adaptation of Forward Heuristic with Moving Breakpoints [START_REF] Rebennack | Continuous piecewise linear delta-approximations for univariate functions: Computing minimal breakpoint systems[END_REF] is a now exact method for computing optimal non necessarily continuous δ-approximators. Its main limitation is the necessity to solve an SIP, for example with an iterative grid and NLP-based procedure as described in Algorithm 2. The tightest approximator computation given the

2018 (v3) Technical report LAAS-CNRS n • 16358 8/ 29 0 1 2 3 4 5 0 1 2 nnc δ-approximator maximizing interval length (unique) optimal continuous δ-approximator f (x) δ-tube around f (x) x f (x) 0 1 2 3 4 5 6 0 1 2 nnc δ-approximator maximizing interval length (unique) optimal continous δ-approximator f (x) δ-tube around f (x) x f (x) (a) (b)
Figure 3: Maximizing the length of the intervals successively is optimal for non necessarily continuous δ-approximators

Algorithm 1 Optimal non necessarily continuous pwl approximation Input: Function f ; domain D = [X -; X + ]; scalar δ > 0 Output: pwl function g defined by G 1: n g := 0; y = X -; G := ∅; x end := X + 2: while y < X + do 3:

x begin := y 4:

solve the following SIP to obtain a, b, and y: max y //maximize ending breakpoint (13)

s.t. |ax + b -f (x)| ≤ δ, ∀x ∈ [x begin , y] //δ-approximation constraint (14) a ∈ R, b ∈ R, y ∈ [x begin , x end ] //domain definition (15) 5: G := G ∪ ([(a, b), (x begin , y)]) 6: n g := n g + 1 7: end while 2018 (v3) Technical report LAAS-CNRS n • 16358 9/ 29
Algorithm 2 Iterative grid and NLP-based procedure to solve SIP ( 13)-( 15)

Input
s.t. s 0 + f (x begin ) + i(sy+f (y)-f (x begin )-s 0 ) |N I |+1 -f (x i ) ≤ δ, ∀i ∈ N I
//δ-approximation on grid(17)

x i = x begin + i |N I |+1 (y -x begin ), ∀i ∈ N I //grid computation (18) s 0 ∈ [-δ, δ], s y ∈ [-δ, δ], y ∈ [x begin , x end ],
and

x i ∈ R, ∀i ∈ N I . //domain definition (19) 5: compute the gap := max x∈[x begin ,y] |ax + b -f (x)| 6: end while
optimal number breakpoints computed with Algorithm 1 would also require the solution of NLPs [START_REF] Rebennack | Continuous piecewise linear delta-approximations for univariate functions: Computing minimal breakpoint systems[END_REF]). Section 4 presents, among other things, an exact method applicable to convex and concave functions that does not need to solve any SIP or any NLP. The method also has an additive worst case guarantee for functions that are neither convex nor concave but that can be decomposed into convex or concave pieces, if the decomposition is part of the input.

4 From pwl δ-approximation to pwl δ-bounding A benefit of computing δ-approximators optimal in terms of number of line-segments is to minimize the number of additional binary variables added to obtain the MILP resulting from the replacement of nonlinear terms with their pwl approximators. A drawback of applying δapproximation is that the optimal solution of the resulting MILP can be infeasible with respect to the original MINLP. In this case the solution may not be of interest for the practitioners who formulated the original problem. A usual alternative is to use over-and under-estimators defined as follows.

Definition 4.1 (δ-underestimator): A pwl function l : D → R is a δ-underestimator of a function f : D → R with δ ∈ R + iff l(x) ≤ f (x) ≤ l(x) + δ, ∀x ∈ D. Definition 4.2 (δ-overestimator): A pwl function l : D → R is a δ-overestimator of a func- tion f : D → R with δ ∈ R + iff l(x) ≥ f (x) ≥ l(x) -δ, ∀x ∈ D.
In the context of energy optimization for hybrid electric vehicles for example, replacing each nonlinear energy loss function with its pwl δ-overestimator yields a MILP solution where a sufficient or excess amount of energy is produced at each time-step. Since excess energy can

(v3)

Technical report LAAS-CNRS n • 16358 10/ 29 be dissipated as heat (in a resistance inserted in the braking system or in mechanical brakes present for security reasons) therefore the MILP solution is applicable on a test bench or on the real world hybrid electric vehicle considered. Having ensured the feasibility of the solution obtained, it is of interest to provide an estimate of the quality of the solution with respect to the optimum of the original MINLP. This can be done using a pwl δ-underestimator of the energy loss function and solving the resulting MILP problem to obtain a lower bound of the MINLP. Each nonlinear function is therefore bounded from above and below with two pwl functions, yielding two MILPs whose optimal solutions cost verify:

z MILP ≤ z MINLP ≤ z MILP .
The resulting solution method is denoted PLB+MILP. Given an optimal continuous pwl δ-approximator g of a function f , an optimal continuous pwl 2δ-over-(resp. under-) estimator can be obtained from shifting g by δ (resp. -δ), as stated by [START_REF] Rebennack | Continuous piecewise linear delta-approximations for univariate functions: Computing minimal breakpoint systems[END_REF] i.e. g(x) + δ (resp. g(x) -δ) is a 2δ-over (resp. under-) estimator of f . This result can be extended to non necessarily continuous pwl approximators and estimators. In addition, in the case of energy optimization for example, the MINLP can be neither convex nor concave even though the individual nonlinear energy loss functions or energy demand/cost conversion functions are convex or concave. Yet, if a nonlinear continuous function f : D → R is convex (resp. concave) over D with a derivative efficiently computable at any point of D, then taking advantage of the following equation, it is possible to compute an optimal δ-underestimator (resp. δ-overestimator) without solving any SIP or even NLP:

f is convex over D ⇔ f (y) ≥ f (x) + f (x)(y -x), ∀x, y ∈ D (20)
Shifting can then be applied to obtain an optimal δ-overestimator (resp. δ-underestimator).

Let us focus, for example, on the computation of the optimal δ-underestimator of a nonlinear continuous convex function f with a derivative efficiently computable. The reasoning hereafter can be adapted to compute the optimal δ-overestimator of concave functions. Thanks to Theorem 3.3 the objective at each iteration i of the algorithm is to maximize x max i for a given x min i . Each line-segment i is tangent to f (otherwise it can be replaced with a line-segment of identical slope but tangent to f ). Therefore the objective at iteration i is to identify a tangent point q i ∈ [x min i , x max i ] that defines the slope a i = f (q i ) and the y-intercept b i = f (q i ) -f (q i )q i of line-segment i so that x max i is maximized given x min i . Given q i , the error in function of y defined as f (y)-(f (q i )y + f (q i )q i ) increases when y decreases if y ≤ q i , or when y increases if y ≥ q i . Therefore, for a given q i , checking whether the δ-approximation constraint is verified at points x min i and x max i is sufficient to ensure the validity of the constraint on the entire interval

[x min i , x max i ].
In this context, the following proposition can be enunciated:

Proposition 4.3 Finding the line-segment that maximizes x max i for a given x min i is equivalent to solving sequentially two separate problems: (1) maximizing q i for the given x min i , and then (2) maximizing x max i for the previously computed q i . Proof Let x max i (q (k) i ) be the maximal value possible of x max i for a given q

(k) i . It suffices to prove that x max i (q (2) i ) ≥ x max i (q (1) i ) ⇒ q (2) i ≥ q (1)
i . This is done using equation ( 20) which leads to 0 ≤ f (y) -(a

(2) i y + b (2) i ) ≤ f (y) -(a (1) i y + b (1) i ), ∀y ≥ q (2) i .

(v3)

Technical report LAAS-CNRS n • 16358 11/ 29

Algorithm 3 Optimal non necessarily continuous pwl δ-underestimator Input: Convex function f ; domain D = [X -; X + ]; scalar δ > 0 Output: pwl function g defined by G 1: n g := 0; y = X -; G := ∅ 2: while y < X + do 3:

x begin := y 4:

solve the following problem to get q and set a := f (q), b := f (q) -f (q)q (PQ)

       max q //maximize tangent point s.t. f (x begin ) -(f (q) + f (q)(x begin -q)) ≤ δ //δ-approximation constraint q ∈ [x begin , X + ] //domain definition 5:
solve the following problem to obtain y

(PY)        max y //maximize ending breakpoint s.t. f (y) -(ay + b) ≤ δ //δ-approximation constraint y ∈ [q, X + ] //domain definition 6: G := G ∪ ([(a, b), (x begin , y)]) 7:
n g := n g + 1 8: end while Algorithm 4 Solve problem (PQ) Input: Convex function f ; x begin ;X + ; δ > 0; α ∈ N Output: tangent point q, slope a, y-intercept b 1: {q, a, b, 0, 0}=Dichotomy_q(f ,x begin , X + ,δ,α,x begin , X + , 0, 0, 0) {use Algorithm 6}

Algorithm 5 Solve problem (PY) Input: Convex function f ;X + ; δ > 0; α ∈ N; tangent point q, slope a, y-intercept b Output: next breakpoint y 1: {q, 0, 0}=Dichotomy_y(f ,x begin , X + ,δ,α,q, X + , 0, 0, 0) {use Algorithm 7}

2018 (v3) Technical report LAAS-CNRS n • 16358 12/ 29
Algorithm 6 Dichotomy_q Input: Convex function f ; x begin ;X + ; δ > 0; α ∈ N; B min , B max ; q; a; b Output: tangent point q, slope a, y-intercept b, lower limit B min , upper limit

B max 1: if B max -B min ≤ 10 -α then 2: q = B min ; a = f (q); b = f (q) -f (q)q 3:
return q, a, b 4: else {y, B min , B max }=Compute_y(f , X + ,δ,α,q, a, b,B min , B max ) 11: end if

5: ã = f (x begin ); b = f (x begin ) -f (q i )q i 6: if f (x begin ) -(ãx begin + b) ≤ δ then 7: B min = B max ; B max = X + 8:

(v3)

Technical report LAAS-CNRS n • 16358 13/ 29

Algorithm 3 summarizes the resulting procedure that computes an optimal δ-underestimator of a nonlinear convex function f with a derivative efficiently computable at any point of [X -, X + ]. Problems (PQ) and (PY) can be solved to optimality with dichotomy search methods as illustrated by Algorithms 6 and 7 where parameter α is the number of significant decimals for q i , x min i , and x max i . Notice that there is no need to solve an SIP or NLP. It can be noted that any two consecutive line-segments

([a i , b i ], [x min i , x max i ]) and ([a i+1 , b i+1 ], [x min i+1 = x max i , x max i+1 ]) of the pwl function obtained verify a i = a i+1 and a i x max i + b i ≤ a i+1 x min i+1 + b i+1 . Their supporting straight lines intersect on a point C of coordinates (x C , y C ), such that y C = a i x C + b i = a i+1 x C + b i+1 and x min i ≤ x C ≤ x max i and f (x C ) -y C ≤ δ.
Consequently, the non necessarily continuous pwl δ-under-estimator computed with Algorithm 3 can be converted into a continuous pwl δ-under-estimator with the same number of line-segments.

Algorithm 3 can be adapted to compute directly an optimal δ-overestimator of a convex function: a tangent point q i in this case defines the slope a i = f (q i ) and the y-intercept b i = f (x min i ) -f (q i )x min i . It can also be adapted for the pwl bounding of concave functions with derivatives efficiently computable. Let us consider the case of continuous functions that are neither concave nor convex but that can be decomposed into p pieces, each piece being either concave or convex. If the decomposition is part of the input, then pwl bounding the function can be done by bounding each piece separately using the algorithms previously described, then aggregating the pwl functions obtained. To that end, let k i ∈ D, ∀i ∈ {1..p+1} be a set of points, such that (i) on each interval

[k i , k i+1 ] the function f (x) is either convex or concave, (ii) k 1 = X -, (iii) k p+1 = X + , and (iv) k i < k i+1 , ∀i ∈ {1..p}. For each interval [k i , k i+1 ]
, let g i be the optimal pwl over-(resp. under-) estimator of f computed with the algorithms from Section 4. Let g be an over (resp. under) estimator of function f on interval D, resulting from the union of pwl functions g i , i.e. g = ∪ i∈{1..p} g i . Let n g i (resp. n g ) be the number of line-segments of g i (resp. g). If n * is the number of line-segments of any optimal pwl over-(resp. under-) estimator of f on interval D, then n g = p i=1 n g i and the following equation is verified, providing an additive worst case guarantee:

n * ≤ n g ≤ n * + p -1 (21)
5 Drawbacks and limitations of pwl δ-bounding

When applying PLB+MILP, choosing relevant δ values is a challenging issue. The chosen value should be orders of magnitude smaller than the resulting solution cost to provide an acceptable precision level. After solving the MILP, verifying whether the chosen tolerance value δ was sufficiently small in the view of the resulting solution costs is straightforward. If it was not the case then the value of δ is decreased before a new round of pwl bounding then MILP solution. Such iterative procedure is not satisfactory because the number of iterations is unknown a priori and the pwl bounding and/or the MILP solution may require significant computational efforts during each single iteration. An alternative is to identify a target precision ∆ of the MILP solution with regards to the optimal MINLP solution value and precompute a priori the corresponding δ values in function of the chosen ∆. When nonlinearity occurs in the constraints, the link between δ and ∆ is not obvious. But even when nonlinearity occurs only in the objective function, δPLB+MILP with a target precision ∆ presents significant limitations described in the remainder of this Section.

(v3)

Technical report LAAS-CNRS n • 16358 14/ 29

Data dependence leading to multiple δ values

Let us consider the problem (1)-( 5) for a specific HEV, which, for simplicity, is assumed to have an ideal supercapacitator (f 2 (x) = x, f 3 (x) = x). Therefore the nonlinearity in the problem comes from the objective function which is composed of I univariate positive nonlinear terms f 1 (x). Ensuring a final precision of at least ∆ requires to bound f 1 with a tolerance δ = ∆/I. Even if the same target precision ∆ is chosen, different power profiles translate into different δ values. Two power profiles with the same time horizon I, may also require different δ values because different power profiles may translate into solution costs which can differ significantly in order of magnitude, requiring different ∆ values, and, therefore, different δ values.

To summarize, for any new power profile provided, the pwl bounding of the same function f 1 may need to be redone with a value of δ suitable to the new data set, even though the nonlinear terms themselves remained unchanged. In addition, since the δ errors on the nonlinear terms are additive, a longer time horizon means a higher number of univariate terms (f 1 (x i )) in equation ( 1), which means that a smaller tolerance δ may be required. This translates into an increase of the number pieces for the pwl functions, and, therefore, more binary variables which adds to the complexity of a MILP that was already penalized by the fact that long time horizons meant more decision variables x 1 i , x 2 i , x 3 i . These all contribute to a substantial reduction of the size of instances that can be solved efficiently.

Solution dependence leading to unknown δ values

There exists several cases of MINLP where knowing ∆ is not sufficient to infer the corresponding δ values. Let us consider, for example, problem (CF) modeling a scheduling problem with a single energy source. There are a set of activities A, each activity a having a release date r a , a due date d a , a duration p a , and an instantaneous energy demand b a . A constant term α at with a ∈ A and t ∈ T is equal to 1 if t ∈ [r a , d a [ and 0 otherwise. The efficency function of the energy source used to satisfy the total demand of activities scheduled at each time period is denoted ρ, i.e. a cost or energy consumption of ρ(x) produces an amount of usable energy x. Therefore ρ is defined on [0, a∈A b a ] and verifies ρ(0) = 0. The goal is to schedule the tasks so as to minimize the total energy cost. A resulting mathematical model requires binary decision variable x at that is equal to 1 iff activity a is active at instant t ∈ T . Continuous variables w t represent the total energy demand at instant t ∈ T . Nonlinearity comes from the objective function which is comprised of |T | univariate nonlinear terms. The resulting problem is: Let (CF) ∆ refer to the MILP derived from problem (CF) with a target tolerance ∆. Identifying a relevant δ value for obtaining (CF) ∆ is not straightforward. Indeed, in the general case, the duration of the schedule for the optimal solution is not known a priori. The only information available is the time horizon |T |, but this value may be far from the real ending time of the optimal schedule. As an illustration, consider an instance with a time horizon |T | = 100 having an optimal solution with an ending time of 10 (i.e. an optimal solution x * that verifies a∈A x * a,10 ≥ 1 and a∈A t>10 x * at = 0). In this case only 10 terms are active in the objective function, therefore bounding each term with a tolerance δ = ∆/10 should have been sufficient to achieve the requested global tolerance ∆. Instead, because the duration is unknown a priori, each term of the objective function has to be bounded with precision δ = ∆/100. In summary, a pwl bounding tolerance 10 times smaller than necessary would be requested, leading to pwl bounding functions with a higher number of line-segments, leading to a higher number of binary variables, leading to MILPs more time consuming to solve than necessary given the tolerance target ∆.

(CF) min t∈T ρ(w t ) //

Using -relative tolerance

To counter the drawbacks identified in Section 5, we propose to perform the upper and lower bounding of energy conversion expressions using non necessarily continuous pwl functions with a relative -tolerance. Definition 6.1 On a compactum D, pwl bounding a function f : D → R * with a relative tolerance value ∈]0, 1] consists in identifying two pwl functions (f , f ) that verify: The resulting PLB+MLP method shares the advantages of the δPLB+MILP method but, in addition, thanks to the use of relative tolerance, guarantees can be obtained on the quality of resulting bounds under conditions expressed in Proposition 6.2. Models and algorithms from Sections 3 and 4 can be adapted to ensure satisfaction of the relative tolerance constraints ( 27)-( 29) instead of the absolute δ tolerance constraints expressed in Definitions 4.1 and 4.2. Proposition 6.2 Let (P) be a (MI)NLP linearly constrained and with an objective function decomposable into a sum of univariate positive linear or nonlinear functions, i.e. verifying the following equation where g ki : D → R + , ∀k, ∀i are linear or nonlinear functions:

f (x) ≤ f (x) ≤ f (x), ∀x ∈ D (27) |f (x) -f (x)| ≤ |f (x)|, ∀x ∈ D (28) |f (x) -f (x)| ≤ |f (x)|, ∀x ∈ D (
(P) min or max z (P ) = g(x) = k i g ki (x i ) s.t. Ax ≤ B (30) 2018 (v3) Technical report LAAS-CNRS n • 16358 16/ 29
Let P (resp. P ) be the MILP resulting from the replacement of each nonlinear term g ki with its over-estimator g ki (resp. under-estimator g ki ). Then the solution values of the optimal solutions z (P ) , z (P ) , and z (P ) of the corresponding problems verify the following equations:

(1 -)z (P ) ≤ max z (P ) , z (P )

1 + ≤ z (P ) ≤ min z (P ) , z (P ) 1 - ≤ (1 + )z (P ) (31)
Proof It results from the combination of equations ( 27) with z (P ) ≤ z (P ) ≤ (1 + )z (P ) and

(1 -)z (P ) ≤ z (P ) ≤ z (P )

Let f (x) (resp. f (x) ) be an optimal -over-(resp. under-) estimator of f . Contrary to the case of absolute tolerance, if relative tolerance is used then f (x) -(resp. f (x) + ) is not necessarily an under-(resp. over-) estimator of f (x). Therefore is not possible to deduce an optimal -over-(resp. under-) estimator of f from applying a shift on an optimal -under-(resp. over-) estimator of f . Each estimator has to be computed separately.

Computational results

Models and algorithms run on an Intel(R) Xeon(R) CPU E3-1271 v3 computer with 32 GB RAM using a single core. PLB heuristic algorithms based on a dichotomic search (Algorithm 3) are coded in MATLAB R2017b 64bits. Their parameter α set to 6 specifies the number of significant decimals. Obviously the computing times would be significantly reduced if all algorithms were implemented in C or C++, but the findings, analysis, and conclusions would not differ. MILPs are solved using ILOG CPLEX 12.6.3 with a time limit of 3600 seconds. We also implemented the exact PLB Algorithm 1, by using GAMS 24.9.2 r64480 with LindoGlobal on neos servers (single core) [START_REF] Czyzyk | The neos server[END_REF] to solve the NLP subproblems that Algorithm 2 generates.

The computational evaluation was done in three parts: first, a comparison between continuous pwl δ-approximation and non necessarily continuous (nnc) pwl δ-approximation was done, for various nonlinear continuous functions. Then, the full PLB+MILP solution method was compared to MINLP solvers as well as state-of-the-art upper and lower bounding methods for solving the MINLP problem of energy optimization in HEV. Finally, the impact of the tolerance type (relative versus absolute) was evaluated.

Continuous pwl δ-approximation versus nnc pwl δ-approximation of nonlinear continuous functions

Let us consider the nine nonlinear continuous expressions summarized in Table 1. For each function and four different values of δ ∈ {0.100, 0.050, 0.010, 0.005}, [START_REF] Rebennack | Continuous piecewise linear delta-approximations for univariate functions: Computing minimal breakpoint systems[END_REF] reported the number of breakpoints of their continuous pwl δ-approximators. Their numerical accuracy was set to 10 -5 , meaning that equations ( 12) could be violated by at most 10 -5 . Their models and algorithms were implemented in GAMS 23.6 and solved using LindoGlobal 23.6.5. Their computations were performed on an Intel(R) i7 using a single core

2018 (v3) Technical report LAAS-CNRS n • 16358 17/ 29 ref function X -X + (I)
x 2 -3.5 3.5

(II) ln(x) 1 32 (III) sin(x) 0 2π (IV) tanh(x) -5 5 (V) sin(x) x 1 12 (VI) 2x 2 + x 3 -2.5 2.5 (VII) e -x sin(x) -4 4 (VIII)
e -100(x-2) 2 0 3 (IX) 1.03e -100(x-1.2) 2 + e -100(x-2) 2 0 3

Figure 4: For function f (x) = 2x 2 +x 3 and interval [X -, X + ] = [-2.5, 2.6], optimal continuous upper bounding leads to 8 line-segments whereas optimal nnc continuous upper bounding leads to 7 line-segments when δ = 0.5 with 2.93 GHz and 12.0 GB RAM on a 64-bit Windows 7 operating system. A time limit of 1800 seconds was set for every MINLP.

Table 2 shows the comparison between the continuous pwl δ-approximation and our nnc pwl δ-approximation obtained by computing a 2δ-overestimator before applying a shift of -δ on it. For functions (I) to (V), easier to approximate by Rebennack and Kallrath's algorithms, only the two smallest (and thus more challenging) δ values are considered, whereas for functions (VI) to (IX), harder to approximate, all four δ values were used. The first column contains the reference of the continuous input function considered. The second column specifies the number p of the subintervals of convexity and concavity for the nonlinear univariate input functions. The third column corresponds to the absolute tolerance value δ. In the remainder of the table, n * is the optimal number of line-segments, n -(resp. n + ) is the lower (resp. upper) bound returned when n * could not be found and cpu is the total computing time. Note that [START_REF] Rebennack | Continuous piecewise linear delta-approximations for univariate functions: Computing minimal breakpoint systems[END_REF] reported only the order of magnitude of the computing time for instances that could be solved to optimality. In their scale few means "≥ 1 and ≤ 10". The instances without computing times are the ones for which the 1800 seconds time limit was reached during the solution of one of the MINLPs. Passmark cpu scores can be used to estimate the computing times of our algorithms if they were run on a machine equivalent to the one of [START_REF] Rebennack | Continuous piecewise linear delta-approximations for univariate functions: Computing minimal breakpoint systems[END_REF]. Such conversion is available in Appendix. The computing times reported in Table 2 are the original values.

Results show that nnc heuristic and exact approximation methods require short computing times for all instances, contrary to continuous approximation which fails to converge for half of the instances. Exact nnc approximation yields a smaller number of pieces than continuous approximation 17 times out of 26, and the same number of pieces 9 times. The

(v3)

Technical report LAAS-CNRS n • 16358 19/ 29 heuristic nnc approximation yields an upper bound better or equal to the one from continuous optimization only 11 times out of 26, but its lower bound outperforms or equals the one from continuous approximation 19 times out of 26. In particular, on instances for which continuous approximation could not converge, the lower bound from the heuristic nnc improves the one from [START_REF] Rebennack | Continuous piecewise linear delta-approximations for univariate functions: Computing minimal breakpoint systems[END_REF] 10 times out of 12. Figure 4 illustrates the optimal continuous and nnc continuous pwl overestimators for function f (x) = 2x 2 + x 3 and interval [X -, X + ] = [-2.5, 2.6] with δ = 0.5. We can observe that the discontinuity helps to lessen the number of line-segments. In summary, the great advantage of using nnc pwl functions is the drastic reduction of the computing times, with the added benefit of sometimes also reducing the number of pieces.

7.2 Evaluation of the PLB+MILP solution method

Instances and PLB+MILP settings

The solution method PLB+MILP is evaluated on the real-world problem (P) of energy optimization for hybrid electric vehicles modeled with (1)-( 5). Instances derive from the data of [START_REF] Ngueveu | Heuristics and lower bound for energy management in hybrid-electric vehicles[END_REF] representing different driving cycles. The vehicle characteristics are:

P 1 min = 1 kW, P 1 max = 60 kW, P 2 min =-60 kW, P 2 max =60kW , E min =400kWs, E max =1600 kWs, E 0 =900 kWs, f 2 (x) = 1.0753x, f 3 (x) = 0.93x and f 1 (x) = 0 if x ≤ 0 f (R) (x) otherwise where f (R) (x) = 0.0000002x 5 -0.0000274x 4 + 0.00151450x 3 -0.02453270x 2 + 1.92434870x + 5.90568630. Function f (R) (x) is convex on [1; 7
.04029] and concave on [7.04029; 60]. The six different power demand profiles that define P i , ∀i ∈ {1...I} varied in size from I = 40 to I = 14001 . This results into 6 instances regrouped in a class denoted (R). Two other classes of instances denoted (A1) and (A2) were obtained, by replacing f (R) (x) with f (A1) (x) = 0.001x 3 -0.024x 2 + 1.92x + 5.91 or f (A2) (x) = -0.005x 3 + 0.5x 2 -0.8x + 10.0. Function f (A1) (x) results from the truncation of function f (R) (x) to avoid numerical errors in MINLP solvers. It is convex on [1; 8.0] and concave on [8.0; 60]. Function f (A2) (x) is an artificial function that is concave on [1; 33.3333] and convex on [33.3333; 60].

Applying PLB+MILP on problem (P) consists in identifying, for a given value of , two pwl functions f 1 and f 1 verifying equations ( 27)-( 29) to replace f 1 . Two MILPs are obtained: solving (P ) yields a lower bound for (P) whereas solving (P ) yields a feasible solution and an upper bound for (P). A tighter upper bound can be obtained by recomputing the cost of the feasible solution from (P ), using f 1 instead of f 1 . Three different values were tested: 0.01, 0.001, and 0.0001, corresponding to an expected gap of 1%, 0.1% and 0.01% between upper and lower bounds. The time limit was set to 3600 seconds for each MILP.

MINLP solvers and previous state-of-the-art upper and lower bounding methods

To the best of our knowledge the best known solution method for problem (P) was proposed by [START_REF] Gaoua | Energy management for an electric vehicle based on combinatorial modeling[END_REF]. It consists in a reformulation of the problem as a MILP using a discrete set of efficiency points obtained experimentally for the FC. This approach was replicated by [START_REF] Chauvin | Combinatorial approach for sizing and optimal energy management of hev including durability constraints[END_REF] in an iterative procedure based on the discretization of the nonlinear continuous efficiency function into equidistant efficiency points before the solution of the resulting sub-problems with a MILP solver. Applying such method on problem (P) consists in defining a set K FC of efficiency points k denoted ( x k , f 1 ( x k )) with x k ∈ [P 1 min , P 1 max ], defining binary variables y k i equal to 1 iff efficiency point k is applied at time i, replacing x 1 i with k∈K FC y k i in constraints (2) and replacing each nonlinear term f 1 (x 1 i ) of the objective function with the linear term k∈K FC f 1 ( x k )y k i , to obtain a MILP. The resulting solution method is denoted EP+MILP. To ensure a fair comparison between EP+MILP and PLB+MILP, the same number of binary variables is imposed by setting |K FC | = n , and the time limit is set to 3600 seconds.

The only known lower bound for problem (P) was proposed by [START_REF] Ngueveu | Heuristics and lower bound for energy management in hybrid-electric vehicles[END_REF], based on an assumption of ideal conditions: that the FC efficiency remains at its maximum level when it is used, i.e. each term f 1 (x 1 i ) is replaced with γx 1 i where γ = (max y∈[P 1 min ,P 1 max ]

f 1 (y) y ). The resulting relaxation is a MILP whose optimal solution cost is a lower bound of (P).

Problem (P) is a MINLP and can therefore be solved with any commercial or opensource MINLP solver. We used GAMS 24.9.2 r64480 and neos servers [START_REF] Czyzyk | The neos server[END_REF] to test four state-of-the-art solvers: Antigone, Baron, Couenne, and LindoGlobal. The time limit is set to 7200 seconds per instance, and the relative gap tolerance is set to , i.e. for = 0.01, GAMS parameter OPTCR is also set to 0.01. [START_REF] Gaoua | Energy management for an electric vehicle based on combinatorial modeling[END_REF]. (refer to Table 4 for MINLP solvers results)

Results and analysis

Tables 3,4, 5, and 6 report the computational results in an aggregated format where the first column is the instance class, the second column is the number of instances in the class, the third column specifies the relative tolerance , and the remainder of the table uses the following headings:

• n , n : numbers of pieces of the pwl functions f 1 and f 1 where UB is the upper bound obtained and LB+ is the best known lower bound from any of the lower bounding methods tested (therefore UB gap = 0% if UB=LB+). If the solver was not able to produce a solution within the time limit, the number of non-solved instances are reported in parenthesis. The average UB gap values are only computed over the instances for which a solution was produced.

• LB ratio: lower bound ratio, equal to 100 • LB UB+ where LB is the lower bound obtained, UB+ is the best known upper bound from any of the upper bounding methods tested (therefore LB ratio = 100% if LB=UB+). If the solver was not able to produce a valid lower bound within the time limit, the number of non-solved instances are reported in parenthesis. The average LB gap values are only computed over the instances for which a solution was produced.

(v3)

Technical report LAAS-CNRS n • Recomp: cost of the optimal solution of (P ) when recomputed using f 1 instead of f 1 .

Results from Table 3 show that PLB+MILP produces better solutions than EP+MILP in a shorter computing time. Table 4 shows that none of the MINLP solvers is able to complete its computations before the 7200 seconds time limit. With this time limit, Antigone and Baron produce solutions of worse quality than the ones from PLB+MILP. Couenne fails to produce a feasible solution in 32 cases out of the 54. LindoGlobal may produce very good solutions, but, due to its problem size limitations (3000 variables and 2000 constraints), the solver can only be applied on the 9 smallest cases out of 54. On these few cases, the solver produces better solutions than PLB+MILP in 6 cases, with a UB gap was improved by a value between 0% and 0.52%. In the 3 cases where PLB+MILP produces better solutions than the solver, the UB gap was worse by a value between 0.11% and 0.36%.

Regarding the lower bounds. Results from Table 5 show that PLB+MILP require a higher computing time, but produces better lower bounds, than the previous lower bounding procedure from the literature. Table 6 show that MINLP solvers are clearly outmatched. Baron, Couenne, and LindoGlobal fail to return a valid lower bound for several instances, and the few lower bounds returned are poor. Antigone always produces valid lower bounds but of poor quality, with a ratio varying between 21 and 54%. These results also suggest that all the MINLP solvers tested would require a prohibitively high computing time to run to completion.

In summary, PLB+MILP outperformed the best known upper and lower bounding methods from the literature, as well as MINLP solvers. New best known solution values are available in the column UB+ of Tables 7 and8.

Comparison between relative and absolute tolerance for piecewise linear bounding

This section focuses on the impact of the type of tolerance (relative or absolute) used for the first phase of the PLB+MILP solution method. For each pair "instance, ", given the best known solution cost . The value obtained is used to compute the pwl over-estimator f 1 δ and the resulting number of pieces n δ is displayed in Tables 7 and8. Results show that the δ value varies significantly from one instance to another: for example, for = 0.01 the δ values vary between 0.0322 and 0.2538 on instances of class R; thus the number of pieces vary, in this case between 10 and 25. Overall, on average the absolute δ tolerance led to more than twice the number of line-segments than the relative tolerance. PLB+MILP therefore outperforms δPLB+MILP on problem (P), by generating tighter MILPs requiring 20% to 80% less binary variables. Similar results were obtained with instances of class A1 and A2.

Conclusion

The paper presents a two-phase solution method for combinatorial optimization problems involving nonlinear univariate functions such as energy conversion functions. The first phase consists in bounding the nonlinear univariate functions from above and below with two piecewise linear non necessarily continuous functions with a relative tolerance. The second phase consists in solving the two mixed integer linear programs obtained when replacing the nonlinear terms with their piecewise linear overestimators or underestimators. Models and algorithms to perform the piecewise linear bounding are presented. Computational results on an energy optimization problem for hybrid electric vehicles illustrate the efficiency of the

(v3)

Technical report LAAS-CNRS n Passmark cpu score and equivalent computing times

The computing times of approximation algorithms reported in Table 2 were obtained with different computers with different characteristics, therefore Table 10 reports the equivalent cpu times if all procedures were run on a computer with a similar passmark score as the one from [START_REF] Rebennack | Continuous piecewise linear delta-approximations for univariate functions: Computing minimal breakpoint systems[END_REF]. Rebennack and Kallrath computed their continuous pwl approximators using a single core Intel(R) i7 with 2.93 GHz that scores a 5366 passmark cpu score. The heuristic pwl bounding procedure (Algorithm 3) was run on an Intel(R) Xeon(R) CPU E3-1271 v3 computer with 32GB RAM using a single score and having a passmark cpu score of 10086. The exact pwl bounding procedure (Algorithm 1) was run using GAMS 24.9.2 r64480 with LindoGlobal on neos servers (single core) [START_REF] Czyzyk | The neos server[END_REF] The first column of Table 10 contains the reference of the continuous input function considered. The second column reports the absolute tolerance value δ. In the remainder of the table, cpu is the original computing time, eq-cpu is the equivalent cpu time if the procedure was run on a computer with a similar passmark score as the one from [START_REF] Rebennack | Continuous piecewise linear delta-approximations for univariate functions: Computing minimal breakpoint systems[END_REF], # per server reports the number of NLP sub-problems solved on each type of neos server, and |N I | is the value of grid size reached when computing each valid line-segment (Algorithm 2).

(v3)

Technical report LAAS-CNRS n

  Figure 2: Maximization of the first interval (projection on D of the first line-segment)

:

  Function f ; limits x begin , x end ; scalar δ > 0; parameters n init I ≥ 2, ω > 1 Output: slope a, y-intercept b, next breakpoint y 1: initialize gap := +∞; grid size n I := n init I /ω and grid N I = {1 . . . n I } 2: while gap > δ do 3: update grid size n I := n I ω and grid N I := {1 . . . n I } 4: solve the following NLP to obtain y, s 0 , and s y in order to compute the coefficients a := f (y)+sy-x begin -s 0 y-x begin and b := f (x begin ) + s 0 -ax begin : max y //maximize ending breakpoint (16)

BB

  max = B min + 0.5(B max -B min ) , b, B min , B max }=Compute_q(f ,x begin , X + ,δ,α,B min , B max , q, a, b) 12: end if Algorithm 7 Dichotomy_y Input: Convex function f ; X + ; δ > 0; α ∈ N; q; a; b; B min , B max Output: next breakpoint y, lower limit B min , upper limitB max 1: if B max -B min ≤ 10 -α then f (B max ) -(aB max + b) ≤ δ then 6:B min = B max ; B max = X + max = B min + 0.5(B max -B min )

  29) Definition 6.1 generalizes the one of Ngueveu et al. (2016) by taking into account D ⊂ R instead of D ⊂ R + .

Table 1 :

 1 Univariate nonlinear functions tested

	input		continuous approximation		nnc approximation
	function	δ	Rebennack and Kallrath (2015)	Heuristic		Exact
	ref p		n cpu
	(I) 1 0.010 25	Hours 25		30 s 25	19 s
		0.005 35	Hours 35		43 s 35	4 s
	(II) 1 0.010	9	Sec	9		11 s	9	221 s
		0.005 13	Sec 13		16 s 13	172 s
	(III) 2 0.010 13	Sec	13 14 12 s 13	57 s
		0.005 17	Few min	17 18 17 s 17	68 s
	(IV) 2 0.010	9	Few sec	9 10 10 s	9	161 s
		0.005 13	Few min	13 14 15 s 13	128 s
	(V ) 4 0.010	9	Sec	7 10 11 s	8	143 s
		0.005 12	Few min	12 15 15 s 12	181 s
	(VI) 2 0.100 12	Min	11 12 11 s 11	115 s
		0.050 16	Few days	15 16 17 s 15	88 s
		0.010	16 35	34 35 36 s 34	164 s
		0.005	16 48	47 48 59 s 47	195 s
	(VII) 3 0.100	5 15	14 16 17 s 14	226 s
		0.050	5 20	19 21 28 s 19	287 s
		0.010	5 44	43 45 52 s 43	268 s
		0.005	5 62	61 63 72 s 61	869 s
	(VIII) 3 0.100	5	Sec	4	6	4 s	4	74 s
		0.050	5	7	4	6	6 s	5	83 s
		0.010	5 12	10 12 11 s 11	138 s
		0.005	5 15	14 16 16 s 14 1466 s
	(IX) 5 0.100	8	Few days	7 11	8 s	7	77 s
		0.050	8 12	7 11 12 s	9	64 s
		0.010	8 22	19 23 23 s 21	114 s
		0.005	8 29	27 31 27 s 27	784 s

* n -n + cpu n * n -n + cpu n *

Table 2 :

 2

Comparison between nnc δ-approximation and continuous δ-approximation for nine univariate nonlinear continuous functions previously used by

[START_REF] Rebennack | Continuous piecewise linear delta-approximations for univariate functions: Computing minimal breakpoint systems[END_REF]

2018 (v3)

Technical report LAAS-CNRS n • 16358 18/ 29 (a)

Table 3 :

 3 Comparison between PLB+MILP and the best known state-of-the-art upper bounding method EP+MILP from

Table 4 :

 4 Evaluation of MINLP solvers upper bounds for problem (P), for a computation time of 7200s. Results show that these solvers are outperformed by the algorithms of Table3

	Instances	LB from Ngueveu et al. (2017)		PLB+MILP PLB (lower)	MILP
	class	# LB ratio	cpu		n	cpu LB ratio	cpu
	(R)	6 97.81 %	2 s	0.01	6	8 s	99.18 %	20 s
					0.01	19	19 s	99.93 %	31 s
					0.001	56	50 s 99.99 %	58 s
	(A1)	6 97.36 %	2 s	0.01	10	11 s	99.25 %	29 s
					0.001	27	27 s	99.94 %	46 s
					0.0001	82	74 s 99.99 % 216 s
	(A2)	6 78.50 %	2 s	0.01	14	16 s	99.50 %	26 s
					0.001	43	42 s	99.95 %	87 s
					0.0001 134	114 s 99.99 % 405 s

Table 5 :

 5 Comparison between PLB+MILP and the best known state-of-the-art lower bounding method for problem (P) (from[START_REF] Ngueveu | Heuristics and lower bound for energy management in hybrid-electric vehicles[END_REF]). (refer to Table6for MINLP solvers results)

	• UB gap: upper bound gap, equal to 100• UB-LB+ LB+

Table 6 :

 6 Evaluation of MINLP solvers lower bounds for problem (P), for a computation time of 7200s. Results show that these solvers are outperformed by the algorithms of Table5

  UB+, we evaluate what would have been the value of absolute tolerance δ

		Class (R)		using relative tol	using absolute tol	Gap
	instance	I	UB+		n	UB δ = U B• I	n δ 100 n δ -n n
	S_(R)	40	453.34	0.01	6	454.3	0.1133	14 133.33 %
				0.001 19	453.5	0.0113	42 121.05 %
				0.0001 56	453.4	0.0011	133 137.50 %
	I_(R)	561	8740.92	0.01	6	8756.6	0.1558	12 100.00 %
				0.001 19	8742.2	0.0156	36 89.47 %
				0.0001 56	8741.0	0.0016	110 96.43 %
	H_(R)	734 18568.41	0.01	6 18626.0	0.253	10 66.67 %
				0.001 19 18576.9	0.0253	28 47.37 %
				0.0001 56 18569.0	0.0025	89 58.93 %
	U_(R)	811	2607.60	0.01	6	2613.5	0.0322	25 316.67 %
				0.001 19	2608.4	0.0032	78 310.53 %
				0.0001 56	2607.7	0.0003	255 355.36 %
	N_(R) 1200 23114.87	0.01	6 23137.7	0.1926	11 83.33 %
				0.001 19 23119.3	0.0193	33 73.68 %
				0.0001 56 23114.9	0.0019	102 82.14 %
	E_(R) 1400 27065.71	0.01	6 27088.9	0.1933	11 83.33 %
				0.001 19 27075.3	0.0193	33 73.68 %
				0.0001 56 27065.9	0.0019	102 82.14 %
	2018 (v3)		Technical report LAAS-CNRS n • 16358		23/ 29

Table 7 :

 7 Comparison between relative pwl bounding and absolute pwl bounding for instances of class (R) to be used for pwl bounding with absolute tolerance in order to ensure an similar solution cost.

It is computed as: δ = •U B+ I

Table 8 :

 8 Comparison between relative pwl bounding and absolute pwl bounding for instances of classes (A1) and (A2) solution method. Results also show that using non necessarily continuous piecewise linear functions to approximate continuous univariate nonlinear functions leads to a significant reduction of the computing time and a reduction of the number of line-segments in comparison to classical continuous piecewise linear functions.

	2018 (v3)	Technical report LAAS-CNRS n • 16358	25/ 29

  to solve the NLP sub-problems that Algorithm 2 generates. It is not possible to chose a priori a specific neos machine, but it is possible to know a posteriori on which machine a code submitted had been run. Five types of servers (P1, P2, P3, P4 and P5) were available on the neos website, with the following specs. P1 is a Dell PowerEdge R430, Intel Xeon E5-2698 with 2.3GHz, that has a passmark score of 21149. P2, P3 and P4 are Dell PowerEdge R410, Intel Xeon X5660 with 2.8GHz that have a passmark score of 7641. P5 is a Dell PowerEdge R420, Intel Xeon E5-2430 with 2.2GHz that have a passmark score of 6878.

Power profiles derived from real-wold test drive cycles[START_REF] Ngueveu | Heuristics and lower bound for energy management in hybrid-electric vehicles[END_REF]) are available at http://homepages.laas.fr/sungueve/Data/HEV_I_PowerProfiles.zip.

(v3) Technical report LAAS-CNRS n • 16358 20/ 29

(v3) Technical report LAAS-CNRS n • 16358 29/ 29

Acknowledgements

This research benefited from the support of the FMJH Program PGMO, from the support of EDF-Thales-Orange and from the support of the PEPS program from the Cellule Energie of the CNRS.

The author would also like to thank the two anonymous reviewers for their useful comments that helped improving the quality of the paper.

Appendices

Characteristics of the mathematical formulations