Aloïs Duguet

Sandra Ulrich Ngueveu

Piecewise linearization of bivariate nonlinear functions: minimizing the number of pieces under a bounded approximation error

published or not. The documents may come

Piecewise linearization of bivariate nonlinear functions: minimizing the number of pieces under a bounded approximation error

Aloïs Duguet, Sandra Ulrich

Problem Description and State of the Art

Let (P) be the optimization problem of approximating a nonlinear function f of two variables by a piecewise linear (PWL) function g subject to approximation error constraints on domain D represented by functions l and u satisfying l(x, y) ≤ f (x, y) ≤ u(x, y) for all (x, y) ∈ D:

(P)      min n (1)
subject to l(x, y) ≤ g(x, y) ≤ u(x, y) ∀(x, y)

∈ D ⊂ R 2 (2)
g is a PWL function with n pieces (3) Constraints (2) are pointwise approximation constraints, which gives an infinite number of constraints because D is a continuous domain. Piecewise linear approximation are more commonly minimizing an approximation error with a fixed number of pieces [START_REF] Toriello | Fitting piecewise linear continuous functions[END_REF]), but our objective is to minimize the number of pieces of g so that an MILP formulation of g introduces less binary variables [START_REF] Vielma | Mixed-integer models for nonseparable piecewise-linear optimization: Unifying framework and extensions[END_REF]). It can be useful for the approximation of a Mixed Integer Nonlinear Programming problem (MINLP) with a Mixed Integer Linear Programming (MILP) by substituting each nonlinear functions by a PWL one. Moreover, it was shown in [START_REF] Geißler | Mixed Integer Nonlinear Programming, chapter Using piecewise linear functions for solving MINLPs[END_REF]) that in some cases, MINLP can be solved by applying only techniques from MILP.

Heuristics and exact methods exist to approximate univariate nonlinear functions (D ⊂ R) [START_REF] Codsi | Lina: A faster approach to piecewise linear approximations using corridors and its application to mixed-integer optimization[END_REF], [START_REF] Ngueveu | Piecewise linear bounding of univariate nonlinear functions and resulting mixed integer linear programming-based solution methods[END_REF], [START_REF] Rebennack | Piecewise linear function fitting via mixed-integer linear programming[END_REF]). We are interested in bivariate functions (D ⊂ R 2 ) and to the best of our knowledge, only two papers address this case, with heuristics only:

• The authors of [START_REF] Rebennack | Continuous piecewise linear delta-approximations for bivariate and multivariate functions[END_REF]) propose two heuristics to solve problem (P) with continuous PWL functions. The first heuristic is based on an iterative subdivision of the domain D into triangles (2-simplexes) until for each subdomain a linear function satisfying the approximation error has been found. The verification that a given linear function fits the subdomain is made by solving a nonlinear programming problem (NLP). The second heuristic can be used if the contribution of the two variables in the function can be separated in two univariate functions (linearly or nonlinearly). In this case, an algorithm finds the two optimal continuous univariate PWL functions and combine them to build a single two-variable PWL function.

• In [START_REF] Kazda | Nonconvex multivariate piecewise-linear fitting using the difference-of-convex representation[END_REF]) an iterative process attempts to find a continuous PWL function written as a Difference of Convex Continuous PWL functions (DC CPWL) that satisfies the approximation error. The idea is to iteratively solve an MILP relaxation of (P) and then to find lazy constraints to add to the relaxation until a solution found is feasible for (P).

The relaxation consists in replacing the infinite number of constraints (2) with a finite number of them.

After the introduction of definitions used throughout the paper in Section 2, the three key ideas of a framework for piecewise linearization are detailed in Section 3. It is followed by explanations on the instantiation of crucial parts of this framework to create different heuristics in Section 4, and finally, numerical experiments comparing the state of the art to our best heuristics are shown in Section 5.

Definitions

The vocabulary used throughout this work is presented below. They are in part extensions to D ⊂ R 2 of definitions from [START_REF] Codsi | Lina: A faster approach to piecewise linear approximations using corridors and its application to mixed-integer optimization[END_REF]) for D ⊂ R.

Definition 1 (Polytope). A polytope P is the convex hull of some points

X i ∈ R m : P = {x ∈ R m , x = i λ i X i , i λ i = 1, λ i ≥ 0 ∀i}.
Throughout the paper, a piece will refer to a polytope that composes the graph of a PWL function, and 1, n := {1, ..., n}.

Definition 2 (PWL function). Let D be a compact set of R 2 . A function g : D → R is a PWL function with n pieces if and only if there exists

{a i } i∈ 1,n ⊂ R 2 , {b i } i∈ 1,n ⊂ R and a family of polytopes {D i } i∈ 1,n ⊂ R 2 such that D = ∪ i∈ 1,n D i ,
and for i ̸ = j the polytopes D i and D j can only intersect on their boundary and g is defined by g(x, y)

= min{a i • (x, y) T + b i |(x, y) ∈ D i ∀i ∈
1, n }, with • denoting the standard scalar product.

Precautions were taken to allow g to be not necessarily continuous at the boundary of a polytope D i . To prevent g(x, y) to have multiple definitions because {D i } i∈ 1,n can intersect on their boundary, g(x, y) is chosen as the minimum over all possible definitions, so that it is a lower semicontinuous function. Definition 8 (Inner corridor). Let C 0 be a corridor between u 0 and l 0 . Let C be a corridor between u and l. We call C an inner corridor of C 0 if and only if D(C) = D(C 0 ) and l 0 (x, y) ≤ l(x, y) < u(x, y) ≤ u 0 (x, y). Definition 9 (R 2 -corridor fitting problem). The R 2 -corridor fitting problem consists in finding a PWL function g of two variables fitting a corridor C such that its number of pieces is minimized.

Definition 3 (Corridor). Let D be a compact set of R 2 . Let u, l : D → R be two continuous functions verifying u(x, y) > l(x, y), ∀(x, y) ∈ D. Define the set C = {(x, y, z) ∈ R 3 |(x, y) ∈ D, l(x, y) ≤ z ≤ u(x,
(P) is equivalent to an R 2 -corridor fitting problem with corridor C between u and l. Moreover, by extending the previous definitions to dimension m ≥ 1, it is possible to define the R m -corridor fitting problem that refers to the problem with

D(C) ⊂ R m . Definition 10 (Truncated-corridor). Let C 1 , C 2 ∈ R 2 be two corridors both defined by functions u, l : R → R on the intervals [a, b 1 ] and [a, b 2 ]. We call C 2 a truncated-corridor of C 1 if and only if [a, b 2 ] ⊂ [a, b 1 ] (b 2 ≤ b 1 ).
Definition 11 (Maximal linear segment). A maximal linear segment in a corridor C is a linear segment within C that induces a truncated-corridor of maximal domain length.

Definition 12 (Truncated-corridor in direction d). Let C 1 and C 2 be two corridors defined by the same functions u and l with compact corridor domains in

R 2 . Let d ∈ R 2 \ {0}. We call C 2 a truncated-corridor of C 1 in direction d if and only if there exists σ ∈ R for which D(C 2 ) = D(C 1 ) ∩ {(x, y) ∈ R 2 , (x, y) • d ≤ σ}, i.e. D(C 2 ) is the intersection of D(C 1 ) with a half-plane.
Definition 13 (Maximal piece in direction d). A maximal piece in direction d ∈ R 2 \ {0} of a corridor C is a polytope within C that induces a truncatedcorridor of C in direction d that is of maximal domain area.

A Framework for Solving the R 2 -Corridor Fitting Problem

We present in this section a framework to create efficient algorithms for the R 2 -corridor fitting problem. The instantiation chosen for different parts of the framework are described in Section 4 as well as some details on the implementation.

Three key ideas are followed in the framework. The first two are used to avoid drawbacks encountered in [START_REF] Rebennack | Continuous piecewise linear delta-approximations for bivariate and multivariate functions[END_REF]), whereas the third one is meant to render a subproblem more tractable:

• Key idea 1: Pieces should be chosen among general convex polygons instead of triangles to constrain less the pieces chosen, possibly decreasing the number of pieces

• Key idea 2: Choose pieces that are good (ideally optimal) solutions of a maximal piece in direction d problem, to aim for a domain covered by a piece that is "as large as possible" because in [START_REF] Rebennack | Continuous piecewise linear delta-approximations for bivariate and multivariate functions[END_REF]) the size of triangles is fixed which increases the number of pieces necessary

• Key idea 3: Compute a good feasible solution of the maximal piece in direction d problem with a series of LP problems obtained after substituting C with a PWL inner corridor of C

The remainder of this section builds upon these principles.

Key Idea 1: Management of the Corridor Domain

The corridor domain should be tiled with shapes as general as possible provided that they can be formulated in an MILP. Such shapes are polygons, but we further restrict those shapes to convex polygons because formulating a non-convex polygon in an MILP introduces additional binary variables. It is expected that allowing convex polygons instead of only triangles as done in [START_REF] Rebennack | Continuous piecewise linear delta-approximations for bivariate and multivariate functions[END_REF]) will lead to a lower number of pieces.

The procedure that manages pieces and the remaining corridor domain is described in Algorithm 1. At each iteration, one piece is computed for each vertex of D(C) by function compute_piece, and a function score selects the "most" suitable to obtain a PWL function with few pieces, see Section 4.1. Function update_domain(C, p) removes the part of D(C) on which p is defined. This function also divides the new reduced corridor C in two corridors if polygon D(C) only has angles at the vertices greater than 90°to avoid a bad behaviour.

Algorithm 1 Finding a PWL function fitting a corridor C with a low number of pieces

1: function pwl_2d_fitting(C) 2: P ← ∅ ▷ list of chosen pieces 3: Q ← {C} ▷ list of corridors with convex domains not already tiled 4: while Q ̸ = ∅ do 5: C = pop(Q) 6: candidate_pieces ← ∅ 7: for v vertex of D(C) do 8: d ← choose_progress_direction(C, v) 9: candidate_pieces ← candidate_pieces ∪ {compute_piece(C, v, d)}
(M P d )    max σ s.t. αx + βy + γ ∈ C d σ ∀(x, y) ∈ D(C d σ ) α, β, γ, σ ∈ R (4) Where D(C d σ ) = D(C) ∩ {(x, y) ∈ R 2 |(x, y) T • d ≤ σ} is the domain of C d σ , the truncated-corridor of C in direction d. (M P d
) is a generalized semi-infinite programming problem (GSIP). Indeed, the number of pointwise constraints is infinite and depends on variable σ, while the number of variables is 4: a real variable σ for the half-plane intersection as well as 3 real variables (α, β, γ) to describe the linear function g(x, y) = αx + βy + γ.

Key

(M P ′ d )          max σ s.t. αx + βy + γ -ũi(x, y) ≤ 0 ∀(x, y) ∈ D(C d σ ) ∩ D i ũ, ∀i ∈ I αx + βy + γ -lj(x, y) ≥ 0 ∀(x, y) ∈ D(C d σ ) ∩ D j l , ∀j ∈ J α, β, γ, σ ∈ R (5) Note that on each D i ũ, g(x, y) -ũ(x, y) is a linear function, thus it suffices to check g(X) -ũ(x, y) ≤ 0 for each vertex of convex polygonal domain D i ũ to ensure constraint g(x, y) -ũ(x, y) ≤ 0 on D i ũ.
A similar reasoning leads to the same result for constraints involving l. Thus, the following problem is equivalent to (M P ′ d ) but has the advantage of having only a finite number of linear constraints.

(M P ′′ d )          max σ s.t. αx + βy + γ -ũi(x, y) ≤ 0 for each vertex (x, y) of D(C d σ ) ∩ D i ũ, ∀i ∈ I αx + βy + γ -lj(x, y) ≥ 0 for each vertex (x, y) of D(C d σ ) ∩ D j l , ∀j ∈ J α, β, γ, σ ∈ R (6)
Finally, problem (M P ′′ d ) has constraints involving polygon intersections depending nonlinearly on variable σ, thus it is not an LP problem. Parameterizing (M P ′′ d ) with σ, the following LP feasibility problem is obtained which can be repeatedly solved until a satisfactory σ value has been found.

(M P ′′ d,σ )    αx + βy + γ -ũi(x, y) ≤ 0 for each vertex (x, y) of D(C d σ ) ∩ D i ũ, ∀i ∈ I αx + βy + γ -lj(x, y) ≥ 0 for each vertex (x, y) of D(C d σ ) ∩ D j l , ∀j ∈ J α, β, γ ∈ R (7)

Framework Key Points Instantiation

In this section, choices made on key points of the framework of Section 3 are described: the function score of Algorithm 1 in Section 4.1, the choice of a direction d in Section 4.2 and the computation of a PWL inner corridor in Section 4.3. C refer to a corridor between u and l in the remainder of the section.

Scoring the Quality of Pieces

In Algorithm 1, a function score ranks the quality of candidate pieces and the piece with highest score is kept. Two scoring functions have been implemented:

• Area measures the area of the domain covered by piece p

• Partial Derivatives Total Variation (PaD) approximates the sum of the total variation of each partial derivative of u and l on the domain covered by p

The total variation of a function g is a measure of how much that function varies on its domain D. The total variation of g on D is equal to D ||∇g(x)|| 2 dx. The interest of P aD needs the introduction of the pointwise height of a corridor.

Definition 14 (pointwise We call C P H (x, y) = u(x, y) -l(x, y) the pointwise height of corridor C at point (x, y).

Remark 15. The most commonly used type of approximation error for a function f is the absolute error. It induces a corridor C such that l(x, y) = f (x, y)-δ and u(x, y) = f (x, y) + δ with δ > 0, that has constant pointwise height.

Area is a straightforward and simple idea to evaluate the piece quality, it will serve as a reference to evaluate other scoring functions. PaD is thought to be an adaptation for two-variable functions of Theorem 1 of [START_REF] Frenzen | On the number of segments needed in a piecewise linear approximation[END_REF]). Indeed, it states that for a corridor C of constant pointwise height 2δ with D(C) = [a, b], when δ → 0, the minimum number of pieces s(δ) satisfies the asymptotic approximation s(δ) ∼ 1 4δ b a |u ′′ (x, y)|dx. Note first that u ′′ (x, y) = l ′′ (x, y) because it is a corridor with constant pointwise height, and second that the integral computed is the total variation of function u ′ (and l ′ ). It is thus expected that PaD performs better than Area for small values of pointwise height. For large values, the total variation should be less relevant to estimate the difficulty of fitting a piece, thus PaD could be less efficient.

Choose a Progress Direction

Line 8 of Algorithm 1 selects progress direction d knowing starting vertex v. Two options were tested for this choice. The first is a naive option, while the second is meant to take into account the "difficulty" of progressing along the two extremal directions given by the two edges starting at v.

Inner Approximation of a Corridor

In the function compute_piece, the computation of a PWL inner corridor C PWL of C boils down to computing PWL functions ũ and l verifying l(x, y) ≤ l(x, y) ≤ ũ(x, y) ≤ u(x, y), for all (x, y) ∈ D(C).

To compute ũ (a similar method works for l), a basic idea is to divide D(C) into rectangular pieces of same sizes, and then to compute a third coordinate ũ(x, y) to each vertex v = (x, y) of each rectangular piece such that ũ(x, y) ≤ u(x, y). Interval analysis on the gradient of u suffices to compute the values of ũ(x, y) such that ũ is an underestimation of u, as explained in Proposition 16.

Proposition 16. Let u ∈ C 1 defined on D = [a, b] × [c, d] ⊂ R 2 . Let l x = b -a and l y = d -c. Let ∇u be the gradient of u. Let [D low x , D high x ] × [D low y , D high y ] be such that ∇u(x, y) ∈ [D low x , D high x ] × [D low y , D high y ] for all (x, y) ∈ D. Let (M x , M y ) = ( a+b 2 , c+d 2 ). Define: u - (a,c) := u(M x , M y ) -D high x • l x -D high y • l y (8) u - (b,c) := u(M x , M y ) + D low x • l x -D high y • l y (9) u - (b,d) := u(M x , M y ) + D low x • l x + D low y • l y (10) u - (a,d) := u(M x , M y ) -D high x • l x + D low y • l y (11)
If a linear function f satisfies:

f (a, c) ≤ u - (a,c) , f (b, c) ≤ u - (b,c) , f (b, d) ≤ u - (b,d) and f (a, d) ≤ u - (a,d) (12)
Then f (x, y) ≤ u(x, y) for all (x, y) ∈ D.

Proof. Let f be a linear function satisfying the 4 inequalities ( 12). Let M = (M x , M y , f (M x , M y )) be the point on the surface defined by u corresponding to the middle of D. Let (x 0 , y 0 ) ∈ D. We have:

u(M x , M y ) -D high x • (M x -x 0 ) -D high y • (M y -y 0 ) ≤ u(x 0 , y 0 ) if x 0 ≤ M x , y 0 ≤ M y u(M x , M y ) + D low x • (x 0 -M x ) -D high y • (M y -y 0 ) ≤ u(x 0 , y 0 ) if x 0 ≥ M x , y 0 ≤ M y u(M x , M y ) + D low x • (x 0 -M x ) + D low y • (y 0 -M y ) ≤ u(x 0 , y 0 ) if x 0 ≥ M x , y 0 ≥ M y u(M x , M y ) -D high x • (M x -x 0 ) + D low y • (y 0 -M y ) ≤ u(x 0 , y 0 ) if x 0 ≤ M x , y 0 ≥ M y because [D low x , D high x ] × [D low y , D high y
] are bounds of ∇u on domain D. Now, define a PWL function f PWL with four rectangle pieces with vertices positionned at {(a, c), (b, c), (b, d), (a, d), ( a+b 2 , c+d 2 )} and height the left-hand sides of ( 12) as well as u(M x , M y ) respectively. In particular, f PWL ≤ u on D.

In addition, direct computations show that f (x, y) ≤ f PWL (x, y) for (x, y) ∈ {(a, c), (b, c), (b, d), (a, d), ( a+b 2 , c+d 2 )}. Finally, as f and f PWL are linear on the four pieces domain of f PWL , we have f ≤ f PWL ≤ u on D.

To build a PWL inner corridor exploiting Proposition 16 in our algorithm, a method called efficiency refinement is used. It is described after the introduction of the bounding efficiency η.

Definition 17 (bounding efficiency η). Let C be a corridor with D(C) a polygonal domain of R 2 . Let C PWL be a PWL inner corridor of C. We say that C PWL achieves a bounding efficiency for C of η ∈ [0, 1] if the pointwise height (PH) ratio (C PWL ) P H (x,y)

C P H (x,y)
is greater or equal to η for each (x, y) ∈ D(C).

Proposition 18. Let C be a corridor with D(C) a polygonal domain of R 2 and let η ∈ [START_REF] Kazda | Nonconvex multivariate piecewise-linear fitting using the difference-of-convex representation[END_REF]1]. If C has constant pointwise height, then a PWL inner corridor C PWL has a bounding efficiency for C of η if the pointwise height ratio

(C PWL ) P H (x,y) C P H (x,y)
is greater or equal to η for each (x, y) vertex of a piece domain of ũ or l.

Proof. C has constant pointwise height. Thus for each piece of C PWL , the minimum pointwise height ratio is on an extreme point of the piece domain, that is to say on a vertex of the piece domain, which is lower bounded by η by hypothesis.

The efficiency refinement procedure builds a PWL inner corridor C PWL of C achieving a bounding efficiency of η if it has a constant pointwise height, but without this property, it only checks that the pointwise height ratio at the vertices of each piece is η.

After having found a rectangle D containing D(C), it creates an initial PWL corridor likely unvalid (l | ≤ u) with only one rectangular piece on D, and then iteratively refines the pieces that do not satisfy a bounding efficiency of η at the 4 vertices into 4 new pieces until each piece satisfies the efficiency, which implies the validity of the PWL inner corridor as well.

Parameter η needs to be adjusted depending on the quality of PWL inner corridor C PWL wanted. To produce a really good approximation of corridor C, η near 1 shall be used, but the number of pieces forming C PWL increases consequently, thus increasing the computation time of (M P ′′ d,β ) to be solved later on.

Numerical Experiments

The performance of our framework is compared to the state of the art. Our best heuristic (found with experiments in appendix [START_REF] Duguet | Appendix of "piecewise linearization of bivariate nonlinear functions: minimizing the number of pieces under a bounded approximation error[END_REF])) called DN99 uses scoring function PaD, starting vertex set to med and η = 0.99. Also, to illustrate the effects of parameter η, results of a heuristic which uses the x 2 -y 2 [0.5, 7.5] × [0.5, 3.5] N3 xsin(y)

[1.0, 4.0] × [0.05, 3.1] L2 x 2 + y 2 [0.5, 7.5] × [0.5, 3.5] N4 sin(x) x y 2 [1.0, 3.0] × [1.0, 2.0] N5 xsin(x)sin(y) [0.05, 3.1] × [0.05, 3.1] N1 xy [2.0, 8.0] × [2.0, 4.0] N6 (x 2 -y 2 ) 2 [1.0, 2.0] × [1.0, 2.0] N2 xe -x 2 -y 2 [0.5, 2.0] × [0.5, 2.0] N7 e -10(x 2 -y 2 ) 2 [1.0, 2.0] × [1.0, 2.0]
same parameters as DN99 but with η = 0.95, are shown. This second heuristic is called DN95. Heuristics from the state of the art are described in articles [START_REF] Codsi | Lina: A faster approach to piecewise linear approximations using corridors and its application to mixed-integer optimization[END_REF], [START_REF] Kazda | Nonconvex multivariate piecewise-linear fitting using the difference-of-convex representation[END_REF], [START_REF] Rebennack | Continuous piecewise linear delta-approximations for bivariate and multivariate functions[END_REF]). Recall that [START_REF] Rebennack | Continuous piecewise linear delta-approximations for bivariate and multivariate functions[END_REF]) proposes two heuristics: one based on triangulation (RK2D), and one based on a decomposition of the bivariate function into a sum of one variable functions to approximate separately (RK1D). Heuristic LinA2D is based on the same decomposition as RK1D, but uses library LinA from [START_REF] Codsi | Lina: A faster approach to piecewise linear approximations using corridors and its application to mixed-integer optimization[END_REF]) to compute the approximation of univariate functions. Obviously LinA2D outperforms RK1D since LinA computes optimal non necessarily continuous univariate PWL functions instead of optimal continuous univariate PWL functions for RK1D. The heuristic of [START_REF] Kazda | Nonconvex multivariate piecewise-linear fitting using the difference-of-convex representation[END_REF]) is denoted KL2D.

Those 5 heuristics are compared on the benchmark instances of [START_REF] Rebennack | Continuous piecewise linear delta-approximations for bivariate and multivariate functions[END_REF]), which consist of 45 instances obtained from 9 functions and 5 different absolute approximation errors for each function. An absolute approximation error of δ for function f means that the corresponding corridor C is between functions u and l satisfying l(x, y) = f (x, y) -δ and u(x, y) = f (x, y) + δ. The first two functions of the benchmark are linearly separable and refered to as L1,L2. The remaining seven are not linearly separable and thus refered to as N1,...,N7. The expression and domain of each function are described in Table 1.

For our heuristics as well as LinA2D, JuMP (v. 0.21.8) is used as modeling language, Gurobi (v. 9.1) is the (MI)LP solver and a CPU of 4.4 GHz using a single core and 32GB RAM. Moreover, heuristic RK2D uses the modeling language GAMS (v. 23.6), the global optimization solver LindoGlobal (v. 23.6.5), a CPU intel i7 with a single core, 2.93 GHz and 12 GB RAM. Finally, KL2D uses Pyomo (v. 5.6.4) as modeling language, CPLEX (v. 12.8) and 10 threads to solve MILPs, COUENNE (v.0.5.8) and 1 thread to solve NLPs, a CPU with 3.6 GHz and 32 GB RAM. As for the time limit, LinA2D and RK2D do not have one, KL2D allows 3600 seconds for each MILP problem, while our heuristics allow 3600 seconds for each LP problem. The differences in computation time between heuristics are in several orders of magnitude, therefore these differences are mainly due to the algorithms differences and are only marginally impacted by the differences in hardware.

Results are shown in Table 2. For each of the five heuristics we present the number of pieces n obtained by the heuristic as well as its computation time in seconds. Bold integer highlight the best solution found for each instance, i.e. the ones with the minimum number of pieces. "TO" means that the heuristic has stopped because of a time out, thus without any valid solution. In this case, "-" means that no solutions were found.

In terms of minimum number of pieces, DN99 performs the best with 28 best solutions out of 45 instances, followed by KL2D with 19 out of 45, DN95 with 15 out of 45, LinA2D with 8 out of 45 and finally RK2D with 1 out of 45. As for the computation time, the hardware and software used for the different heuristics are of different quality, thus it will not be a precise tool to compare the time needed for each heuristics. However, it can be said that LinA2D takes less than 1 second to execute, RK2D and DN95 take seconds to minutes, while KL2D and DN99 take seconds to hours.

A more in-depth analysis shows that LinA2D is the best only for instances with linearly separable functions (L1 and L2), but it does really poorly on the other functions. KL2D is the best 19 times out of the 24 instances where it terminates with a solution, which shows a clear limit to the size of the instance it can tackle, due to the increasing size of MILPs it solves. The effect of the change of parameter η from 95% to 99% in heuristics DN95 and DN99 is on average a decrease of 9.1% of the number of pieces of the solutions, at the cost of an increase of 301% of the computation time. It is an open question what value of parameter η should be taken in order to achieve a certain trade off between computation time and minimisation of the number of pieces.

Conclusion

We introduced a framework to create linearization algorithms based on the solution of the R 2 -corridor fitting problem. Convex polygons were used to tile the domain instead of triangles, a scoring function ranking candidate pieces was developed and a good feasible solution of a GSIP was computed via a series of LP feasibility problems. Finally, numerical experiments showed that our heuristics outperforms the state of the art on not linearly separable functions. Further work could attempt to diminish the relatively high computation time of the heuristics, search for better instanciation of the different parts of the framework or apply those heuristics to approximate real-world MINLPs to show their practical usage.

  y)} as the corridor between u and l. If D ⊂ R 2 , we call the area of D the domain area of C. A similar definition can be made for D an interval [a, b], in which case we call b -a the domain length of C. Definition 4 (Corridor domain). Let C be a corridor, C ⊂ R 3 . The domain of corridor C noted D(C) is the projection of C on its two first coordinates, which is also the domain on which u and l need to be defined. Definition 5 (Piece within a corridor). A polytope P ⊂ R 3 is within a corridor C if and only if there exists a linear function g : D ⊂ D(C), such that P = {(x, y, g(x, y)), (x, y) ∈ D} and P ⊂ C. Definition 6 (Fitting). A PWL function g fits a corridor C if and only if the pieces of g (polytopes {P } i∈ 1,n of the graph of g) are within C and g is defined on the entire domain D(C). Definition 7 (PWL corridor). A corridor C is called a PWL corridor if and only if u and l defining C are both PWL functions.

  Idea 2: the Maximal Piece in Direction d Problem We chose to find a new piece by covering an area starting from point v and extending as far as possible in direction d. This direction d points to the interior of D(C) when starting from v and is computed via function choose_progress_direction of Algorithm 1. The hypothesis of starting from a vertex instead of any point of the border of the polygon D(C) is made. Computing the piece consists in solving a maximal piece in direction d problem (M P d ):

•

  the direction going along the bisector of the two edges of D(C) having v as starting point, denoted bd for Bisector Direction • Compute two maximal linear segments starting from v and following each edge of D(C) having v as endpoint. The direction orthogonal to the line joining the two ends of the maximal linear segments is chosen as the progress direction d. It is denoted med for Mean progress along Edges Direction

Idea 3: Computing a Feasible Solution of a Max- imal Piece in Direction d Problem

  on which ũ and l are linear, indexed by I and J respectively. Then, (M P d ) is a relaxation of the following problem because C is replaced by an inner corridor.

	As we want a computationally cheap solution to the GSIP (M P d ) because of
	the high number of times such a problem has to be solved, a feasible solution
	of (M P d ) is computed via a series of LP problems, as described below, using a
	PWL inner corridor of the original corridor C, because it allows to replace the
	infinite number of nonlinear constraints by a finite number of linear constraints
	while ensuring feasibility.
	Let corridor C PWL ũ and l for readability; note (D i d σ be a PWL inner corridor of C d σ with associated functions ũ) i∈I and (D j l ) i∈J the subdomains of corridor
	C PWL	d σ

Table 1

 1

		: Expression and domain of benchmark functions
	ref expression	Domain	ref	expression	Domain
	L1

Table 2 :

 2 Comparison with state of the art heuristics

			DN95		DN99	RK2D	LinA2D		KL2D
	ref	δ	n	time	n	time	n	time	n	time	n	time
		1.5	8	20.6	7	44.7	16	30.8	6	0.0	6	87.1
		1.0	10	28.2	9	67.6	20	84.4	8	0.0	6	148.9
	L1	0.5	22	54.5	20	175.1	48 150.4	15	0.0	-	TO
		0.25	46	110.9	42	258.2	80 272.6	21	0.0	-	TO
		0.1	113 356.8 106 762.9 224 380.6	60	0.0	-	TO
		1.5	8	29.6	6	39.1	24	26.8	6	0.0	-	TO
		1.0	13	39.7	11	101.2	28	7.4	8	0.0	6	1495.5
	L2	0.5	26	47.6	27	117.9	84	38.0	15	0.0	-	TO
		0.25	59	114.0	53	279.7 121 35.8	21	0.0	-	TO
		0.1	146 287.2 136 642.5 351 171.7	60	0.0	-	TO
		1.0	4	9.5	3	21.2	4	0.8	6	0.0	4	18.5
		0.5	10	22.4	6	52.4	12	72.4	8	0.0	6	416.8
	N1 0.25	16	70.2	14	116.1	20	4.7	18	0.0	12 14762.5
		0.1	32	115.8	28 324.5	59	59.3	45	0.0	-	TO
		0.05	68	293.2	56 564.6	94	45.3	91	0.0	-	TO
		0.1	1	4.3	1	4.1	2	0.3	4	0.0	1	6.0
		0.05	2	14.6	2	47.8	6	18.7	9	0.0	2	13.4
	N2 0.03	4	24.6	4	57.0	10	12.7	12	0.0	4	39.0
		0.01	11	87.9	11 450.4	31	54.6	30	0.0	-	TO
		0.001 105 1041.9	96 3832.1 350 652.6 238	0.0	-	TO
		1.0	3	10.4	3	17.1	5	1.0	12	0.0	-	TO
		0.5	4	10.0	4	37.4	8	13.1	16	0.0	-	TO
	N3 0.25	7	19.5	6	42.2	16	30.0	22	0.0	8	342.6
		0.1	18	57.5	16 306.1	44	74.6	68	0.0	-	TO
		0.05	34	193.0	30 681.2	85 141.9 120	0.0	-	TO
		0.5	2	15.5	2	23.2	2	1.8	3	0.0	2	13.3
		0.25	3	9.2	4	57.8	4	1.0	8	0.0	4	19.6
	N4	0.1	7	41.2	6	230.0	9	25.8	21	0.0	6	66.9
		0.05	14	161.7	11 632.6	23	14.4	27	0.0	12	7246.5
		0.03	21	272.7	19 1066.4 40 161.4 44	0.0	-	TO
		1.0	2	8.3	1	5.0	6	7.7	20	0.0	1	6.3
		0.5	8	33.5	5	96.5	6	1.3	42	0.0	4	86.9
	N5 0.25	12	49.9	14	242.7	21	30.8	72	0.0	5	2091.2
		0.1	26 103.5	26 646.9	96	73.0 168	0.0	-	TO
		0.05	56	313.4	51 954.9 274 305.5 340	0.0	-	TO
		1.0	3	11.7	3	51.4	6	22.8	9	0.3	3	23.3
		0.5	5	22.7	5	78.6	9	15.6	9	0.0	4	127.2
	N6 0.25	9	42.7	8	188.2	12	22.9	16	0.0	6	305.5
		0.1	16 136.7	16 803.9	40 202.8 49	0.0	-	TO
		0.05	34	516.9	28 1986.8 87 174.1 81	0.0	-	TO
		1.0	1	8.1	1	23.8	2	0.8	4	0.0	1	2.9
		0.5	2	32.8	2	234.8	4	66.6	4	0.0	2	9.5
	N7 0.25	4	145.7	4 1370.6	6	4.4	9	0.0	4	35.5
		0.1	6	462.4	7	6210.7 84 231.5 16	0.0	5	377.7
		0.05	10 2050.0	-	TO	86	57.8	36	0.0	-	TO


Victor Fiorese
email: victor.fiorese1@st.com

Joao Carlos

Azevedo Goncalves

Simon Bouvot

Emmanuel Dubois

Christophe Gaquière

Guillaume Ducournau

François Danneville

Sylvie Lépilliet

Daniel Gloria

Christophe Gaquiere

A 140 GHz to 170 GHz active tunable noise source development in SiGe BiCMOS 55 nm technology

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

In the 5G development context, the increased need for accurate measurements at high frequencies (HF) has brought new concerns regarding testing setups. Indeed, the operating frequencies of recently developed transistors have exceeded 370 GHz for f [START_REF] Chevalier | SiGe BiCMOS Current Status and Future Trends in Europe[END_REF], introducing the development of new mmW circuits. The design of such systems relies on accurate and reliable models of MOS or bipolar transistors that are obtained by small signal, large signal, and HF noise characterization. The usual frequency range of such noise models is limited to 110 GHz. Beyond this frequency limit, measurements are challenged by important interconnect losses and the difficulty to respect the minimum detectable signal (MDS) condition which is related to the noise receiver noise figure (NF ). Moreover, by operating above 110 GHz comes the difficulty to generate a sufficient Excess Noise Ratio (ENR) that is directly used during the Y-factor measurement method. This method is commonly used for noise measurement and consists in determining the ratio of measured noise power when the NS is in ON state P and when it is in OFF state P . The NS is based on a B55 noise diode [START_REF] Azevedo Goncalves | Millimeter-Wave Noise Source Development on SiGe BiCMOS 55-nm Technology for Applications up to 260 GHz[END_REF] biased in avalanche regime for the ON state. Such noise sources can generate around 17 dB of ENR at 130 GHz [START_REF] Goncalves | Développement de bancs de caractérisation pour la mesure de bruit et la détection de puissance entre 130 GHz et 260 GHz[END_REF]. Output noise power can prove to be insufficient to comply with noise receiver's sensitivity [START_REF] Bouvot | Contribution au BIST in-situ : Intégration sur silicium d'un banc de caractérisation en bruit en bande D[END_REF]. Indeed, the usual NF value for offwafer noise receiver in these frequency ranges is 10 dB [START_REF] Deng | Contribution à la caractérisation et la modélisation jusqu'à 325 GHz de transistors HBT des technologies BiCMOS[END_REF]. As suggested in [START_REF] Goncalves | Développement de bancs de caractérisation pour la mesure de bruit et la détection de puissance entre 130 GHz et 260 GHz[END_REF], using a biased state instead of an unbiased state for noise source low state can help to validate the MDS condition. The drawback of such an approach is a decrease in the ENR maximum reachable levels, given that the delta between the two equivalent noise temperatures T and T corresponding respectively to the ON and OFF NS states is smaller. Another solution can be the approach consisting in integrating a low noise amplifier (LNA) to have a constant output NS matching regardless of bias condition and higher ENR levels. Several works leveraging photodiodes [START_REF] Vidal | Broadband Photonic Microwave Noise Sources[END_REF], [START_REF] Oeuvrard | On Wafer Silicon Integrated Noise Source Characterization up to 110 GHz based on Germanium-on-Silicon Photodiode[END_REF], NMOS transistors [START_REF] Alimenti | Avalanche Microwave Noise Sources in Commercial 90-nm CMOS Technology[END_REF] or p-i-n diodes [START_REF] Alimenti | Millimeter-Wave Avalanche Noise Sources Based on p-i-n Diodes in 130 nm SiGe BiCMOS Technology: Device Characterization and CAD Modeling[END_REF] showed conclusive results up to 40 GHz. In this work, the addressed frequency range is above 140 GHz. Limitations of HF noise measurement are detailed in section II while section III describes the solution proposed to ensure the MDS validation and a sufficient ENR level at the same time. Section IV describes the methodology used for ENR extractions. Finally, section V provides some results and prospects.

II. HIGH FREQUENCY NOISE MEASUREMENTS CHALLENGES OVERVIEW

A. DUT Noise factor measurement with Y-factor method

The HF noise calibration is performed by connecting the noise source to the input of the noise receiver and making Yfactor measurement [START_REF] Bouvot | Contribution au BIST in-situ : Intégration sur silicium d'un banc de caractérisation en bruit en bande D[END_REF]. During the calibration, the noise figure (F ) and gain (G ) of the receiver are extracted allowing the calibration in plane B. Once the noise receiver is calibrated, the DUT can be characterized when inserted between noise source and noise receiver as shown in Figure 1. The noise figure of the DUT ( F " can be determined using the Y-factor method [START_REF] Technologies | Fundamentals of RF and Microwave Noise Figure Measurements[END_REF]. To perform noise characterization above 110 GHz with sufficient accuracy, the MDS condition must be validated. This step follows the noise receiver calibration and is necessary to perform a noise characterization with sufficient responsivity and accuracy. This condition described in (1) is based on the equivalent noise temperature T " and the available power gain G " of the DUT. With k ; being the Boltzmann constant, T < ' 290K, T @A the equivalent temperature of the noise source, ∆ B the frequency bandwidth, T the ambient temperature, close to T < . The MDS condition ensures that the output noise power signal of the noisy system (P " must exceed a certain value depending on the noise receiver NF in order to be detectable with a margin of 3 dB. Usually, at frequencies higher than 110 GHz, 10 dB is a relevant value for NF . In the case of a B55 high speed HBT noise characterization, the DUT contribution to G " is not significant since this kind of component has roughly 1 dB of gain at 200 GHz. In the case of the Y-factor method, the two NS operating modes (ON and OFF) must validate the MDS criteria, the OFF state being the most challenging one.

B. The ENR level requirement

A noise source is always associated to an ENR level, defined as the noise power that a noise source can generate compared to the noise floor power of -174 dBm/Hz. This notion is introduced by the commonly used method of noise factor measurement involving the two operating states ON and OFF of a noise source, assigned to two noise temperatures T and T

. From these two temperatures, a definition of the ENR can be introduced by ( 2), with T0=290 K.

CDE FG =10IJK & LMN O& PMQR & S .
(

) 2
Commercial packaged noise sources have shown ENR levels of 12 dB up to 220 GHz [11] with a 2 dB flatness. Other works have developed on-wafer noise sources based on SiGe BiCMOS 55 nm noise diodes, showing monitorable ENR levels up to 260 GHz [START_REF] Ghanem | Modeling and Analysis of a Broadband Schottky Diode Noise Source Up To 325 GHz Based on 55-nm SiGe BiCMOS Technology[END_REF]. Such noise sources generate noise when biased in avalanche regime. By changing the biasing current value, several ENR levels can be reached, from 6 to 20 dB in the 130-170 GHz frequency range. Following the trend of these works, a new idea of an active noise source composed of a diode and a LNA is proposed. Such a solution can show a variable ENR level while respecting the MDS condition.

III. TEST STRUCTURES AND CHARACTERIZATION SETUP DESCRIPTION

A. Noise source and LNA description

The used NS is a SiGe BiCMOS 55 nm diode showing monitorable ENR levels up to 17 dB at 260 GHz when biased in avalanche regime [START_REF] Azevedo Goncalves | Millimeter-Wave Noise Source Development on SiGe BiCMOS 55-nm Technology for Applications up to 260 GHz[END_REF]. The used LNA design and characterization well described in [START_REF] Bouvot | Contribution au BIST in-situ : Intégration sur silicium d'un banc de caractérisation en bruit en bande D[END_REF], has a 15 dB flat gain between 130 GHz and 160 GHz. It is structured as a classic cascode structure on two stages. Each stage is composed of two HBT NPN with an emitter width and length of respectively 0.2 µm and 5.56 µm. Since ENR measurements were performed between 130 GHz and 170 GHz, this LNA appears to be a relevant candidate for this study and is used as a buffer stage to enhance ENR levels. A layout view of this LNA is proposed in Figure 2. This LNA has been characterized in Sparams and noise in D band (140-170 GHz), showing a NF lower than 8 dB between 135 GHz and 170 GHz as described in Figure 3. The lower NF level corresponds to the NF of the LNA de-embedded of the RF pads.

B. ENR extraction when using a noise source and LNA setup

In order to assess the interest of adding a LNA to increase the ENR value, complementary extractions were performed in the setup illustrated on Figure 4: This setup has been used to extract ENR levels in several planes A and B to analyse the ENR level contribution of each element. To be able to extract ENR levels in plane A, measurements of the NS and LNA were performed in ON and OFF states. Besides this noise power level aspect, another advantage of using the chosen LNA is to greatly improve the output matching, regardless of NS bias current, due to its IV. METHODOLOGY In order to calculate the ENR, calibrated noise powers were measured in ON and OFF states at the input noise receiver plane [C]. To take into account the mismatch between the components, reflection coefficients and available gains calculation of the several blocks was made. Figure 2 shows the several parameters calculated to extract the ENR in planes A and B. At first, calculation of each reflection coefficient and available gain is done. This is possible due to the knowledge of the noise receiver and noise source reflection coefficients Γ and Γ Aab Ac for each bias condition. Moreover, probe S-parameters are given by probe supplier and output pad is not de-embedded in this study. The first step consists in calculating Γ de fg and G de fg with the knowledge of Γ Aab Ac and Γ Aa , determined during Sparameters measurements. During this development, the noise receiver is assimilated to be a 50 Ω load termination. This permits to calculate the mismatch factor (M) [START_REF] Miller | Noise Standards, Measurements, and Receiver Noise Definitions[END_REF], [START_REF] Dunleavy | Characterization and Applications of On-Wafer Diode Noise Sources[END_REF] used to transform effective noise powers into available noise powers, to be able to extract ENR levels. With this knowledge, the transform of relative measured effective noise powers (P gBB ) compared to -174 dBm is done to obtain equivalent available noise temperatures at noise receiver input plane (T ,; ). (4)

Then, the equivalent available noise temperatures at probe input plane (T ,; is calculated thanks to the knowledge of G de fg as shown in equation ( 5 In this study, several bias states were used for NS and LNA. Extractions of ENR levels at plane A (see Figures 1 and4) expressed in equation ( 2) have been done. T ,c, is a bias state where NS is OFF and LNA is ON, T ,c, is a bias state where NS is biased with currents varying from 0.25 mA to 25 mA and LNA is ON. It is to be highlighted that unlike conventional ENR definition, for which T is usually equal to the room temperature, used T values in our study are varying between 39000 K and 73000 K (see Figure 6). V. RESULTS Extracted ENR levels are represented in Figure 7, showing highly monitorable powers up to 37 dB. To our knowledge, these are the first results of such monitorable ENR levels at these frequencies using SiGe technology. Considering the MDS aspect, the use of NS in low state and LNA in high state shows a validation of the MDS condition with a 10 dB margin in the 140-170 GHz frequency band. When NS and LNA are set in low state, MDS condition can obviously not be satisfied as depicted in Figure 8. A mmW and sub-mmW noise sources state-of-the-art is presented Table 1.

Table 1. mmW and sub-mmW noise sources state-of-the-art VI. CONCLUSION A new approach of in-situ active noise source integrated on silicon using a D-band LNA has been described. This way of operating induces higher monitorable ENR levels while respecting the MDS condition. A noise source low biasing LNA high biasing proves to be sufficient to validate this condition with a 10 dB margin. This active noise source shows a constant output matching in both high and low states of biasing. In the case of multi-impedance noise characterization of a DUT, an impedance tuner may be used in source-pull configuration. As a perspective, further studies will bring the noise measurement setup approach of this active noise source followed by an impedance tuner and a DUT in order to evaluate the multiimpedance method.

Technology

Fig. 1 .

 1 Fig. 1. Classical setup for NF measurement leveraging Y-factor method

Fig. 2 .

 2 Fig. 2. Layout view of the LNA used in the present active noise source [4].

Fig. 3 .

 3 Fig. 3. NF measured in D-band of the 2 stage LNA [4].

Fig. 4 .

 4 Fig. 4. Block diagram of the noise source and LNA ENR extraction setup

5 .Fig. 5 .

 55 Fig. 5. Output matching comparison of NS (solid lines) and NS+LNA (dot lines) between 140 GHz and 170 GHz for 3 NS biasing.

Fig. 6 .

 6 Fig. 6. Used values for T ,c, K in probe input plane (A plane).

Fig. 7 .

 7 Fig. 7. Extracted ENR dB in probe input plane (A plane).

  with

Fig. 8 .

 8 Fig. 8. MDS condition compared to Pnoise eff (dBm) in B plane