This smells like a problem that might be NP-hard, so you'll probably have to
use one of the standard methods for dealing with intractability: e.g., devise
a heuristic, or formulate it so you can apply an existing solver.
I suggest that you try formulating it as a mixed integer linear programming
(ILP) problem, and if that's too slow, try simulated annealing. I'll sketch
both approaches.
Integer programming
We'll think of our task as selecting a subset of the edges, so that they can
be pieced together into paths as required. Introduce a zero-or-one variable
for each edge to indicate whether or not it is selected. In other words, for
each pair of vertices $i,j$, introduce a zero-or-one variable $x_{i,j}$, with
the intended meaning that $x_{i,j}=1$ implies the edge $(i,j$) is included.
We'll add a bunch of inequalities that capture the requirement that the
selected edges form a bunch of paths of the desired form:
Each vertex (other than $r$) has exactly one edge out of it:
$\sum_j x_{i,j} = 1$ for each $i \ne r$.
Each vertex (other than $r$) has at most one edge into it:
$\sum_i x_{i,j} \le 1$ for each $j \ne r$.
$r$ has no edges coming out of it: $x_{r,j} = 0$ for all $j$.
For each vertex $i$ introduce a variable $y_i$ (not restricted to integer
values), with the intended meaning that $y_i$ captures the distance from
$r$ to $i$ along the harness that connects the two. We'll require that
no harness is longer than the limit: $y_i \le b$ for all $i$. Also
$y_r = 0$.
Finally, we'll ensure the consistency of the $x$'s and $y$'s by adding
$$y_i \ge y_j + \ell(i,j) - (b + \ell(i,j)) (1-x_{i,j})$$
for all $i,j$, where $\ell(i,j)$ is the length of the edge $(i,j)$.
(Notice how this works: if edge $(i,j)$ is not selected, i.e., $x_{i,j}=0$,
this inequality implies only $y_i \ge 0$; but if edge $(i,j)$ is
selected, this inequality ensures that $i$'s distance from $r$ is
$j$'s distance from $r$ plus the length of the edge between $i$ and $j$.)
Now ask an ILP solver to minimize the objective function
$\sum_{i,j} \ell(i,j) x_{i,j} $, subject to these inequalities.
The solution will correspond the optimal solution to the original problem.
You'll have to test to see whether existing off-the-shelf ILP solvers
(e.g., CPLEX) are able to solve this ILP instance efficiently. You end
up with an ILP instance with about $|V|^2$ variables and inequalities. If it takes too long, some solvers will have an option to specify a "timeout", so they give you the best solution they've found so far within the specified time, even if it isn't necessarily the globally optimal solution.
Some minor optimizations are possible. For instance, we can omit any
edge $(i,j)$ whose length exceeds $b$ (eliminating the variable $x_{i,j}$
from the linear program, as it is forced to be zero). Also, for any pair
of edges $(i,j),(k,l)$ that "cross" (i.e., their line segments intersect),
we know they can't both be simultaneously selected, i.e., we can add the
inequality $x_{i,j}+x_{k,l}\le 1$. You probably don't want to add all
such inequalities, but you could try adding a small number (e.g., for
all pairs of short edges). You can also add the inequality $y_i \ge \ell(r,i)$ for all $i$, in case that helps the solver.
This approach is also fairly flexible: for instance, we can limit the number of
paths that are allowed to be used (by adding an inequality to upper-bound
$\sum_i x_{i,r}$) or make other minor adjustments. You can limit the number of edges allowed in each path by adding variables $z_i$ that count the number of hops between $i$ and $r$ along the harness connecting the two, and adding the inequality $z_i \ge z_j + 1 - K (1-x_{i,j})$ where $K$ is a suitably large constant (larger than the maximum number of edges per path) and then adding inequalities to upper-bound all of the $z_i$'s.
Simulated annealing
Another option is to use some form of local search, such as simulated
annealing. The idea is that we start from a candidate solution and at
each step try making a small, randomly selected, local change and check
whether it is an improvement or not.
Part of the art of simulated annealing is to choose a set of allowable
changes; at each step, you'll randomly select from this set to obtain
your candidate change. In this case, one reasonable class of "changes"
is: Pick a pair of vertices (probably ones that aren't too far apart),
and swap which paths they are in. For example, suppose I've selected
vertex $b$ and vertex $e$, where the path containing $b$ looks like
$\cdots \to a \to b \to c \to \cdots$ and the path containing $e$ looks
like $\cdots \to d \to e \to f \to \cdots$. Then after the swap operation
you'll end up with the two paths $\cdots \to a \to e \to c \to \cdots$
and $\cdots \to d \to b \to f \to \cdots$. You can then evaluate whether
this increased or decreased the objective function. You might want to
experiment with the distribution on pairs of vertices (do you try to
make nearby vertices more likely to be selected than far-apart ones?).
There may be other operations that would be useful to add and make
simulated annealing work better. Selecting an appropriate set of
operations and choosing the parameters for simulated annealing is largely an art, and may require some experimentation. It's hard to predict how well this will work without trying it.
