Finding stages in routing algorithm
I don't know whether there is an efficient (polynomial-time) algorithm for this problem.   I wouldn't be surprised if it is NP-hard.

One approach would be to formulate this as an instance of integer linear programming, and then apply an off-the-shelf ILP solver to find the best solution it can within some fixed amount of time.  This gets awfully ugly.  You could introduce a bunch of zero-or-one variables to represent the solution, write down constraints (inequalities) to capture that they are consistent, express the cost of a solution as a linear combination of these variables, and then use the ILP solver to minimize the cost.  However, writing down the set of inequalities might get super messy.

Alternatively, it might be possible to apply heuristics for weighted set cover or weighted exact cover to this problem.  Given any candidate combination of 3..5 terminals and a candidate location for the intermediate point, you can easily compute the cost of connecting them.  Consider this as a "set" that covers those terminals and that intermediate node, and whose cost is equal to the computed cost.  There are $n_1 \times ({n_2 \choose 3} + {n_2 \choose 4} + {n_2 \choose 5})$ many such "sets", where $n_1$ is the number of internal nodes and $n_2$ is the number of terminals.  This is $O(n^6)$, where $n$ is the number of nodes in the graph.  Now the goal is to find a cover, i.e., a subset of "sets" that covers each terminal exactly once and that covers each internal node at most once.

In principle you could then try applying any algorithm for exact cover.  However, I suspect this will be far too slow, as there will be a few hundred billion possible "sets" or so.

An alternative might be to apply a heuristic or approximation algorithm for set cover.  For instance, consider the greedy algorithm for set cover: at each step, it chooses a set that minimizes the cost of the set divided by the number of terminals that are newly covered by that set and haven't been covered yet (best "bang-for-the-buck").  This gives a solution that is not too much more expensive than the optimal solution.

Moreover, you can implement the greedy algorithm efficiently, without enumerating all sets.  At each step, we have already covered some of the terminals and some of the internal nodes, and we want to choose the next set with the best bang-for-the-buck ratio.  This seems feasible.  One naive way to do that is to iterate over all possible choices for the location of the intermediate point, and for the number $k$ of terminals (3..5), and then find the $k$ closest uncovered terminals to that intermediate point; each gives you a candidate set.  Also find the closest uncovered terminal to $s$; this gives you another candidate set (with no intermediate point).  Among these candidates, take the one with the lowest bang-for-the-buck ratio and add it to your solution so far.  If you precompute all-pairs shortest paths (e.g., with the Floyd-Warshall algorithm), then each step can be implemented in $O(n^2)$ time.  You need to do at most $n$ steps, so the total running time will be $O(n^3)$.

This might give a highly suboptimal solution, but it is at least one solution, and it should be relatively efficient with your parameter settings.

Alternatively, it might be possible to use local search or simulated annealing to optimize your solution.  The key requirement is to have a way, given a current solution, to find a way to perturb the solution in a random way.  The algorithm then tries this many times, hopefully gradually decreasing the cost of the solution over time.  In your situation, one way to random perturb a candidate solution is to either (1) randomly pick three terminals that currently aren't part of any pack and merge them into a pack, (2) randomly pick a terminal that currently isn't part of any pack and a pack of size 3..4, and merge them into a single pack, or (3) randomly split a pack of size 4..5 into a pack of size 3..4 and a terminal that isn't part of any pack.  You could then try simulated annealing together with such a choice of perturbation (neighborhood relation).  There might be smarter ways to choose the perturbation that makes those algorithms more effective.

Maybe someone else will be able to come up with a better approach.