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)
