As you say, four points (without their associated "$t$" values) are not enough to define a planar cubic Bezier curve.
In fact, 6 points are needed. But, on the other hand, the 6 points have to satisfy certain conditions, or else the required curve does not exist. In your case, you know that the curve exists, because it was used to generate the points in the first place.
The full story is given in this paper: J. Kozak, M. Krajnc, Geometric interpolation by planar cubic polynomial curves, Comput. Aided Geom. Des., 24 (2007), pp. 67-78. You can get a copy of the paper at the authors website.
Another approach (instead of the clever techniques suggested in the paper) is to use brute-force numerical methods. Start with 6 points $X_1, \ldots, X_6$ and 6 parameter values $t_1, \ldots, t_6$. We can assume that $t_1 = 0$ and $t_6 = 1$, but $t_2, t_3, t_4, t_5$ are unknown.We want to find control points $P_1, P_2, P_3, P_4$ for our curve. Obviously we can set $P_1 = X_1$ and $P_4 = X_6$, but $P_2$ and $P_3$ are unknown. We can (numerically) find values of $t_2, t_3, t_4, t_5$ and $P_2$ and $P_3$ that minimize the quantity
$$ \sum_{i=2}^{i=5}\Vert C(t_i) - X_i \Vert^2 $$
where $C(t)$ is the Bezier curve with control points $P_1, P_2, P_3, P_4$. This is a fairly nasty non-linear problem, but a good numerical algorithm will usually converge to the desired solution, given a decent starting point.
Another numerical approach: with the same notation we used above, solve the equations:
$$ C(t_i) = X_i \quad (i = 2,3,4,5)$$
Since the $X_i$ are 2D points, there are 8 equations here. We have 8 unknowns -- $t_2, t_3, t_4, t_5$ and the $x$ and $y$ coordinates of $P_2$ and $P_3$. The equations are non-linear, but a good numerical root-finding package should be able to handle them.
A much better approach is to get the Bezier curves from the CAD package. If the design guys won't give you the curves, you should charge them more money for your services, since they are making your life so much more difficult. Their choice -- curves or cash :-)
Solved it myself, as below. I don't feel that this solution is "right" somehow, because it treats the interior points $\mathbf B$ and $\mathbf C$ very asymmetrically, so improvements are invited. But, anyway ...
Let's use the notation from the last few paragraphs of the question.
It's clear that $\mathbf P(0) = \mathbf A$, and $\mathbf P(1) = \mathbf D$, so these two points are interpolated, already, and we only have to worry about the other two points, $\mathbf C$ and $\mathbf D$.
First, we find numbers $h$ and $k$ such that
$$\mathbf C = \mathbf A + h(\mathbf B - \mathbf A) + k(\mathbf D - \mathbf A)
= (1-h-k)\mathbf A + h\mathbf B + k\mathbf D$$
This is possible provided that $\mathbf A$, $\mathbf B$, $\mathbf D$ are not collinear. Then, since $\mathbf P(v) = \mathbf C$, we have
$$(1-h-k)\mathbf A + h\mathbf B + k\mathbf D =
(1-v)^2 \mathbf A + 2v(1-v) \mathbf P + v^2 \mathbf D $$
But, since $\mathbf P(u) = \mathbf B$, we know that
$$ \mathbf B = (1-u)^2 \mathbf A + 2u(1-u) \mathbf P + u^2 \mathbf D$$
Substituting for $\mathbf B$ on the left-hand side, and equating coefficients of $\mathbf A$, $\mathbf P$, $\mathbf D$ gives
$$(1-v)^2 = 1 - h - k +h(1-u)^2 $$
$$2v(1-v) = 2hu(1-u)$$
$$v^2 = hu^2 + k $$
We can easily eliminate $v$ from these last three equations using the fact that $[2v(1-v)]^2 = 4[v^2][(1-v)^2]$. We get:
$$[2hu(1-u)]^2 = 4[hu^2 + k][1 - h - k +h(1-u)^2]$$
After a little algebra, this reduces to:
$$h(1-h)u^2 - 2hku + k(1-k) = 0$$
So, we solve this quadratic for $u$, and then get the unknown interior control point $\mathbf P$ from
$$ \mathbf P = \frac{\mathbf B - (1-u)^2 \mathbf A -u^2 \mathbf D}{2u(1-u)}$$
The number of solutions depends on the number of real solutions of the quadratic. Quite often, you can draw two different quadratic Bezier curves through the four given points, as explained in the paper cited in the question.
Best Answer
The locus
After doing some computations in the spirit (and building on the code) of my vertex anser, I found the locus of $P$ to be a quadric
$$ a\,x^2 + b\,y^2 + c\,xy + d\,x + e\,y + f = 0 $$
with the following parameters:
\begin{align*} a&=4\left(-A_y\,B_y + B_y^2 + A_y\,C_y - B_y\,C_y\right)\\ b&=4\left(-A_x\,B_x + B_x^2 + A_x\,C_x - B_x\,C_x\right)\\ c&=4\left(A_y\,B_x + A_x\,B_y - 2\,B_x\,B_y - A_y\,C_x + B_y\,C_x - A_x\,C_y + B_x\,C_y\right)\\ d&=4\left(A_y\,B_x\,B_y - A_x\,B_y^2 + A_y\,B_y\,C_x - B_y^2\,C_x - 2\,A_y\,B_x\,C_y + A_x\,B_y\,C_y + B_x\,B_y\,C_y\right)\\ e&=4\left(-A_y\,B_x^2 + A_x\,B_x\,B_y + A_y\,B_x\,C_x - 2\,A_x\,B_y\,C_x + B_x\,B_y\,C_x + A_x\,B_x\,C_y - B_x^2\,C_y\right)\\ f&=A_y^2\,B_x^2 - 2\,A_x\,A_y\,B_x\,B_y + A_x^2\,B_y^2 - 2\,A_y^2\,B_x\,C_x + 2\,A_x\,A_y\,B_y\,C_x \\&\quad{} - 2\,A_y\,B_x\,B_y\,C_x + 2\,A_x\,B_y^2\,C_x + A_y^2\,C_x^2 - 2\,A_y\,B_y\,C_x^2 + B_y^2\,C_x^2 \\&\quad{} + 2\,A_x\,A_y\,B_x\,C_y + 2\,A_y\,B_x^2\,C_y - 2\,A_x^2\,B_y\,C_y - 2\,A_x\,B_x\,B_y\,C_y \\&\quad{} - 2\,A_x\,A_y\,C_x\,C_y + 2\,A_y\,B_x\,C_x\,C_y + 2\,A_x\,B_y\,C_x\,C_y - 2\,B_x\,B_y\,C_x\,C_y \\&\quad{} + A_x^2\,C_y^2 - 2\,A_x\,B_x\,C_y^2 + B_x^2\,C_y^2 \end{align*}
Assumption confirmed
Using this formulation, I could verify your assumption: the conic is indeed a hyperbola, with the lines $AB$ and $CB$ as asymptotes. The way to check this is by ensuring that the lines are tangents, and that they touch the conic at infinity.
Special location
In a comment below, you mention that you might assume $B_x=B_y=0$. With that your conic will become
\begin{align*} a&=4\,A_y\,C_y\\ b&=4\,A_x\,C_x\\ c&=-4\left(A_y\,C_x + A_x\,C_y\right)\\ d&=0\\ e&=0\\ f&=A_y^2\,C_x^2 - 2\,A_x\,A_y\,C_x\,C_y + A_x^2\,C_y^2 \end{align*}
The formula is indeed a lot easier, so it might be a good starting point for a geometric description of the curve. For example, it is readily apparent that this hyperbola will be symmetric around the origin. But we already knew $B$ to be the center due to the asymptotes. You can even take this one step further and, via an affine transformation, consider only the case where
\begin{align*} A &= \begin{pmatrix}1\\0\end{pmatrix} & B &= \begin{pmatrix}0\\0\end{pmatrix} & C &= \begin{pmatrix}0\\1\end{pmatrix} \end{align*}
You end up with the hyperbola
$$4xy=1$$
From this you can e.g. conclude that the midpoint between $A$ and $C$ will also lie on that hyperbola. It will of course correspond to a Bézier curve which passes through $B$ for some $t<0$ or $t>1$. You can take that midpoint and reflect it in $B$ to obtain a point on the segment of the hyperbola you get for $0\le t\le 1$. This construction is invariant under affine transformations, so it still holds for the general case. Two asymptotes plus one point on the hyperbola amounts to five real degrees of freedom, so this should be enough to uniquely define your hyperbola.
Summary
To sum it up: the locus of the apex $P$ is the unique hyperbola with asymptotes $AB$ and $CB$ which passes through the midpoint between $A$ and $C$. The portion of it which corresponds to parabolas where $B$ lies between $A$ and $C$, i.e. is obtained in the Bézier curve for $0<t<1$, is the component of the hyperbola which does not contain that midpoint. It does contain the point obtained by reflecting that midpoint in $B$.
The code
Here is the sage code I used to obtain this representation:
Originally I had more complicated code which did not rely on the interpretation as a Bézier curve. The result was the same, though.
Generalizing to rational case
A parabola has four real degrees of freedom. If you choose a non-rational Bézier curve, you have two real degrees of freedom, but with these you not only specify a parabola but also select a start point and an end point on that parabola, so the degrees of freedom match. A conic in general has five real degrees of freedom. So if you want it to pass through three given points, that still leaves a two-parameter family of corresponding conics. Therefore your locus will not be a single curve, but either the whole plane or some portion of it.
You can define a conic using five points through which it should pass. In addition to your $A,B,C$ you might use two more points, which you move close to the end points $A$ and $C$. By making them arbitrary close (i.e. computing some limit), you can use these control points to exactly and arbitrarily determine the direction of the tangents in $A$ and $C$. Therefore you can choose any point $P$ in the plance, and find a conic through $A,B,C$ which will have $P$ as its apex. In this sense, the whole plane will be your locus.
I'm not completely sure whether a rational Bézier curve can be defined in such a way that it passes through infinity, but I believe that to be the case. If not, then there might be cases where the resulting conic would be a hyperbola, and $A,B,C$ are not all three in the same of its components. This might result in a restriction to part of the plane.