GAMS [ Home | Downloads | Documentation | Solvers | APIs | Tools | Model Libraries | Resources | Sales | Support | Contact Us | Search ]

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.

$offtext

$ondollar
* not all platforms support quad-precision, check that here
$set HAVEQUAD false
$if DEG==%system.buildcode%  $set HAVEQUAD true
$if LEG==%system.buildcode%  $set HAVEQUAD true
$if VS8==%system.buildcode%  $set HAVEQUAD true
$if WEI==%system.buildcode%  $set HAVEQUAD true

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

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

parameters
  omega(i) 'frequency value'
  g(i)     'audiogram measurement' /
     20  1
    125  1
    250  1
    500  1
    750  1.51991
   1000  1.51991
   1500  2.31013
   2000  3.51119
   3000  12.3285
   4000  18.7382
   6000  50
   8000  71
  16000  100
  32000  100
  /
  gamma(i) 'g-squared'
  ;
omega(i) = i.val * 2*pi;
gamma(i) = sqr(g(i));

free variable
  beta
  z
  ;
positive variables
  u1
  u2
  v1
  v2
  U(i)
  V(i)
  x(i)    'beta * U(i)'
  ;
equations
  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
scalars
  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

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

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 eq %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
  /;