Models and Algorithms for Optimal Piecewise-Linear Function Approximation

Piecewise-linear functions can approximate nonlinear and unknown functions for which only sample points are available. This paper presents a range of piecewise-linear models and algorithms to aid engineers to find an approximation that fits best their applications. The models include piecewise-linear functions with a fixed and maximum number of linear segments, lower and upper envelopes, strategies to ensure continuity, and a generalization of these models for stochastic functions whose data points are random variables. Derived from recursive formulations, the algorithms are applied to the approximation of the production function of gas-lifted oil wells.


Introduction
Piecewise-linear functions are widely used to approximate functions for which only sample points are known and to model nonlinear functions.In petroleum engineering, the fluid flow from an oil well and the pressure drop in a pipeline can be approximated with a piecewise-linear function [1].In optimization, nonlinear problems can be recast as a mixedinteger linear programming (MILP) problem, which is then solved with MILP algorithms [2][3][4].
Although many works in the literature are dedicated to piecewise linearization, few of them address the problem of generating piecewise-linear models which is often delegated to the scientist or engineer.Some techniques are based on point clustering, such as -means, which has found applications in nonlinear dynamic systems and control [5].This work attempts to mitigate undesirable effects of clustering, particularly the influence of outliers and the synthesis of subpotimal approximations.Akin to the previous technique, a fuzzy clustering strategy is developed in [6] to represent a nonlinear system with a piecewise-linear model.However, fuzzy clustering is dependent on the initial parameters and can reach a local minimum.In order to reduce the chance of yielding a local minimum, genetic algorithms have been proposed for clustering [7].Other works rely on local searches that introduce linear segments iteratively until satisfying a stopping criterion [8,9].
However, all of the cited approaches do not offer a certificate of optimality and do not allow the user to specify desirable properties for the resulting approximation.To this end, this paper develops formal models and algorithms for generating piecewise-linear functions that fit best the needs of the end user.The models consider a fixed and maximum number of linear segments, lower and upper envelope approximations, constraints to ensure continuity, and stochastic generalizations of these models in which the data points are random variables.The underlying approximation problems for these models are cast in a recursive form that can often be efficiently solved with dynamic programming strategies.To the best of the knowledge of the authors, this is the first work on the synthesis of piecewise-linear functions based on formal models and algorithms which considers an extensive range of conditions and approximation metrics.Software tools are obtained from these strategies to assist

Basic Models and Algorithms
Let  = { 1 , . . .,   } be a set of points of a function Φ : R → R such that   = (  ,   ),   = Φ(  ), and  1 <  2 < ⋅ ⋅ ⋅ <   .Φ may represent any function, hereafter referred to as generating function, which can be known explicitly, implemented by a simulator, or obtained from a real-world process such as the output flow from an oil well as a function of the compressed gas injection.The problem is to find a set  = { 1 , . . .,   } of linear segments that best approximates Φ.Assuming that the segment breakpoints are points in , a piecewise-linear model is a sequence  = ⟨(1), (2), . . ., ( + 1)⟩ of breakpoints such that An approximation for () is obtained by minimizing the squared error between the linear segment   and the points in (), for which an analytical solution yields the intercept   and the slope   .Table 1 presents a set of sample points (the sample points were obtained in an arbitrary manner by perturbing the oil production function  oil ( inj ) of an oil well induced by the lift-gas injection rate  inj ; the purpose of these sample points is only to illustrate the behavior of the models and algorithms to be developed in this work) that will serve to illustrate models and algorithms.

Fixed Number of Linear Segments.
To find the piecewiselinear approximation with a fixed number  of linear segments that minimizes the squared error, the following problem is solved: This problem can be recast in a recursive form and solved with a dynamic programming (DP) algorithm.For an introduction to dynamic programming in discrete domains, the interested reader can refer to the text books by Kleinberg and Tardos [10] and Cormen et al. [11].The book by Bertsekas [12] is recommended for an advanced treatment of DP with particular emphasis on continuous domains.Let (, ) be the minimum cost to approximate the points { 1 ,  2 , . . .,   } with exactly  segments.Then, The optimal sequence  = ⟨(1), (2), . . ., ( + 1)⟩ is obtained recursively: ( + 1) = ; () is the index  that minimizes ((+1), ) in (2); (−1) is the index  that minimizes ((), −1); and so forth until ((2), 1) is reached, at which point (1) = 1.Algorithm 1 outlines the DP algorithm PL-Fixed for computing (, ).The algorithm records in table (, ) the indices  that minimize (2), which are the optimal breakpoints.PL-Fixed is correct because it solves recursive equation (2).It runs in Θ( 3 ) time since  < .The notation Θ(ℎ()) = {() : ∃ 0 ,  1 ,  2 > 0, such that  1 ℎ() ≤ () ≤  2 ℎ() for all  ≥  0 } represents the set of all nondecreasing functions () that grow at the same asymptotic rate of ℎ().Therefore, the time taken by algorithm PL-Fixed to solve an instance of size  is approximately  3 .A precise and extensive discussion on the running-time complexity of algorithms is found in [11].
Recursive equation ( 3) leads to a DP algorithm PL-Max, which is omitted here in the interest of brevity.PL-Max(, , ) runs in Θ( 3 ) provided that  ≤ .This algorithm records the index  corresponding to the optimal breakpoint of (3) in a table (, ).This table is key to determine the optimal number  ⋆ of segments and construct the optimal sequence ⟨(1), (2), . . ., ( ⋆ + 1)⟩.

Continuous Approximations
The preceding models may produce discontinuous piecewiselinear functions.To enforce continuity either because the actual function is known to be continuous or because continuity is desirable for numerical optimization, some adaptations are proposed for the basic models.

Ensuring Continuity at the Breakpoints.
Continuity is enforced by approximating the set  , = {  ,  +1 , . . .,   } with segments going through   and   .The recursive formulations ensuring continuity with a fixed and maximum number of segments are identical to the preceding ones, except that  , is approximated with the segment passing through   and   .The theorem below is a consequence of this observation.
Theorem 3. Suppose that the coefficients of the linear segment approximating  , are  , = (  −   )/(  −   ) and  , =   −  ,   .Then, recursive formulations ( 2) and ( 3) produce optimal piecewise-linear approximation functions for  that are continuous at the breakpoints using exactly and at most  linear segments, respectively.
Applying PL-Max to the sample instance with a maximum of  = 5 segments and cost  = 0.01, but forcing the segments approximating  , to go through the end points, yields the piecewise-linear function with  ⋆ = 4 segments, which is illustrated in Figure 2.This approximation is given by the sequence ⟨(1), . . ., (5)⟩ = ⟨1, 2, 4, 6, 8⟩ which, coincidentally, is identical to the one obtained without ensuring continuity.However, its approximation cost (8, 5) = 0.046556 exceeds the optimal cost 0.044208 of the discontinuous approximation, which is expected because a discontinuous model encompasses its continuous counterpart.This property is formally stated by the proposition below.
Proposition 4. The approximation cost induced by the schemes that enforce continuity at the breakpoints, with a fixed and maximum number of segments, is not lower than the cost induced by the schemes that do not enforce continuity.

Rejecting Discontinuous Approximations.
The basic models with fixed and maximum number of segments can be modified to ensure continuity by rejecting discontinuous approximations.The recursive equation for generating a piecewise-linear function with a fixed number  > 1 of segments is Algorithm 2 gives the dynamic programming solution of recursive equation ( 4), for a fixed number  of segments.The algorithm records in table (, ) the points where the linear segments intercept, which are used to determine the breakpoints of the piecewise-linear function.The algorithm runs in O( 3 ) time provided that  = O().

Generating Continuous Approximations.
Continuity may be ensured by solving a least squares problem at the moment that the recursive equation is evaluated.For a fixed number of segments  > 1, recursive equation (2) becomes while for a maximum number of segments (3) becomes where (, ) ⋆ is obtained by solving the following problem: where  ,−1 and  ,−1 are the slope and intercept of the rightmost segment of the optimal approximation of  1, with exactly or at most  − 1 segments, depending on which recursion is being solved.(, ) ⋆ = ∞ if the problem is infeasible, for example, when  = .The computation of (, ) ⋆ entails minimizing the squared error of the approximation of {  ,  +1 , . . .,   } with a linear segment, while ensuring that this segment intercepts the preceding one, the rightmost segment of the piecewise-linear approximation of { 1 , . . .,   }, in the interval [ max{1,−1} ,  min{+1,} ].The chief difficulty lies in the nonlinearity of (7b).However, the optimal solution is obtained by solving two quadratic programs.First, notice that variable  is eliminated by replacing the constraints (7b)-(7c) with the following nonlinear inequality:

(9b)
The end result is that (, ) ⋆ can be computed by solving two convex, quadratic programs (QPs) and choosing the least costly solution.The QPs minimize the objective (7a), with the first being subject to constraint (9a), whereas the second is constrained by (9b).Dynamic programming algorithms are obtained to solve recursive equations ( 5) and ( 6) using a standard QP solver.The algorithms run in O( 3 ) time, with  being an upper bound on the cost for solving QPs with 2 variables and 3 constraints.The application of the DP algorithm for the sample instance with a fixed number of  = 3 segments yields the sequence ⟨(1), . . ., (4)⟩ = ⟨1, 5, 6, 8⟩.The breakpoints of the segments are defined by their intercepts, which are [ 1 , 0.3120] = [0.0114,0.3120] for the first segment, [0.3120, 0.3215] for the second, and finally [0.3215,  8 ] = [0.3215,0.9104] for the third.The approximation cost is (8, 3) = 2.3147 × 10 −2 .

Lower and Upper Envelopes
The preceding models can be modified to obtain piecewiselinear lower and upper envelopes of the generating function Φ.For the upper envelope, the slope  , and intercept  , should minimize the squared error of the line  =  ,  +  , approximating  , = {  ,  +1 , . . .,   }, while keeping these points below the line.Formally, this is achieved by solving the QP:  , : min The basic algorithms can compute the optimal upper envelope, with a fixed and maximum number of segments, by solving  , to obtain the coefficients of the segment approximating  , .The theorem below formalizes these results for the upper envelope approximation.Algorithms are developed similarly for the lower envelope by changing the direction of inequality (10b).

Theorem 5. Recursive equations (2) and (3) yield the optimal piecewise-linear upper envelope approximation of 𝑃 with exactly and at most 𝑡 linear segments, respectively, having the squared approximation error as the objective, provided that the coefficients for the linear segments are computed by solving problem 𝐿 𝑖,𝑗 given by (10a) and (10b).
The running time of the resulting algorithm is O( 3 ()) where () is an upper bound on the cost to solve the quadratic program  , with 2 variables and  constraints.
For the algorithm that computes linear segments by solving QP problem (7a), (7b), and (7c), we only need to introduce the linear inequalities that ensure bounding, such as the constraints (10b) for an upper envelope approximation.

Dealing with Random Data
Uncertainty is addressed by modeling each point   ∈  as a random variable p  = (x  , y  ) characterized by a probability density function   (, ).Two models are proposed, one minimizing expected squared error and the other using chance constraints [13].

Minimizing Expected Squared
Error.Uncertainty may be treated by minimizing the expected squared error of the piecewise-linear approximation.For the points in the set P , = {p  , p +1 , . . ., p  }, this approach finds a linear segment given by the slope  , and intercept  , that solves the problem: ) is convex on the variables ( , ,  , ) since the term ( ,  +  , − ) 2 is convex for fixed (, ), generalizing an infinite weighted sum [14,Sec. 3.2.1].This means that  , is a convex function and so is problem  , .Hence, an optimal solution to  , can be found with numerical optimization methods, such as the interior-point method, using numerical differentiation procedures or directly when a closed form is known for the gradient.
The basic algorithms of Section 2 yield the optimal approximation, for a fixed and maximum number of segments, by solving  , and using the obtained coefficients for the linear segment approximating P , .This result is formalized in the theorem below, which can be easily proven.Theorem 6. Recursions ( 2) and ( 3) yield the optimal piecewise-linear approximation of a set P = {p 1 , . . ., p  } of random points with exactly and at most  segments, respectively, having the expected squared approximation error as the objective if the coefficients for the linear segments are obtained by solving  , given by ( 11a) and (11b).
The resulting algorithms run in O( 3 ()) time, with () being an upper bound on the running time to solve  , with 2 variables and  integration terms.

Lower and Upper
Envelopes.Instead of minimizing only the expected squared error, approximations for the lower and upper envelope can be obtained by introducing probabilistic constraints to problem  , .An algorithm for the lower envelope is developed below, which is easily adapted for the upper envelope.
To find the lower envelope of P , = {p  , p +1 , . . ., p  }, the constraint that the random point p  = (x  , y  ) is above the where P(⋅) denotes the probability operator.The DP algorithms can find the lower envelope by changing how the slope  , and intercept  , of the linear segments are computed.The problem of finding the segment coefficients that minimize the expected squared error, subject to the probabilistic constraints for the lower envelope, is while being subject to P ( ,  +  , ≤ y  ) ≥ ,  ∈   ,  = , . . ., .
By expressing probabilistic constraint (14c) in terms of integrals and the joint probability density function   , a deterministic equivalent to this problem is obtained: where   , is the probability of p  lying below the linear segment that approximates P , .This problem may be posed in terms of the probability   , of p  being above the linear segment.
Problem  le , is a nonlinear program for which algorithms reach a locally optimal solution, such as sequential quadratic programming (SQP).
The basic algorithms can find the optimal approximation, for a fixed and maximum number of segments, by solving  le , and using the obtained coefficients in the approximation of P , .This result is formalized in the theorem below, whose proof is similar to that of Theorem 1.
Theorem 7. Recursions ( 2) and ( 3) yield the optimal piecewise-linear approximation of a set P = {p 1 , . . ., p  } of random points with exactly and at most  segments, respectively, having the expected squared approximation error as the objective and subject to chance constraints to probabilistically ensure the lower envelope, if the coefficients for the linear segments correspond to the optimal solution of problem  le , .

Computational Analysis
In order to keep the analysis brief, the experiments were restricted to the following approaches: (1) PL-Max, the basic model of Section 2.2 which enforces a maximum number of linear segments.
(2) PL-Max-Cont, the version of the basic model that imposes continuity by solving quadratic programs according with the developments of Section 3.3.
(3) PL-Max-Bound, the version of PL-Max which minimizes the number of linear segments while ensuring a bound on approximation error.
The algorithms were implemented in C++ in a workstation equipped with an Intel i7 CPU, 2.67 GHz, with 8 GB of RAM, and running the Linux operating system version Ubuntu 2.6.35-32.
The CFSQP software package and its embedded QP solver were used to solve optimization problems [15].CFSQP follows the sequential quadratic programming (SQP) framework which solves a sequence of QPs approximating the original problem but having the Lagrangian as the objective to capture the curvature.The gradient of the objective and constraints were provided analytically.The Hessian of the objective was computed numerically using the BFGS strategy.
The approaches were put to the test in the piecewise-linear approximation of the production function of gas-lifted oil wells.As an oil field matures, the reservoir internal pressure reduces to a level that is not sufficiently high to sustain natural production, requiring the application of artificial lifting techniques.Among these techniques, gas-lift is widely used by the industry for its robustness and relatively low operating costs [16].Gas-lift works by injecting high pressure gas at the bottom of the production tubing, reducing the fluid density due to gasification and thereby lifting the fluid to the surface.
For the computational analysis, the generation function Φ corresponds to the total mass flow   (  inj ) of oil, water, and gas given as a function of the lift-gas rate   inj injected in an oil well .The relation   (  inj ) is referred to as well production function.This work adopted the oil production system presented in [4], which consists of subsea gas-lifted wells, pipeline-risers, and topside processing equipment such as compressors and separators.The production system is modeled in Pipesim, a commercial multiphase flow simulator widely adopted by the oil industry.
A detailed piecewise-linear model was obtained by uniformly sampling the simulation model of   for a sufficiently high number of   inj values.In our experiments, 50 breakpoints produced a highly accurate and detailed description of   .However, simpler and sufficiently accurate piecewiselinear approximations can be obtained by a better distribution of the breakpoints within the function domain, a task to be performed by the algorithms proposed in this work.
Table 4 shows the running time of the three implemented strategies for piecewise-linear approximation of 10 oil-well production functions, each originally represented with 50 breakpoints.For PL-Max and PL-Max-Cont, the maximum number of segments was set to  = 10 and the cost per segment to  = 0.For PL-Max-Bound,  = 10 and the approximation bound was  = 10.The CPU times are consistent with the time complexity of the algorithms: PL-Max is faster than PL-Max-Bound due to the extra constraints of the latter; in turn, PL-Max-Bound is faster than PL-Max-Cont because he latter needs to solve QP problems.
The behavior of the PL-Max strategy is shown in Table 5 for a varying cost  and a maximum number of segments  = 15.The algorithm was applied to the production function of well #1.The squared error grows as  increases, corroborating the theoretical behavior.
The behavior of PL-Max-Cont strategy is illustrated in Table 6 for a varying number of maximum segments.The algorithm was also applied to the production function of well #1.As expected, the squared error decreases as the limit  increases.
The computational experiments show that piecewiselinear functions respecting the desired criteria can be obtained with fewer breakpoints than the original representation.Such simpler piecewise-linear models can have a major impact on the production optimization of large oil fields [17].

Summary and Conclusion
Applications of piecewise-linear models abound in science and engineering.In petroleum engineering, the performance curve of an oil well can be approximated with a piecewiselinear function.Owing to their importance, this paper developed models for piecewise linearization of nonlinear functions that trade off model complexity and approximation quality.The models and algorithms serve as decision-support tools for engineers to find piecewise-linear models that fit best their needs.
Several models were developed for the synthesis of piecewise-linear approximations under a variety of conditions.Given a potentially large set  = { 1 , . . .,   } of graph points generated by a function Φ, not necessarily known explicitly, the problem consists in finding a reduced subset of points that induces a more compact, yet sufficiently accurate piecewise-linear approximation of Φ.To this end, a number of criteria were proposed in order to formally state the approximation problems for which algorithms have been designed.The main contributions of this paper consist of the following formulations and algorithms for piecewiselinear approximation: (1) A model and dynamic programming algorithm for finding the optimal approximation using a fixed number of linear segments, which minimizes a quadratic error with respect to the point set .
(2) A model and DP algorithm to find the optimal approximation using a maximum number of segments, which minimizes the quadratic approximation error and a fixed cost per segment.
(3) Variations of the two preceding models and algorithms, referred to as basic approaches, to ensure a continuous piecewise-linear approximation.
(4) Variations of the basic approaches that yield optimal upper and lower envelopes for Φ.
(5) Extensions of the preceding approaches to handle uncertain data, whereby the points of Φ are random variables.
The developed approaches can be combined or extended to create models and algorithms that are tailored to the user's needs, such as the synthesis of a piecewise-linear function approximating the upper envelope, with continuity and in the presence of random data points.The proposed models measured error with the  2 -norm; however other norms can be easily adopted such as the  1 -and  ∞ -norm by solving suitable linear and quadratic programs.

Figure 2 :
Figure 2: Continuous piecewise-linear approximation of the sample instance with a maximum of 5 segments.

Figure 3 :
Figure 3: Continuous, upper envelope approximation with 3 linear segments obtained with the rejection strategy.

Table 2 :
Tables with values (, ) and (, ) for the sample instance.

Table 4 :
Running time of the algorithms (in seconds).

Table 5 :
Behavior of the PL-Max strategy.

Table 6 :
Behavior of the PL-Max-Cont strategy.