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.
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
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
/;