$title Piecewise Linear Approximations (IMSL,SEQ=59)

$onText
This sample problem finds the best piecewise linear approximation in terms
of the sum of absolute deviations from the sampled observation. The
calculation of the interpolation weights relies on equal intervals of the
approximation function. The sine function is also implemented as a power
series expansion to demonstrate certain language features. The problem is
solved in its primal and dual version.

IMSL, LP/PROTRAN - A Problem Solving System for Linear Programming .
Tech. rep., IMSL INC, 1984.

Keywords: linear programming, piecewise linear approximation, interpolation, mathematics
$offText

Set
   n 'x-coordinate labels for data'          / d-00*d-60 /
   m 'x-coordinate labels for approximation' / a-00*a-10 /;

Parameter
   y(n)   'data values'
   t(n)   'x-coordinate values for data'
   s(m)   'x-coordinate values for approximation'
   w(m,n) 'interpolation weights';

Scalar
   k     'n / m'
   deltn 'data increments'
   deltm 'approximation increment';

t(n)  = (ord(n) - 1)/(card(n) - 1);
s(m)  = (ord(m) - 1)/(card(m) - 1);
k = (card(n) - 1) / (card(m) - 1);
abort$[0 <> mod(k,1)] 'k should be integer', k;
deltm = 1/(card(m) - 1);
deltn = 1/(card(n) - 1);

w(m+floor((ord(n)-1)/k),n)$(ord(m) = 1) = 1 - mod(ord(n)-1,k)/k;
w(m+1,n)$w(m,n) = 1 - w(m,n);

* The function sin(x) is evaluated for x between 0 and pi.
* The sets l and r are only used for the power series approximation of sin(x).
Set
   l 'length of power series'               / 0*6  /
   r 'set needed for factorial calculation' / 0*14 /;

abort$(card(l) > 2*card(r) - 1) "incorrect approximation sets", l, r;

y(n) = sum(l, power(-1,ord(l) - 1)*power(t(n)*pi,2*ord(l) - 1)/prod(r$(ord(r) <= 2*ord(l) - 1), ord(r)));
y(n) = round(y(n),6);

display y, t, s, deltm, deltn, w;

Parameter test(n,*) 'comparison of approximating sin()';
test(n,"power-ser") = y(n);
test(n,"sinus-fun") = sin(t(n)*pi);
test(n,"error")     = test(n,"sinus-fun") - test(n,"power-ser");

display test;

Variable
   ym(m) 'approximation values'
   dp(n) 'positive deviation'
   dn(n) 'negative deviation'
   tdev  'total deviation'
   tdual 'total dual value'
   z(n)  'dual values of deviations';

Positive Variable dp, dn;

Equation
   ddev(n)  'deviation definitions'
   ddual(m) 'dual definition'
   dtdev    'total dev definition'
   dtdual   'total dual definition';

ddev(n)..  sum(m, w(m,n)*ym(m)) - y(n) =e= dp(n) - dn(n);

ddual(m).. sum(n, w(m,n)*z(n)) =e= 0;

dtdev..    tdev  =e= sum(n, dp(n) + dn(n));

dtdual..   tdual =e= sum(n, y(n)*z(n)) + 1;

Model
   primal / ddev, dtdev   /
   dual   / ddual, dtdual /;

solve primal using lp minimizing tdev;

Parameter
   prep(n,*) 'primal solution report'
   primaldev 'sum of absolute deviations from primal solution';

prep(n,"t")   = t(n);
prep(n,"y")   = y(n);
prep(n,"dev") = dp.l(n) - dn.l(n);
primaldev     = sum(n, abs(prep(n,"dev")));

display prep, primaldev;

z.lo(n) = -1;
z.up(n) =  1;

solve dual using lp maximizing tdual;

Parameter
   drep(n,*) 'dual solution report'
   dualdev   'sum of absolute deviations from dual solution';

drep(n,"t")  = t(n);
drep(n,"y")  = y(n);
drep(n,"dev")= sum(m, -ddual.m(m)*w(m,n)) - y(n);
dualdev      = sum(n, abs(drep(n,"dev")));

display drep, dualdev;