imsl.gms : Piecewise Linear Approximation
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.
Reference:
- IMSL Inc, LP/PROTRAN - A Problem Solving System for Linear Programming . Tech. rep., IMSL, 1984.
Small Model of Type: LP
<![CDATA[
$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.
$Offtext
Sets n x-coordinate labels for data / d-00 * d-60 /
m x-coordinate labels for approximation / a-00 * a-10 /
Parameters y(n) data values
t(n) x-coordinate values for data
s(m) x-coordinate values for approximation
w(m,n) interpolation weights
Scalars deltn data increments
deltm approximation increment ;
t(n) = (ord(n)-1)/(card(n)-1); s(m) = (ord(m)-1)/(card(m)-1) ;
deltm = 1/(card(m)-1); deltn = 1/(card(n)-1);
w(m+floor(t(n)/deltm),n)$(ord(m) eq 1) = 1 - mod(t(n),deltm)/deltm;
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).
Sets l length of power series / 0*6 /
r set needed for factorial calculation / 0*14 / ;
Abort$(card(l) gt 2*card(r)-1) "incorrect approximation sets",l,r;
$Offdigit
Scalar pi / 3.1415926538979 /;
$Ondigit
y(n) = sum(l, power(-1,ord(l)-1)*power(t(n)*pi,2*ord(l)-1) / prod(r$(ord(r) le 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;
Variables 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 variables dp, dn;
Equations 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;
]]></plaintext></body></div>