Finding bounds of a set defined by two affine functions and bounded variables

Question

Suppose I have two affine functions $f:\mathbb{R}^n\to \mathbb{R}$

$f(\vec{x}) = a_1\cdot x_1 + \dots + a_n \cdot x_n +q$
and
$h(\vec{x}) = b_1 \cdot x_1 + \dots + b_n \cdot x_n + p$

where $x_i$ are bounded varibles (for each $x_i$ there is a Interval with $x_i\in[l_i, u_i]$) . What im trying to find are the bounds of the set $\{ f(x) \,| \,h(x)=0\}$.

First i was trying to isolate one variable (for instance $x_1$) in $h(\vec{x})=0$
$\Leftrightarrow 0 = b_1 \cdot x_1 + \dots + b_n \cdot x_n + p $
and then replace it in $f(\vec{x})$.
I then used interval arithmetic to replace the variables by their bounds to get the bounds of the set I defined above. But the result I get depends on which variable i'm isolating as they are not part of the resulting equation. So this can't be the solution.

Then I thought about interpreting it geometrically and calculate the intersection of two hyperplanes, but then I can't seem to figure out how to calculate the bounds.

I would really appreciate if someone could give me a hint here on how to correctly calculate the bound of the set $\{ f(x) \,| \,h(x)=0\}$.

Thank you for your time and help!

Answer 1

Below is an algorithmic approach. It's exponental in the dimension $n$, so that it's only feasible for small dimensions. What it does is

• Iterate over the $2^n$ corners of an $n$-dimensinal cube, and for each corner iterates over adjacent edges. There's a total of $n\cdot2^{n-1}$ edges.

• The hyper-cube (resp. its corners) are mapped to the domain of $h$ which is the Cartesian product over all $[\ell_i, u_i]$. Then $h$ is evaluated at the edges (resp. their end-points) and checked whether $h(x) = 0$ on the edge by means of a sign-flip, i.e. one has found $x_1$, $x_2$ such that $h(x_1)h(x_2)\leqslant 0$.

• If $h$ flips sign from $x_1$ to $x_2$, then calculate $t\in[0,1]$ such that $h(x=tx_1 + (1-t)x_2) = 0$. This is straight forward as $h$ is linear and edges are lines: $t = h(x_2) \,/\, (h(x_2)-h(x_1))$.

• Compute $f(x) = f(tx_1 + (1-t)x_2) = tf(x_1) + (1-t)f(x_2)$ and accumulate it into $f_\min$ and $f_\max$.

The code uses implicitly that the domain $D$ restricted to $h(x)=0$ is a convex polytope, that $h$ is linear and that $f$ is linear. In particular, I made the assumption that if $H=\{x\mid h(x)=0\}$ is non-empty, then $f|_H$ attains its extrema on the corners of $H$.

As input I used the values from your other comment; the output is f_min, f_max = 3.0 7.5.

#!/usr/bin/env python3

# f(x) = a[] * x[] + q
# h(x) = b[] * x[] + p

a, q = [2,-3], 2
b, p = [-2,6], 1

l = [0,0]
u = [5,2]

dim = len(a)
assert len(a)==len(b) and len(b)==len(l) and len(l)==len(u)

def prod (x,y):
    """ Return scalar product of 2 vectors X and Y."""
    assert len(x) == len(y)
    s = 0
    for i in range(len(x)):
        s += x[i] * y[i]
    return s

def f(x): return q + prod(a,x)

def h(x): return p + prod(b,x)

def get_xs (lo, hi, c):
    """C in [0, 2**dim) represents a corner of a dim-dimensiomal cube.
       Map corners like 001 to lists [lo[0], lo[1], hi[2]] etc., i.e. if
       the i-th bit of C is 0, then the i-th list element is lo[i];
       otherwise the i-th list element is hi[i]."""
    xs = []
    for i in range (dim):
        xs.append (lo[i] if c % 2 == 0 else hi[i])
        c //= 2
    return xs

# Get corners in R^n that represent the domain for h().
xs = tuple (get_xs (l, u, bits) for bits in range (2**dim))

# Pre-compute h(x) at all corners x.
h_x = tuple (h(x) for x in xs)

# Pre-compute f(x) for all corners x.
f_x = tuple (f(x) for x in xs)

# Now investigate all vertices of the n-cube.

f_min, f_max = None, None

for x1 in range (2**dim):
    for b in range (dim):
        # Corner x2 is going to be adjacent to corner x1: Flip one bit.
        # Only use edges with x2 > x1 so that no edge is used twice.
        x2 = x1 ^ (1 << b)
        if x2 <= x1:
            continue

        # Evaluate h() at the ends of the edge.
        h1, h2 = h_x[x1], h_x[x2]
        
        # If h() does not change sign, then h() != 0 on the edge.
        # Skip such edges.
        if (h1 > 0 and h2 > 0) or (h1 < 0 and h2 < 0):
            continue

        # Work out t in [0,1] such that h(x) = 0,  x = t*x1 + (1-t)*x2.
        v = h2 - h1    # = b * (x2 - x1)
        if v == 0 and h2 != 0:
            continue
        t = h2 / v if v != 0 else 0

        # x = t*x1 + (1-t)*x2
        # Compute f(x) = f(t*x1 + (1-t)*x2) = t*f(x1) + (1-t)*f(x2).
        fx = t * f_x[x1] + (1-t) * f_x[x2]
        f_min = fx if f_min is None else min (f_min, fx)
        f_max = fx if f_max is None else max (f_max, fx)

print ("f_min, f_max =", f_min, f_max)