$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

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);

$onEcho > quadminos.o10
feasibility tolerance 1e-14
$offEcho

$onEcho > quadminos.o11
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
  /;