$\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }$

Purpose
The function $f : \B{R}^3 \rightarrow \B{R}$ defined by $$f( x_0, x_1 ) = ( x_0^2 + x_1^2 ) / 2 + | x_0 - 5 | + | x_1 + 5 |$$ For this case, the abs_min_quad object should be equal to the function itself. In addition, the function is convex and abs_min_quad should find its global minimizer. The minimizer of this function is $x_0 = 1$, $x_1 = -1$.

Source


namespace {
{     size_t n = x.size();
size_t s = u.size();
for(size_t j = 0; j < n; j++)
xu[j] = x[j];
for(size_t j = 0; j < s; j++)
xu[n + j] = u[j];
return xu;
}
}
{     bool ok = true;
//
//
//
size_t level = 0;     // level of tracing
size_t n     = 2;     // size of x
size_t m     = 1;     // size of y
size_t s     = 2 ;    // number of data points and absolute values
//
// record the function f(x)
for(size_t j = 0; j < n; j++)

// create its abs_normal representation in g, a
f.abs_normal_fun(g, a);

// check dimension of domain and range space for g
ok &= g.Domain() == n + s;
ok &= g.Range()  == m + s;

// check dimension of domain and range space for a
ok &= a.Domain() == n;
ok &= a.Range()  == s;

// --------------------------------------------------------------------
// Choose the point x_hat = 0
d_vector x_hat(n);
for(size_t j = 0; j < n; j++)
x_hat[j] = 0.0;

// value of a_hat = a(x_hat)
d_vector a_hat = a.Forward(0, x_hat);

// (x_hat, a_hat)
d_vector xu_hat = join(x_hat, a_hat);

// value of g[ x_hat, a_hat ]
d_vector g_hat = g.Forward(0, xu_hat);

// Jacobian of g[ x_hat, a_hat ]
d_vector g_jac = g.Jacobian(xu_hat);

// trust region bound
d_vector bound(n);
for(size_t j = 0; j < n; j++)
bound[j] = 10.0;

// convergence criteria
d_vector epsilon(2);
double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
epsilon[0]   = eps99;
epsilon[1]   = eps99;

// maximum number of iterations
s_vector maxitr(2);
maxitr[0] = 10; // maximum number of abs_min_quad iterations
maxitr[1] = 35; // maximum number of qp_interior iterations

// set Hessian equal to identity matrix I
d_vector hessian(n * n);
for(size_t i = 0; i < n; i++)
{     for(size_t j = 0; j < n; j++)
hessian[i * n + j] = 0.0;
hessian[i * n + i] = 1.0;
}

// minimize the approxiamtion for f (which is equal to f for this case)
d_vector delta_x(n);
}