qfilter.gms : Audio filter design using quad-precision MINOS

Description

```Given:
A set i of audio frequencies,
parameters omega(i) and gamma(i), and
positive variables u1, u2, v1, v2, and intermediate variables
U(i) = 1 + u1*sqr(omega(i)) + u2*power(omega(i),4)
V(i) = 1 + v1*sqr(omega(i)) + v2*power(omega(i),4)

NLP1: Find
min  beta
s.t. abs[V(i)/U(i) - gamma(i)] <= beta   for all i

Multiply through by U(i) and replace the abs[] < beta*U with:
-beta*U <= V - gamma*U <= beta*U

If we make beta a constant we have an LP. Binary search on beta will
give us a solution without requiring a nonlinear solver.  The scaling
of the powers of omega is such that double precision is not enough to
accurately solve this model, and the typical rules of thumb about
convergence tolerances and small values (i.e. what is small enough to
treat as zero) do not apply.  Using quadruple precision we can get
some meaningful results, but we must adjust convergence tolerances
in MINOS and zero tolerances in GAMS to do so.
```

Reference

• GAMS Development Corporation, Formulation and Language Example.

Small Model of Type : GAMS

Category : GAMS Model library

Main file : qfilter.gms

``````\$title Audio Filter Design using quad-precision MINOS (QFILTER,SEQ=405)

\$onText
Given:
A set i of audio frequencies,
parameters omega(i) and gamma(i), and
positive variables u1, u2, v1, v2, and intermediate variables
U(i) = 1 + u1*sqr(omega(i)) + u2*power(omega(i),4)
V(i) = 1 + v1*sqr(omega(i)) + v2*power(omega(i),4)

NLP1: Find
min  beta
s.t. abs[V(i)/U(i) - gamma(i)] <= beta   for all i

Multiply through by U(i) and replace the abs[] < beta*U with:
-beta*U <= V - gamma*U <= beta*U

If we make beta a constant we have an LP. Binary search on beta will
give us a solution without requiring a nonlinear solver.  The scaling
of the powers of omega is such that double precision is not enough to
accurately solve this model, and the typical rules of thumb about
convergence tolerances and small values (i.e. what is small enough to
treat as zero) do not apply.  Using quadruple precision we can get
some meaningful results, but we must adjust convergence tolerances
in MINOS and zero tolerances in GAMS to do so.

Michael A. Saunders, private communication, filter.pdf: 22 Aug 2014.

Keywords: nonlinear programming, GAMS language features, audio filter design
\$offText

\$onDollar
* not all platforms support quad-precision, check that here

\$if not %HAVEQUAD%==true \$log Aborting: system.buildcode = %system.buildcode% not known to support quadminos

Set i 'audio frequencies' / 20,   125,  250,  500,   750, 1000,  1500
2000, 3000, 4000, 6000, 8000, 16000, 32000 /;

Parameter
omega(i) 'frequency value'
gamma(i) 'g-squared'
g(i)     'audiogram measurement'
/  20 1,  750 1.51991, 3000 12.3285, 16000 100
125 1, 1000 1.51991, 4000 18.7382, 32000 100
250 1, 1500 2.31013, 6000 50
500 1, 2000 3.51119, 8000 71 /;

omega(i) = i.val*2*pi;
gamma(i) = sqr(g(i));

Free     Variable beta, z;
Positive Variable u1, u2, v1, v2, U(i), V(i), x(i) 'beta * U(i)';

Equation
zDef
xDef(i)  'define x := beta*U(i) - the only nonlinear constraint'
absUp(i) 'upper half of absolute value constraint'
absLo(i) 'lower half of absolute value constraint'
UDef(i)
VDef(i);

zDef..     z    =e= beta;

xDef(i)..  x(i) =e= beta*U(i);

absUp(i).. V(i) - gamma(i)*U(i) =l=  x(i);

absLo(i).. V(i) - gamma(i)*U(i) =g= -x(i);

UDef(i)..  U(i) =e= 1 + u1*sqr(omega(i)) + u2*power(omega(i),4);

VDef(i)..  V(i) =e= 1 + v1*sqr(omega(i)) + v2*power(omega(i),4);

Model m / all /;

* compute initial values as in the reference
Scalar
a1 / 2.6e-5  /
b1 / 2.6e-4  /
a2 / 3.5e-10 /
b2 / 3.5e-8  /;

u1.l = sqr(a1) - 2*a2;
v1.l = sqr(b1) - 2*b2;
u2.l = sqr(a2);
v2.l = sqr(b2);

U.l(i) = 1 + u1.l*sqr(omega(i)) + u2.l*power(omega(i),4);
V.l(i) = 1 + v1.l*sqr(omega(i)) + v2.l*power(omega(i),4);

feasibility tolerance 1e-14
\$offEcho

feasibility tolerance 1e-17
\$offEcho

Scalar
dTol   'convergence tol' / 1e-5  /
delta  'error bound'
n      'iteration count' /   0   /
loBeta                   / 274.5 /
upBeta                   / 275   /;

option lp = quadminos, sysOut = off, solPrint = silent;
m.optFile   = 10;
m.tolProj   = 1e-24;
m.holdFixed = 1;
m.tryLinear = 1;

File fp / '' /;
fp.nr =  2;
fp.nw = 18;
fp.nd = 10;
fp.nz = 1e-30;
put fp;

* try left endpoint: should be infeasible
beta.fx = loBeta;
x.l(i)  = beta.l*U.l(i);

solve m using nlp min z;
abort\$[m.modelStat <> %modelStat.infeasible%] 'Expected left endpoint to be infeasible', loBeta;

* try right endpoint: should be feasible
beta.fx = upBeta;
x.l(i)  = beta.l*U.l(i);

solve m using nlp min z;
abort\$[m.modelStat <> %modelStat.optimal%] 'Expected right endpoint to be feasible', upBeta;

* with warm start we can solve with smaller tolerance
m.optFile = 11;
x.l(i)    = beta.l*U.l(i);

solve m using nlp min z;
abort\$[m.modelStat <> %modelStat.optimal%] 'Expected right endpoint to be feasible with small tol', upBeta;

while{1,
delta = upBeta - loBeta;
break\$(delta <= dTol);

n = n + 1;
beta.fx = loBeta + 0.5*delta;
x.l(i)  = beta.l*U.l(i);
putClose ' '
/ 'n = ', n:0:0, '  delta = ', delta, '   beta = ', beta.l /;

solve m using nlp min z;
abort\$[(m.modelStat <> %modelStat.optimal%) and
(m.modelStat <> %modelStat.infeasible%)] 'Unexpected model status', m.modelStat, beta.l;
if{(m.modelStat = %modelStat.optimal%),
upBeta = beta.l;
else
loBeta  = beta.l;
beta.fx = upBeta;
};
};

if{(m.modelStat <> %modelStat.optimal%),
x.l(i) = beta.l*U.l(i);
putClose ' '
/ 'Final re-solve at right end-point' /;

solve m using nlp min z;
abort\$[(m.modelStat <> %modelStat.optimal%)] 'Expected feasible beta on exit', m.modelStat, beta.l;
};

put ' '
/ 'Successful termination'
/ 'iteration count: ', n:10:0
/ '           dTol: ', dTol
/ '          delta: ', delta
/ '           beta: ', beta.l
/ '         loBeta: ', loBeta
/ '             u1: ', u1.l
/ '             u2: ', u2.l
/ '             v1: ', v1.l
/ '             v2: ', v2.l
/;
``````