! Fit a piecewise linear function of one variable
  to an arbitrary function, y = f(x).
  The objective is: minimize weighted sum of maximum error plus average error.
  The original function is presented as a sequence of (x,y) pairs.
  Having a piecewise linear approximation is useful if you want to model
  the problem with a Mixed Integer Linear Model, e.g., using SOS2 sets;
! Keywords: Piecewise linear, SOS2, Approximation, Linear approximation,
   Graph, Chart;
SETS:
    SAMPLE: X, Y, ERRUP, ERRDN, ZX, PCONS, PSLOPE, YFIT;
ENDSETS
DATA: NSEG = 7; ! Number line segments in the fit; YLB = -999999; ! Lower bound on y fitted, e.g., 0 if should have y >= 0; ERRUPUB = 999999; ! Upper bound on up error, = 0 to give a lower envelope; ERRDNUB = 999999; ! Upper bound on down error, = 0 to give an upper envelope; ALPHA = .05; ! Weight, in [0, 1], applied to average abs(error) vs. max abs(error); TIMEUB = 300; ! Time limit for guaranteed optimal; ! The points must be strictly increasing in X; ! Case 1, a function with an interior min and max; ! X, Y = -1.5 1.5625 -1.425 1.062187891 -1.35 0.67650625 -1.275 0.391406641 -1.2 0.1936 -1.125 0.070556641 -1.05 0.01050625 -0.975 0.002437891 -0.9 0.0361 -0.825 0.102000391 -0.75 0.19140625 -0.675 0.296344141 -0.6 0.4096 -0.525 0.524719141 -0.45 0.63600625 -0.375 0.738525391 -0.3 0.8281 -0.225 0.901312891 -0.15 0.95550625 -0.075 0.988781641 0 1 0.075 0.988781641 0.15 0.95550625 0.225 0.901312891 0.3 0.8281 0.375 0.738525391 0.45 0.63600625 0.525 0.524719141 0.6 0.4096 0.675 0.296344141 0.75 0.19140625 0.825 0.102000391 0.9 0.0361 0.975 0.002437891 1.05 0.01050625 1.125 0.070556641 1.2 0.1936 1.275 0.391406641 1.35 0.67650625 1.425 1.062187891 1.5 1.5625 1.575 2.192250391; ! Case 2; ! X, Y = 0 1.6 1 0.3 2 0.0 3 0.2 4 0.7 5 0.9 6 1.0 7 0.9 8 0.6 9 0.2 10 0.0; ! Case 3: A decline curve, y = (a - x)/(b + c*x); ! X, Y = 0.025 13 0.05 9.5 0.075 7.4 0.1 6 0.125 5 0.15 4.25 0.175 3.666666666666667 0.2 3.2 0.225 2.818181818181818 0.25 2.5 0.275 2.230769230769231 0.3 2 0.325 1.8 0.35 1.625 0.375 1.470588235294118 0.4 1.333333333333333 0.425 1.210526315789474 0.45 1.1 0.475 1 0.5 0.9090909090909091 0.525 0.8260869565217391 0.55 0.75 0.575 0.6800000000000001 0.6 0.6153846153846154 0.625 0.5555555555555556 0.65 0.5 0.675 0.4482758620689655 0.7 0.4 0.725 0.3548387096774194 0.75 0.3125 0.775 0.2727272727272727 0.8 0.2352941176470588 0.825 0.2 0.85 0.1666666666666667 0.875 0.1351351351351351 0.9 0.1052631578947368 0.925 0.07692307692307693 0.95 0.05 0.975 0.02439024390243903 1 0 ; ! SINE curve in[0, PI]; X Y = 0.0000000 0.0000000 0.0981748 0.0980171 0.1963495 0.1950903 0.2945243 0.2902847 0.3926991 0.3826834 0.4908739 0.4713967 0.5890486 0.5555702 0.6872234 0.6343933 0.7853982 0.7071068 0.8835729 0.7730105 0.9817477 0.8314696 1.0799225 0.8819213 1.1780972 0.9238795 1.2762720 0.9569403 1.3744468 0.9807853 1.4726216 0.9951847 1.5707963 1.0000000 1.6689711 0.9951847 1.7671459 0.9807853 1.8653206 0.9569403 1.9634954 0.9238795 2.0616702 0.8819213 2.1598449 0.8314696 2.2580197 0.7730105 2.3561945 0.7071068 2.4543693 0.6343933 2.5525440 0.5555702 2.6507188 0.4713967 2.7488936 0.3826834 2.8470683 0.2902847 2.9452431 0.1950903 3.0434179 0.0980171 3.1415927 0.0000000 ; ! Case: SINE curve in[0, 2*PI]; ! X Y = 0.0000000 0.0000000 0.0981748 0.0980171 0.1963495 0.1950903 0.2945243 0.2902847 0.3926991 0.3826834 0.4908739 0.4713967 0.5890486 0.5555702 0.6872234 0.6343933 0.7853982 0.7071068 0.8835729 0.7730105 0.9817477 0.8314696 1.0799225 0.8819213 1.1780972 0.9238795 1.2762720 0.9569403 1.3744468 0.9807853 1.4726216 0.9951847 1.5707963 1.0000000 1.6689711 0.9951847 1.7671459 0.9807853 1.8653206 0.9569403 1.9634954 0.9238795 2.0616702 0.8819213 2.1598449 0.8314696 2.2580197 0.7730105 2.3561945 0.7071068 2.4543693 0.6343933 2.5525440 0.5555702 2.6507188 0.4713967 2.7488936 0.3826834 2.8470683 0.2902847 2.9452431 0.1950903 3.0434179 0.0980171 3.1415927 0.0000000 3.2397674 -0.0980171 3.3379422 -0.1950903 3.4361170 -0.2902847 3.5342917 -0.3826834 3.6324665 -0.4713967 3.7306413 -0.5555702 3.8288160 -0.6343933 3.9269908 -0.7071068 4.0251656 -0.7730105 4.1233404 -0.8314696 4.2215151 -0.8819213 4.3196899 -0.9238795 4.4178647 -0.9569403 4.5160394 -0.9807853 4.6142142 -0.9951847 4.7123890 -1.0000000 4.8105638 -0.9951847 4.9087385 -0.9807853 5.0069133 -0.9569403 5.1050881 -0.9238795 5.2032628 -0.8819213 5.3014376 -0.8314696 5.3996124 -0.7730105 5.4977871 -0.7071068 5.5959619 -0.6343933 5.6941367 -0.5555702 5.7923115 -0.4713967 5.8904862 -0.3826834 5.9886610 -0.2902847 6.0868358 -0.1950903 6.1850105 -0.0980171 6.2831853 -0.0000000 ; ENDDATA SUBMODEL FITIT: ! Decision variables: PSLOPE(i) = slope of line for point i to i+1, PCONS(i) = constant term of line for point i to i+1, ZX(i) = 1 if there is a breakpoint at i; NP = @SIZE(SAMPLE); ! Points in representation of original function; MIN = ALPHA*AVGERR ! Min weighted combination of total and max error; + (1-ALPHA)*MAXERR; AVGERR = ( (ERRUP(1) + ERRDN(1))*(X(2) - X(1)) + (ERRUP(NP) + ERRDN(NP))*(X(NP) - X(NP-1)) + @SUM( SAMPLE(i) | 1 #LT# i #AND# i #LT# NP: ( ERRUP(i) + ERRDN(i))*(X(i+1) - X(i-1))) )/(2*(X(NP) - X(1))); ! Limit on number of nontrivial transitions/breakpoints; [MXBP] @SUM( SAMPLE(i): ZX(i)) <= NSEG - 1; @FOR(SAMPLE(I): @FREE( PCONS(i)); @FREE( PSLOPE(i)); @FREE( YFIT(i)); @BIN( ZX(i)); ! Fitted or predicted value at point i from piecewise linear (PL) approximation; YFIT(I) = PCONS(i) + PSLOPE(i)*X(i); YFIT( i) >= ylb; ! Allow a lower bound, e.g., 0, on the fitted value; ! Error at point i; ERRUP(I) - ERRDN(i) = YFIT(i) - Y(I); MAXERR >= ERRUP(I) + ERRDN(i); ! Max error at eany point; ! Upper bounds on the types of errors; ERRUP(i) <= ERRUPUB; ERRDN(i) <= ERRDNUB; ); @FOR( SAMPLE(i) | i #GT# 1: ! The PL curve must be continuous at i; YFIT( i) = PCONS(i-1) + PSLOPE(i-1)*X(i); ! Slope from i-1 to i must equal slope from i to i + 1 if no break at i; PSLOPE(i) >= PSLOPE(i-1) - ML*ZX(i); PSLOPE(i) <= PSLOPE(i-1) + ML*ZX(i); ); ENDSUBMODEL
CALC: @SET( 'TERSEO',2); ! Output level (0:verb, 1:terse, 2:only errors, 3:none); @SET( 'IPTOLR', .9); ! Set ending relative optimality tolerance; @SET( 'TIM2RL', TIMEUB); ! Time in seconds to apply optimality tolerance; ! ML = estimate of max difference of two slopes; ML = @MAX( SAMPLE(i) | i #GT# 1: 2*@ABS((Y(i)-Y(i-1))/(X(i)-X(i-1)))); ! @GEN( FITIT); @SOLVE( FITIT); ! Generate a graph/chart of the original and fitted curve; @CHARTCURVE( 'Fitting a Piece-wise Linear Curve,'+@FORMAT(NSEG,"2.0f")+' segments. Max error= '+@FORMAT(maxerr,"10.5f"), 'x-axis', 'y-axis', 'Original', x, y, 'Fitted', x, yfit); @WRITE(' With ', NSEG,' segments: Errors of approximation: MAX= ', @FORMAT( MAXERR,'9.5f'),', Avg= ',@FORMAT( AVGERR,'9.5f'), @NEWLINE(1)); @WRITE(' The break points of the piecewise linear fitted curve are:', @NEWLINE(1), ' X YFIT ', @NEWLINE(1)); @WRITE( ' ', @FORMAT(x(1),'12.6f'), ' ', @FORMAT(yfit(1),'12.6f'), @NEWLINE(1)); @FOR( sample(i) | ZX(i) #GT# 0.5: @WRITE( ' ', @FORMAT(x(i),'12.6f'), ' ', @FORMAT(yfit(i),'12.6f'), @NEWLINE(1)); ); @WRITE( ' ', @FORMAT(x(NP),'12.6f'), ' ', @FORMAT(yfit(NP),'12.6f'), @NEWLINE(1)); ENDCALC