# Models and Algorithms for Optimal Piecewise-Linear Function Approximation

Authors: Eduardo Camponogara; Luiz Fernando Nazari
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

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.


---

## Body

## 1. 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 mixed-integer linear programming (MILP) problem, which is then solved with MILP algorithms [2–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 asK-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 scientists and engineers in the synthesis of piecewise-linear approximations of nonlinear models and processes.Sections2 through 5 develop a range of models and algorithms for piecewise-linear approximation. Section 6 reports computational results on the approximation of well production in gas-lifted oil fields. Section 7 offers a summary and suggests directions for future research.


## 2. Basic Models and Algorithms

LetP={p1,…,pn} be a set of points of a function Φ:R→R such that pi=(xi,yi), yi=Φ(xi), and x1<x2<⋯<xn. Φ may represent any function, hereafter referred to asgenerating 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 R={r1,…,rT} of linear segments that best approximates Φ. Assuming that the segment breakpoints are points in P, a piecewise-linear model is a sequence I=〈i(1),i(2),…,i(T+1)〉 of breakpoints such that(i)
i(1)=1<i(2)<⋯<i(T+1)=n;(ii)
each subsetP(t)={i(t),i(t)+1,…,i(t+1)}⊆P is approximated by the linear segment rt={(x,y):y=atx+bt,xi(t)≤x≤xi(t+1)} for t=1,…,T.An approximation forP(t) is obtained by minimizing the squared error between the linear segment rt and the points in P(t), for which an analytical solution yields the intercept bt and the slope at. Table 1 presents a set of sample points (the sample points were obtained in an arbitrary manner by perturbing the oil production function qoil(qinj) of an oil well induced by the lift-gas injection rate qinj; 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.Table 1
Set of sample points.
i
x
i
y
i
1
0.0114
0.1857
2
0.0861
0.4502
3
0.1721
0.5584
4
0.2928
0.5857
5
0.3130
0.7020
6
0.3212
0.8338
7
0.4899
0.8713
8
0.8807
0.9104


### 2.1. Fixed Number of Linear Segments

To find the piecewise-linear approximation with a fixed numberT of linear segments that minimizes the squared error, the following problem is solved:(1)minIE=∑t=1TEt=∑t=1T∑k=itit+1yk-atxk-bt2.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 F(j,t) be the minimum cost to approximate the points {p1,p2,…,pj} with exactly t segments. Then,(2)Fj,t=mint≤i≤jFi,t-1+Ei,j,where E(i,j) is the minimum squared error of approximating Pi,j={pi,pi+1,…,pj} with just one segment. Notice that E(i,i)=∞, F(j,t)=∞ for all j≤t, and F(j,1)=E(1,j) for all j>1. From recursive equation (2), clearly F(n,T)=minE.Theorem 1.
Recursive equation (2) yields the piecewise-linear approximation of P of minimum squared error having exactly t linear segments.Proof.
The proof is by induction inj and t. For the basis with j=1, F(j,t)=∞ is the optimal cost since the approximation of one point is undefined for all t, and, with t=1, F(j,t)=E(1,j) for all j>1. For the induction, take any j>1 and t>1. Any approximation of {p1,…,pj} with t segments must approximate {p1,…,pi} with t-1 segments and {pi,…,pj} with one segment. The point pi is chosen to minimize the aggregated cost F(i,t-1) and the cost E(i,j) to approximate {pi,…,pj} with one segment. Since F(i,t-1) is by induction the minimum cost to approximate {p1,…,pi} with t-1 segments, F(j,t) is the minimum cost.The optimal sequenceI=〈i(1),i(2),…,i(T+1)〉 is obtained recursively: i(T+1)=n; i(T) is the index i that minimizes F(i(T+1),T) in (2); i(T-1) is the index i that minimizes F(i(T),T-1); and so forth until F(i(2),1) is reached, at which point i(1)=1. Algorithm 1 outlines the DP algorithm PL-Fixed for computing F(j,t). The algorithm records in table I(i,t) the indices i that minimize (2), which are the optimal breakpoints.Algorithm 1:PL-Fixed(P,T).
(1)for  i=1 to n  do
(2)E(i,i)≔∞;
(3)for  j=i+1 to n  do
(4)ai,j≔(j-i+1)∑k=ik=j‍xkyk-(∑k=ik=j‍xk)(∑k=ik=j‍yk)(j-i+1)∑k=ik=j‍xk2-(∑k=ik=j‍xk)2;
(5)bi,j≔∑k=ik=j‍yk-ai,j∑k=ik=j‍xk(j-i+1);
(6)Ai,j≔ai,j; Bi,j≔bi,j;
(7)Ei,j≔∑k=ijyk-ai,jxk-bi,j2;
(8)for  j=1 to n  do
(9)for  t=j to T  do
(10)Fj,t≔∞; I(j,t)≔0;
(11)F(j,1)≔E(1,j); I(j,1)≔1;
(12)for  t=2 to min⁡{j-1,T}  do
(13)F(j,t)≔∞; I(j,t)≔0;
(14)for  i=t to j  do
(15)if  F(j,t)>F(i,t-1)+E(i,j)  then
(16)F(j,t)≔F(i,t-1)+E(i,j);
(17)I(j,t)≔i;PL-Fixed is correct because it solves recursive equation (2). It runs in Θ(n3) time since T<n. The notation Θ(h(n))={g(n):∃n0,c1,c2>0, such that c1h(n)≤g(n)≤c2h(n) for all n≥n0} represents the set of all nondecreasing functions g(n) that grow at the same asymptotic rate of h(n). Therefore, the time taken by algorithm PL-Fixed to solve an instance of size n is approximately n3. A precise and extensive discussion on the running-time complexity of algorithms is found in [11].Running the algorithm on the sample instance withT=3 yields the values F(i,j) and I(i,j) depicted in Table 2. The optimal sequence is recursively obtained from the table as follows: i(4)=8; i(3)=I(8,3)=5; i(2)=I(5,2)=2; and i(1)=I(2,1)=1, as indicated in boldface. The optimal sequence is 〈1,2,5,8〉 according to which the first segment approximates the points {p1,p2}, the second segment approximates {p2,…,p5}, and the third segment approximates {p5,…,p8}. The approximation cost is F(8,3)=1.6648×10-2. Figure 1 depicts this approximation.Table 2
Tables with valuesF(j,t) and I(j,t) for the sample instance.
j
∖
t
F
(
j
,
t
) : I(j,t)
1
2
3
1
∞ : 1
∞ : 0
∞ : 0
2
0 : 1
∞ : 0
∞ : 0
3
5.5459
e
-
03 : 1
0 : 2
∞ : 0
4
2.2611
e
-
02 : 1
1.7718
e
-
03 : 2
0 : 3
5
2.3016
e
-
02 : 1
5.5550
e
-
03 : 2
1.7718
e
-
03 : 4
6
3.5038
e
-
02 : 1
2.3016
e
-
02 : 5
4.0758
e
-
03 : 4
7
4.4083
e
-
02 : 1
2.4617
e
-
02 : 2
1.3596
e
-
02 : 5
8
1.2685
e
-
01 : 1
3.4109
e
-
02 : 5
1.6648
e
-
02 : 5Figure 1
Piecewise-linear approximation for the sample instance with 3 segments.


### 2.2.  Maximum Number of Linear Segments

A variation of the previous model is obtained by limiting rather than fixing the number of linear segments. Each segment incurs a costc to trade off model simplicity and quality of approximation. This problem is also cast in a recursive form, with F(j,t) being the minimum cost to approximate {p1,p2,…,pj} with at most t≥1 segments. Formally,(3)Fj,t=min1≤i≤jFi,t-1+Ei,j+c,where E(i,j) is defined as above, F(1,t)=∞, and F(2,t)=E(1,2)+c for all t.Theorem 2.
Recursive equation (3) yields the piecewise-linear model of P with at most t segments, whose objective is the combination of the approximation error and a cost c per segment.Recursive equation (3) leads to a DP algorithm PL-Max, which is omitted here in the interest of brevity. PL-Max(P,T,c) runs in Θ(n3) provided that T≤n. This algorithm records the index i corresponding to the optimal breakpoint of (3) in a table I(j,t). This table is key to determine the optimal number T⋆ of segments and construct the optimal sequence 〈i(1),i(2),…,i(T⋆+1)〉.The optimal number of segments and sequence of breakpoints are determined from the tables produced by PL-Max(P,T,c) as follows: i(T+1)=n; i(T)=I(n,T) is the index that minimizes F(n,T) in (3); i(T-1)=I(i(T),T-1) is the index that minimizes F(i(T),T-1); and so forth until reaching the first breakpoint where i(T+1-k)=1 for some k. Then, T⋆=k and the optimal breakpoint sequence is 〈i(T-k+1),i(T-k+2),…,i(T+1)〉.With a limitT=5 and cost c=0.01, PL-Max produces an optimal approximation that uses 4 out of the 5 segments available, inducing a total cost of F(8,5)=0.044208. The optimal sequence is 〈i(1),…,i(5)〉=〈1,2,4,6,8〉 with four segments approximating the point sets {p1,p2}, {p2,p3,p4}, {p4,p5,p6}, and {p6,p7,p8}.


## 2.1. Fixed Number of Linear Segments

To find the piecewise-linear approximation with a fixed numberT of linear segments that minimizes the squared error, the following problem is solved:(1)minIE=∑t=1TEt=∑t=1T∑k=itit+1yk-atxk-bt2.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 F(j,t) be the minimum cost to approximate the points {p1,p2,…,pj} with exactly t segments. Then,(2)Fj,t=mint≤i≤jFi,t-1+Ei,j,where E(i,j) is the minimum squared error of approximating Pi,j={pi,pi+1,…,pj} with just one segment. Notice that E(i,i)=∞, F(j,t)=∞ for all j≤t, and F(j,1)=E(1,j) for all j>1. From recursive equation (2), clearly F(n,T)=minE.Theorem 1.
Recursive equation (2) yields the piecewise-linear approximation of P of minimum squared error having exactly t linear segments.Proof.
The proof is by induction inj and t. For the basis with j=1, F(j,t)=∞ is the optimal cost since the approximation of one point is undefined for all t, and, with t=1, F(j,t)=E(1,j) for all j>1. For the induction, take any j>1 and t>1. Any approximation of {p1,…,pj} with t segments must approximate {p1,…,pi} with t-1 segments and {pi,…,pj} with one segment. The point pi is chosen to minimize the aggregated cost F(i,t-1) and the cost E(i,j) to approximate {pi,…,pj} with one segment. Since F(i,t-1) is by induction the minimum cost to approximate {p1,…,pi} with t-1 segments, F(j,t) is the minimum cost.The optimal sequenceI=〈i(1),i(2),…,i(T+1)〉 is obtained recursively: i(T+1)=n; i(T) is the index i that minimizes F(i(T+1),T) in (2); i(T-1) is the index i that minimizes F(i(T),T-1); and so forth until F(i(2),1) is reached, at which point i(1)=1. Algorithm 1 outlines the DP algorithm PL-Fixed for computing F(j,t). The algorithm records in table I(i,t) the indices i that minimize (2), which are the optimal breakpoints.Algorithm 1:PL-Fixed(P,T).
(1)for  i=1 to n  do
(2)E(i,i)≔∞;
(3)for  j=i+1 to n  do
(4)ai,j≔(j-i+1)∑k=ik=j‍xkyk-(∑k=ik=j‍xk)(∑k=ik=j‍yk)(j-i+1)∑k=ik=j‍xk2-(∑k=ik=j‍xk)2;
(5)bi,j≔∑k=ik=j‍yk-ai,j∑k=ik=j‍xk(j-i+1);
(6)Ai,j≔ai,j; Bi,j≔bi,j;
(7)Ei,j≔∑k=ijyk-ai,jxk-bi,j2;
(8)for  j=1 to n  do
(9)for  t=j to T  do
(10)Fj,t≔∞; I(j,t)≔0;
(11)F(j,1)≔E(1,j); I(j,1)≔1;
(12)for  t=2 to min⁡{j-1,T}  do
(13)F(j,t)≔∞; I(j,t)≔0;
(14)for  i=t to j  do
(15)if  F(j,t)>F(i,t-1)+E(i,j)  then
(16)F(j,t)≔F(i,t-1)+E(i,j);
(17)I(j,t)≔i;PL-Fixed is correct because it solves recursive equation (2). It runs in Θ(n3) time since T<n. The notation Θ(h(n))={g(n):∃n0,c1,c2>0, such that c1h(n)≤g(n)≤c2h(n) for all n≥n0} represents the set of all nondecreasing functions g(n) that grow at the same asymptotic rate of h(n). Therefore, the time taken by algorithm PL-Fixed to solve an instance of size n is approximately n3. A precise and extensive discussion on the running-time complexity of algorithms is found in [11].Running the algorithm on the sample instance withT=3 yields the values F(i,j) and I(i,j) depicted in Table 2. The optimal sequence is recursively obtained from the table as follows: i(4)=8; i(3)=I(8,3)=5; i(2)=I(5,2)=2; and i(1)=I(2,1)=1, as indicated in boldface. The optimal sequence is 〈1,2,5,8〉 according to which the first segment approximates the points {p1,p2}, the second segment approximates {p2,…,p5}, and the third segment approximates {p5,…,p8}. The approximation cost is F(8,3)=1.6648×10-2. Figure 1 depicts this approximation.Table 2
Tables with valuesF(j,t) and I(j,t) for the sample instance.
j
∖
t
F
(
j
,
t
) : I(j,t)
1
2
3
1
∞ : 1
∞ : 0
∞ : 0
2
0 : 1
∞ : 0
∞ : 0
3
5.5459
e
-
03 : 1
0 : 2
∞ : 0
4
2.2611
e
-
02 : 1
1.7718
e
-
03 : 2
0 : 3
5
2.3016
e
-
02 : 1
5.5550
e
-
03 : 2
1.7718
e
-
03 : 4
6
3.5038
e
-
02 : 1
2.3016
e
-
02 : 5
4.0758
e
-
03 : 4
7
4.4083
e
-
02 : 1
2.4617
e
-
02 : 2
1.3596
e
-
02 : 5
8
1.2685
e
-
01 : 1
3.4109
e
-
02 : 5
1.6648
e
-
02 : 5Figure 1
Piecewise-linear approximation for the sample instance with 3 segments.


## 2.2.  Maximum Number of Linear Segments

A variation of the previous model is obtained by limiting rather than fixing the number of linear segments. Each segment incurs a costc to trade off model simplicity and quality of approximation. This problem is also cast in a recursive form, with F(j,t) being the minimum cost to approximate {p1,p2,…,pj} with at most t≥1 segments. Formally,(3)Fj,t=min1≤i≤jFi,t-1+Ei,j+c,where E(i,j) is defined as above, F(1,t)=∞, and F(2,t)=E(1,2)+c for all t.Theorem 2.
Recursive equation (3) yields the piecewise-linear model of P with at most t segments, whose objective is the combination of the approximation error and a cost c per segment.Recursive equation (3) leads to a DP algorithm PL-Max, which is omitted here in the interest of brevity. PL-Max(P,T,c) runs in Θ(n3) provided that T≤n. This algorithm records the index i corresponding to the optimal breakpoint of (3) in a table I(j,t). This table is key to determine the optimal number T⋆ of segments and construct the optimal sequence 〈i(1),i(2),…,i(T⋆+1)〉.The optimal number of segments and sequence of breakpoints are determined from the tables produced by PL-Max(P,T,c) as follows: i(T+1)=n; i(T)=I(n,T) is the index that minimizes F(n,T) in (3); i(T-1)=I(i(T),T-1) is the index that minimizes F(i(T),T-1); and so forth until reaching the first breakpoint where i(T+1-k)=1 for some k. Then, T⋆=k and the optimal breakpoint sequence is 〈i(T-k+1),i(T-k+2),…,i(T+1)〉.With a limitT=5 and cost c=0.01, PL-Max produces an optimal approximation that uses 4 out of the 5 segments available, inducing a total cost of F(8,5)=0.044208. The optimal sequence is 〈i(1),…,i(5)〉=〈1,2,4,6,8〉 with four segments approximating the point sets {p1,p2}, {p2,p3,p4}, {p4,p5,p6}, and {p6,p7,p8}.


## 3. Continuous Approximations

The preceding models may produce discontinuous piecewise-linear 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.


### 3.1. Ensuring Continuity at the Breakpoints

Continuity is enforced by approximating the setPi,j={pi,pi+1,…,pj} with segments going through pi and pj. The recursive formulations ensuring continuity with a fixed and maximum number of segments are identical to the preceding ones, except that Pi,j is approximated with the segment passing through pi and pj. The theorem below is a consequence of this observation.Theorem 3.
Suppose that the coefficients of the linear segment approximatingPi,j are ai,j=(yj-yi)/(xj-xi) and bi,j=yj-ai,jxj. Then, recursive formulations (2) and (3) produce optimal piecewise-linear approximation functions for P that are continuous at the breakpoints using exactly and at most t linear segments, respectively.Applying PL-Max to the sample instance with a maximum ofT=5 segments and cost c=0.01, but forcing the segments approximating Pi,j to go through the end points, yields the piecewise-linear function with T⋆=4 segments, which is illustrated in Figure 2. This approximation is given by the sequence 〈i(1),…,i(5)〉=〈1,2,4,6,8〉 which, coincidentally, is identical to the one obtained without ensuring continuity. However, its approximation cost F(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.Figure 2
Continuous piecewise-linear approximation of the sample instance with a maximum of 5 segments.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.


### 3.2.  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 numbert>1 of segments is(4)Fj,t=mini∈Sj,tFi,t-1+Ei,j,where Sj,t={i:t≤i≤j, the lines approximating P(t-1) and P(t) intercept at point x∈[xmax{1,i-1},xmin{n,i+1}]}, and F(j,t)=∞ if Sj,t=∅. Further, E(i,i)=∞, F(j,t)=∞ for all j≤t, and F(j,1)=E(1,j).Algorithm2 gives the dynamic programming solution of recursive equation (4), for a fixed number T of segments. The algorithm records in table X(j,t) the points where the linear segments intercept, which are used to determine the breakpoints of the piecewise-linear function. The algorithm runs in O(n3) time provided that T=O(n).Algorithm 2:PL-Fixed-Cont-Reject(P,T).
(1)  Same as lines (1) through (7) of PL-Fixed
(2)for  j=1 to n  do
(3)for  t=j to T  do
(4)F(j,t)≔∞; I(j,t)≔0; X(j,t)≔0;
(5)F(j,1)≔E(1,j); I(j,1):=1; X(j,1)≔0;
(6)for  t=2 to min⁡{j-1,T}  do
(7)F(j,t)≔∞; I(j,t)≔0; X(j,1)≔0;
(8)for  i=t to j-1  do
(9)k≔I(i,t-1);
(10)if  k≠0 and Ak,i≠Ai,j  then
(11)x≔(B(k,i)-B(i,j))(A(i,j)-A(k,i));
(12)if  xmax⁡{1,i-1}≤x≤xmin⁡{n,i+1} and F(j,t)>F(i,t-1)+E(i,j)  then
(13)F(j,t)≔F(i,t-1)+E(i,j);  I(j,t)≔i;  X(j,t)≔x;Algorithm2 applied to the sample instance with T=3 yields 〈i(1),…,i(4)〉=〈1,5,6,8〉 as the optimal sequence. However, the breakpoints are not those of the sequence but defined by the intercept of the segments: the first segment is defined for [x1,0.3120]=[0.0114,0.3120], the second for [0.3120,0.3216], and the third for [0.3216,x8]=[0.3216,0.9104]. Notice that x4≤0.3112≤x6 and x5≤0.3216≤x7 as expected. The approximation cost is F(8,3)=2.3148×10-2, which exceeds the cost 1.6648×10-2 of the discontinuous approximation.


### 3.3.  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 segmentst>1, recursive equation (2) becomes(5)Fj,t=mint≤i≤jFi,t-1+Ei,j⋆,while for a maximum number of segments (3) becomes(6)Fj,t=min1≤i≤jFi,t-1+Ei,j⋆+c,where E(i,j)⋆ is obtained by solving the following problem:(7a)Ei,j⋆=minai,j,bi,j,x∑k=ijai,jxk+bi,j-yk2(7b)s.t.:Ai,t-1x+Bi,t-1=ai,jx+bi,j(7c)xmax1,i-1≤x≤xmin⁡i+1,n,where Ai,t-1 and Bi,t-1 are the slope and intercept of the rightmost segment of the optimal approximation of P1,i with exactly or at most t-1 segments, depending on which recursion is being solved. E(i,j)⋆=∞ if the problem is infeasible, for example, when j=i.The computation ofE(i,j)⋆ entails minimizing the squared error of the approximation of {pi,pi+1,…,pj} with a linear segment, while ensuring that this segment intercepts the preceding one, the rightmost segment of the piecewise-linear approximation of {p1,…,pi}, in the interval [xmax{1,i-1},xmin{i+1,n}]. The chief difficulty lies in the nonlinearity of (7b). However, the optimal solution is obtained by solving two quadratic programs. First, notice that variable x is eliminated by replacing the constraints (7b)-(7c) with the following nonlinear inequality:(8)xmax1,i-1≤bi,j-Bi,t-1Ai,t-1-ai,j≤xmini+1,nwhich, for ai,j<Ai,t-1, is transformed into the linear inequalities (9a)xmax1,i-1ai,j+bi,j≥Bi,t-1+xmax⁡1,i-1Ai,t-1,xmini+1,nai,j+bi,j≤Bi,t-1+xmin⁡i+1,nAi,t-1,ai,j<Ai,t-1and the following ones when ai,j>Ai,t-1:(9b)xmax1,i-1ai,j+bi,j≤Bi,t-1+xmax⁡1,i-1Ai,t-1,xmini+1,nai,j+bi,j≥Bi,t-1+xmin⁡i+1,nAi,t-1,ai,j>Ai,t-1.The end result is that E(i,j)⋆ 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(n3g) time, with g 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 T=3 segments yields the sequence 〈i(1),…,i(4)〉=〈1,5,6,8〉. The breakpoints of the segments are defined by their intercepts, which are [x1,0.3120]=[0.0114,0.3120] for the first segment, [0.3120,0.3215] for the second, and finally [0.3215,x8]=[0.3215,0.9104] for the third. The approximation cost is F(8,3)=2.3147×10-2.


## 3.1. Ensuring Continuity at the Breakpoints

Continuity is enforced by approximating the setPi,j={pi,pi+1,…,pj} with segments going through pi and pj. The recursive formulations ensuring continuity with a fixed and maximum number of segments are identical to the preceding ones, except that Pi,j is approximated with the segment passing through pi and pj. The theorem below is a consequence of this observation.Theorem 3.
Suppose that the coefficients of the linear segment approximatingPi,j are ai,j=(yj-yi)/(xj-xi) and bi,j=yj-ai,jxj. Then, recursive formulations (2) and (3) produce optimal piecewise-linear approximation functions for P that are continuous at the breakpoints using exactly and at most t linear segments, respectively.Applying PL-Max to the sample instance with a maximum ofT=5 segments and cost c=0.01, but forcing the segments approximating Pi,j to go through the end points, yields the piecewise-linear function with T⋆=4 segments, which is illustrated in Figure 2. This approximation is given by the sequence 〈i(1),…,i(5)〉=〈1,2,4,6,8〉 which, coincidentally, is identical to the one obtained without ensuring continuity. However, its approximation cost F(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.Figure 2
Continuous piecewise-linear approximation of the sample instance with a maximum of 5 segments.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.


## 3.2.  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 numbert>1 of segments is(4)Fj,t=mini∈Sj,tFi,t-1+Ei,j,where Sj,t={i:t≤i≤j, the lines approximating P(t-1) and P(t) intercept at point x∈[xmax{1,i-1},xmin{n,i+1}]}, and F(j,t)=∞ if Sj,t=∅. Further, E(i,i)=∞, F(j,t)=∞ for all j≤t, and F(j,1)=E(1,j).Algorithm2 gives the dynamic programming solution of recursive equation (4), for a fixed number T of segments. The algorithm records in table X(j,t) the points where the linear segments intercept, which are used to determine the breakpoints of the piecewise-linear function. The algorithm runs in O(n3) time provided that T=O(n).Algorithm 2:PL-Fixed-Cont-Reject(P,T).
(1)  Same as lines (1) through (7) of PL-Fixed
(2)for  j=1 to n  do
(3)for  t=j to T  do
(4)F(j,t)≔∞; I(j,t)≔0; X(j,t)≔0;
(5)F(j,1)≔E(1,j); I(j,1):=1; X(j,1)≔0;
(6)for  t=2 to min⁡{j-1,T}  do
(7)F(j,t)≔∞; I(j,t)≔0; X(j,1)≔0;
(8)for  i=t to j-1  do
(9)k≔I(i,t-1);
(10)if  k≠0 and Ak,i≠Ai,j  then
(11)x≔(B(k,i)-B(i,j))(A(i,j)-A(k,i));
(12)if  xmax⁡{1,i-1}≤x≤xmin⁡{n,i+1} and F(j,t)>F(i,t-1)+E(i,j)  then
(13)F(j,t)≔F(i,t-1)+E(i,j);  I(j,t)≔i;  X(j,t)≔x;Algorithm2 applied to the sample instance with T=3 yields 〈i(1),…,i(4)〉=〈1,5,6,8〉 as the optimal sequence. However, the breakpoints are not those of the sequence but defined by the intercept of the segments: the first segment is defined for [x1,0.3120]=[0.0114,0.3120], the second for [0.3120,0.3216], and the third for [0.3216,x8]=[0.3216,0.9104]. Notice that x4≤0.3112≤x6 and x5≤0.3216≤x7 as expected. The approximation cost is F(8,3)=2.3148×10-2, which exceeds the cost 1.6648×10-2 of the discontinuous approximation.


## 3.3.  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 segmentst>1, recursive equation (2) becomes(5)Fj,t=mint≤i≤jFi,t-1+Ei,j⋆,while for a maximum number of segments (3) becomes(6)Fj,t=min1≤i≤jFi,t-1+Ei,j⋆+c,where E(i,j)⋆ is obtained by solving the following problem:(7a)Ei,j⋆=minai,j,bi,j,x∑k=ijai,jxk+bi,j-yk2(7b)s.t.:Ai,t-1x+Bi,t-1=ai,jx+bi,j(7c)xmax1,i-1≤x≤xmin⁡i+1,n,where Ai,t-1 and Bi,t-1 are the slope and intercept of the rightmost segment of the optimal approximation of P1,i with exactly or at most t-1 segments, depending on which recursion is being solved. E(i,j)⋆=∞ if the problem is infeasible, for example, when j=i.The computation ofE(i,j)⋆ entails minimizing the squared error of the approximation of {pi,pi+1,…,pj} with a linear segment, while ensuring that this segment intercepts the preceding one, the rightmost segment of the piecewise-linear approximation of {p1,…,pi}, in the interval [xmax{1,i-1},xmin{i+1,n}]. The chief difficulty lies in the nonlinearity of (7b). However, the optimal solution is obtained by solving two quadratic programs. First, notice that variable x is eliminated by replacing the constraints (7b)-(7c) with the following nonlinear inequality:(8)xmax1,i-1≤bi,j-Bi,t-1Ai,t-1-ai,j≤xmini+1,nwhich, for ai,j<Ai,t-1, is transformed into the linear inequalities (9a)xmax1,i-1ai,j+bi,j≥Bi,t-1+xmax⁡1,i-1Ai,t-1,xmini+1,nai,j+bi,j≤Bi,t-1+xmin⁡i+1,nAi,t-1,ai,j<Ai,t-1and the following ones when ai,j>Ai,t-1:(9b)xmax1,i-1ai,j+bi,j≤Bi,t-1+xmax⁡1,i-1Ai,t-1,xmini+1,nai,j+bi,j≥Bi,t-1+xmin⁡i+1,nAi,t-1,ai,j>Ai,t-1.The end result is that E(i,j)⋆ 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(n3g) time, with g 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 T=3 segments yields the sequence 〈i(1),…,i(4)〉=〈1,5,6,8〉. The breakpoints of the segments are defined by their intercepts, which are [x1,0.3120]=[0.0114,0.3120] for the first segment, [0.3120,0.3215] for the second, and finally [0.3215,x8]=[0.3215,0.9104] for the third. The approximation cost is F(8,3)=2.3147×10-2.


## 4. Lower and Upper Envelopes

The preceding models can be modified to obtain piecewise-linear lower and upper envelopes of the generating functionΦ. For the upper envelope, the slope ai,j and intercept bi,j should minimize the squared error of the line y=ai,jx+bi,j approximating Pi,j={pi,pi+1,…,pj}, while keeping these points below the line. Formally, this is achieved by solving the QP: (10a)Li,j:minai,j,bi,j∑k=ijai,jxk+bi,j-yk2(10b)s.t.:yk≤ai,jxk+bi,j,pk∈Pi,j.The basic algorithms can compute the optimal upper envelope, with a fixed and maximum number of segments, by solvingLi,j to obtain the coefficients of the segment approximating Pi,j. 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 P with exactly and at most t linear segments, respectively, having the squared approximation error as the objective, provided that the coefficients for the linear segments are computed by solving problem Li,j given by (10a) and (10b).The running time of the resulting algorithm isO(n3g(n)) where g(n) is an upper bound on the cost to solve the quadratic program Li,j with 2 variables and n constraints.Algorithm PL-Fixed applied to the sample instance, with the coefficients obtained by solvingLi,j and a fixed number of segments T=3, yields the optimal sequence 〈i(1),…,i(T+1)〉=〈1,2,5,8〉 with cost F(8,3)=2.8762×10-2. The slopes and intercepts of the segments are (aij,bij)=(3.5413,0.1453) for (1,2); (1.0188,0.3830) for (2,5); and (0.2225,0.7623) for (5,8).The models and algorithms presented in Section3 can be modified to produce continuous lower and upper envelopes. For PL-Fixed-Cont-Reject, it suffices to use the slope ai,j and intercept bi,j obtained by solving Li,j. This algorithm applied to the sample instance with T=3 yields the optimal sequence 〈i(1),…,i(T+1)〉=〈1,5,6,8〉 with cost F(8,3)=4.609×10-2. The domains of the linear segments are [x1,0.3153]=[0.0114,0.3153] for the first segment, [0.3153,0.3222] for the second, and [0.3222,x8]=[0.3222,0.9104] for the third. This approximation is illustrated in Figure 3.Figure 3
Continuous, upper envelope approximation with 3 linear segments obtained with the rejection strategy.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.


## 5. Dealing with Random Data

Uncertainty is addressed by modeling each pointpi∈P as a random variable pi=(xi,yi) characterized by a probability density function fi(x,y). Two models are proposed, one minimizing expected squared error and the other using chance constraints [13].


### 5.1. Minimizing Expected Squared Error

Uncertainty may be treated by minimizing the expected squared error of the piecewise-linear approximation. For the points in the setPi,j={pi,pi+1,…,pj}, this approach finds a linear segment given by the slope ai,j and intercept bi,j that solves the problem: (11a)Fi,j:fi,j=minai,j,bi,j∑k=ijfi,jkai,j,bi,j,where(11b)fi,jkai,j,bi,j=∫x∈Xk∫y∈Ykai,jx+bi,j-y2fkx,ydxdy,(Xk,Yk) is the domain of pk=(xk,yk), fk is the joint probability density function for (xk,yk), and fi,jk is the expected error incurred by pk when approximating Pi,j with a linear segment.The functionfi,jk(ai,j,bi,j) is convex on the variables (ai,j,bi,j) since the term (ai,jx+bi,j-y)2 is convex for fixed (x,y), generalizing an infinite weighted sum [14, Sec. 3.2.1]. This means that fi,j is a convex function and so is problem Fi,j. Hence, an optimal solution to Fi,j 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 Section2 yield the optimal approximation, for a fixed and maximum number of segments, by solving Fi,j and using the obtained coefficients for the linear segment approximating Pi,j. 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={p1,…,pn} of random points with exactly and at most t segments, respectively, having the expected squared approximation error as the objective if the coefficients for the linear segments are obtained by solving Fi,j given by (11a) and (11b).The resulting algorithms run inO(n3g(n)) time, with g(n) being an upper bound on the running time to solve Fi,j with 2 variables and n integration terms.A sample instance is obtained by definingxi and yi as independent random variables, both characterized by Gaussian distributions with mean μix=xi and μiy=yi and standard deviations σix=xi/2 and σiy=yi/2. Table 3 shows the statistics of the random points of the stochastic sample instance. Thus, the joint probability density function is (12)fixi,yi=e-1/2xi-μix/σix2σix2π×e-1/2yi-μiy/σiy2σiy2π.Table 3
Statistics of the random points.
i
μ
i
x
σ
i
x
μ
i
y
σ
i
y
1
0.0114
0.0057
0.1857
0.0929
2
0.0861
0.0431
0.4502
0.2251
3
0.1721
0.0861
0.5584
0.2792
4
0.2928
0.1464
0.5857
0.2929
5
0.3130
0.1565
0.7020
0.3510
6
0.3212
0.1606
0.8338
0.4169
7
0.4899
0.2450
0.8713
0.4357
8
0.8807
0.4404
0.9104
0.4552With a fixed numberT=3 of segments and E(i,j)=f(i,j), PL-Fixed produced the optimal sequence 〈i(1),…,i(T+1)〉=〈1,2,4,8〉 with an expected approximation cost F(8,3)=1.1311. This piecewise-linear function differs from the one with zero variance, for which the optimal sequence is 〈i(1),…,i(T+1)〉=〈1,2,5,8〉 with cost F(8,3)=1.6648×10-2.


### 5.2. 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 problemFi,j. An algorithm for the lower envelope is developed below, which is easily adapted for the upper envelope.To find the lower envelope ofPi,j={pi,pi+1,…,pj}, the constraint that the random point pk=(xk,yk) is above the linear segment ai,jx+bi,j, i≤k≤j, is imposed with a given probability α. Mathematically,(13)Pai,jx+bi,j≤yk≥α,x∈Xk,where P(·) denotes the probability operator. The DP algorithms can find the lower envelope by changing how the slope ai,j and intercept bi,j 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 (14a)fi,jle=minai,j,bi,j∑k=ijfi,jkai,j,bi,jwhile being subject to(14b)fi,jkai,j,bi,j=∫x∈Xk∫y∈Ykai,jx+bi,j-y2fkx,ydxdy,(14c)Pai,jx+bi,j≤yk≥α,x∈Xk,k=i,…,j.By expressing probabilistic constraint (14c) in terms of integrals and the joint probability density function fk, a deterministic equivalent to this problem is obtained: (15a)Fi,jle:fi,jle=minai,j,bi,j∑k=ijfi,jkai,j,bi,jwhile being subject to(15b)fi,jkai,j,bi,j=∫x∈Xk∫y∈Ykai,jx+bi,j-y2fkx,ydxdy,gi,jkai,j,bi,j=∫x∈Xk∫y∈Yk:y≤ai,jx+bi,jfkx,ydxdy≤1-α,k=i,…,j,where gi,jk is the probability of pk lying below the linear segment that approximates Pi,j. This problem may be posed in terms of the probability g¯i,jk of pk being above the linear segment.ProblemFi,jle 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 solvingFi,jle and using the obtained coefficients in the approximation of Pi,j. 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={p1,…,pn} of random points with exactly and at most t 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 Fi,jle.Algorithm PL-Fixed was applied to the stochastic sample instance given in Table3, with a fixed number T=3 of segments and probabilities α=60% and 70%. The problems Fi,jle were solved approximately using an SQP algorithm. For α=60%, the algorithm produced the sequence 〈i(1),…,i(T+1)〉=〈1,2,5,8〉 for the probabilistic lower envelope. The expected cost of this solution is F(8,3)=1.3280, above the expected cost of 1.1311 obtained for the solution without chance constraints. For α=70%, the algorithm produced the sequence 〈1,3,5,8〉 with an expected cost of F(8,3)=1.6369.


## 5.1. Minimizing Expected Squared Error

Uncertainty may be treated by minimizing the expected squared error of the piecewise-linear approximation. For the points in the setPi,j={pi,pi+1,…,pj}, this approach finds a linear segment given by the slope ai,j and intercept bi,j that solves the problem: (11a)Fi,j:fi,j=minai,j,bi,j∑k=ijfi,jkai,j,bi,j,where(11b)fi,jkai,j,bi,j=∫x∈Xk∫y∈Ykai,jx+bi,j-y2fkx,ydxdy,(Xk,Yk) is the domain of pk=(xk,yk), fk is the joint probability density function for (xk,yk), and fi,jk is the expected error incurred by pk when approximating Pi,j with a linear segment.The functionfi,jk(ai,j,bi,j) is convex on the variables (ai,j,bi,j) since the term (ai,jx+bi,j-y)2 is convex for fixed (x,y), generalizing an infinite weighted sum [14, Sec. 3.2.1]. This means that fi,j is a convex function and so is problem Fi,j. Hence, an optimal solution to Fi,j 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 Section2 yield the optimal approximation, for a fixed and maximum number of segments, by solving Fi,j and using the obtained coefficients for the linear segment approximating Pi,j. 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={p1,…,pn} of random points with exactly and at most t segments, respectively, having the expected squared approximation error as the objective if the coefficients for the linear segments are obtained by solving Fi,j given by (11a) and (11b).The resulting algorithms run inO(n3g(n)) time, with g(n) being an upper bound on the running time to solve Fi,j with 2 variables and n integration terms.A sample instance is obtained by definingxi and yi as independent random variables, both characterized by Gaussian distributions with mean μix=xi and μiy=yi and standard deviations σix=xi/2 and σiy=yi/2. Table 3 shows the statistics of the random points of the stochastic sample instance. Thus, the joint probability density function is (12)fixi,yi=e-1/2xi-μix/σix2σix2π×e-1/2yi-μiy/σiy2σiy2π.Table 3
Statistics of the random points.
i
μ
i
x
σ
i
x
μ
i
y
σ
i
y
1
0.0114
0.0057
0.1857
0.0929
2
0.0861
0.0431
0.4502
0.2251
3
0.1721
0.0861
0.5584
0.2792
4
0.2928
0.1464
0.5857
0.2929
5
0.3130
0.1565
0.7020
0.3510
6
0.3212
0.1606
0.8338
0.4169
7
0.4899
0.2450
0.8713
0.4357
8
0.8807
0.4404
0.9104
0.4552With a fixed numberT=3 of segments and E(i,j)=f(i,j), PL-Fixed produced the optimal sequence 〈i(1),…,i(T+1)〉=〈1,2,4,8〉 with an expected approximation cost F(8,3)=1.1311. This piecewise-linear function differs from the one with zero variance, for which the optimal sequence is 〈i(1),…,i(T+1)〉=〈1,2,5,8〉 with cost F(8,3)=1.6648×10-2.


## 5.2. 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 problemFi,j. An algorithm for the lower envelope is developed below, which is easily adapted for the upper envelope.To find the lower envelope ofPi,j={pi,pi+1,…,pj}, the constraint that the random point pk=(xk,yk) is above the linear segment ai,jx+bi,j, i≤k≤j, is imposed with a given probability α. Mathematically,(13)Pai,jx+bi,j≤yk≥α,x∈Xk,where P(·) denotes the probability operator. The DP algorithms can find the lower envelope by changing how the slope ai,j and intercept bi,j 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 (14a)fi,jle=minai,j,bi,j∑k=ijfi,jkai,j,bi,jwhile being subject to(14b)fi,jkai,j,bi,j=∫x∈Xk∫y∈Ykai,jx+bi,j-y2fkx,ydxdy,(14c)Pai,jx+bi,j≤yk≥α,x∈Xk,k=i,…,j.By expressing probabilistic constraint (14c) in terms of integrals and the joint probability density function fk, a deterministic equivalent to this problem is obtained: (15a)Fi,jle:fi,jle=minai,j,bi,j∑k=ijfi,jkai,j,bi,jwhile being subject to(15b)fi,jkai,j,bi,j=∫x∈Xk∫y∈Ykai,jx+bi,j-y2fkx,ydxdy,gi,jkai,j,bi,j=∫x∈Xk∫y∈Yk:y≤ai,jx+bi,jfkx,ydxdy≤1-α,k=i,…,j,where gi,jk is the probability of pk lying below the linear segment that approximates Pi,j. This problem may be posed in terms of the probability g¯i,jk of pk being above the linear segment.ProblemFi,jle 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 solvingFi,jle and using the obtained coefficients in the approximation of Pi,j. 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={p1,…,pn} of random points with exactly and at most t 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 Fi,jle.Algorithm PL-Fixed was applied to the stochastic sample instance given in Table3, with a fixed number T=3 of segments and probabilities α=60% and 70%. The problems Fi,jle were solved approximately using an SQP algorithm. For α=60%, the algorithm produced the sequence 〈i(1),…,i(T+1)〉=〈1,2,5,8〉 for the probabilistic lower envelope. The expected cost of this solution is F(8,3)=1.3280, above the expected cost of 1.1311 obtained for the solution without chance constraints. For α=70%, the algorithm produced the sequence 〈1,3,5,8〉 with an expected cost of F(8,3)=1.6369.


## 6. Computational Analysis

In order to keep the analysis brief, the experiments were restricted to the following approaches:(1)
PL-Max, the basic model of Section2.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 Section3.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 qi(qinji) of oil, water, and gas given as a function of the lift-gas rate qinji injected in an oil well i. The relation qi(qinji) 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 ofqi for a sufficiently high number of qinji values. In our experiments, 50 breakpoints produced a highly accurate and detailed description of qi. However, simpler and sufficiently accurate piecewise-linear 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.Table4 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 T=10 and the cost per segment to c=0. For PL-Max-Bound, T=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.Table 4
Running time of the algorithms (in seconds).
Well
Max
Max-Cont
Max-Bound
1
0.004
1.120
0.148
2
0.008
1.172
0.144
3
0.004
1.156
0.164
4
0.004
1.152
0.168
5
0.008
1.140
0.152
6
0.008
1.272
0.164
7
0.008
1.268
0.188
8
0.008
1.160
0.160
9
0.008
1.268
0.184
10
0.004
1.076
0.160The behavior of the PL-Max strategy is shown in Table5 for a varying cost c and a maximum number of segments T=15. The algorithm was applied to the production function of well #1.Table 5
Behavior of the PL-Max strategy.
c
Segments
Cost
0
15
29.42
10
15
179.42
40
12
565.28
50
11
677.97
100
10
1 197.42
200
8
2 113.22
1000
6
7 639.54The squared error grows asc increases, corroborating the theoretical behavior.The behavior of PL-Max-Cont strategy is illustrated in Table6 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 T increases.Table 6
Behavior of the PL-Max-Cont strategy.
T
Cost
1
784 752.00
5
3 346.49
10
197.42
15
29.42
20
6.44
25
1.68
30
0.47The computational experiments show that piecewise-linear 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].


## 7. 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 piecewise-linear 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 setP={p1,…,pn} 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 piecewise-linear 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 setP.(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 thel2-norm; however other norms can be easily adopted such as the l1- and l∞-norm by solving suitable linear and quadratic programs.



---

*Source: 876862-2015-07-07.xml*