Q:

Location Selection Algorithm in Solar Engineering

This is a practical problem in energy generation with heliostats.
We have a number of heliostats basically forming the shape of a doughnut. The facility needs to deploy hubs on those heliostats. One hub can support 12 harnesses (a kind of cable) and one harness can support 16 heliostats; namely, one hub can support 192 heliostats. The cost of deploying a hub and harness per meter are given. In other words, every harness starts from a hub and can connecting at most 16 heliostats. Deployment of a hub results in a cost. Also, the cost of harnesses result from its distance from the hub to the furthest heliostat on this harness. However, there is no distance from heliostats to its connecting harness. 
Currently my approach is to use k-means to cluster the whole region and run greedy on distance from harness for each split region, though I do not think this is a good approach. 

This is a pattern generated by my algorithm. As you can see, some harness only support 1 or 2 heliostats and result in a waste of harness. 
Our job is to optimize the cost. In the end, the algorithm should return locations of hubs, harness topology patterns, and the total cost. 
Can anyone enlighten me potential algorithms on this?
(Every point on the image is a heliostat, there are 40981 in total)
----------- Updated 14/12/2016 -----------
More details on the question, thanks D.W's suggestion: 

Coordinates of heliostats are given. These coordinates correspond to real-world metrics (meters).
The placement of hubs are not constrained. However, one harness cannot be longer than 45 meters due to voltage requirement. 

----------- Updated 15/12/2016 -----------

The classification on "furthest" in "the cost of harnesses result from its distance from the hub to the furthest heliostat on this harness": Furthest is in the sense of "last", not "furthest" in the sense of "having the largest Euclidean distance from the hub"

A:

The problem looks messy, so there probably isn't a clean algorithm to
produce a globally optimal solution.  Instead, I suggest you use some
combination of heuristics to try to find a nearly-optimal solution
in a reasonable amount of time, perhaps by iteratively making small
incremental improvements.
Any candidate solution is comprised of three elements:

Location of the hubs.
Assignment of each heliostat to a hub.
Routing of the harnesses.

One way to make an incremental improvement is to repeatedly fix two
of these three and optimize the other (holding the other two fixed),
and rotate through which one you change.  You seem to have a reasonable
strategy to obtain an initial candidate solution, and you could iterate from
there.  I'll sketch in more detail each of these three optimizations,
which you could rotate through repeatedly until you cannot make any
further improvement:
Harness optimization. Hold the location of the hubs and the
assignment of which hub each heliostat is associated with fixed
(to be the same as in the current candidate solution).  The goal
is optimize the routing of the harness wiring.  For
candidate algorithms for that, see the other question.  My
suspicion is that this step is the one where there are the greatest
opportunities for improving your solution, and where it's worth
devoting the most effort to.
Hub location optimization. Hold the heliostat-to-hub assignments
and harness routing fixed.  Pick a single hub and the heliostats
associated with it.  Then it's easy to fine-tune the location
of the hub to minimize overall cost, under the assumption that you
don't change the topology/routing of the harnesses.
In particular, each
harness contributes one edge from the hub $r$ to the heliostats
$h_1,h_2,\dots,h_{12}$ it is immediately connected to.  Now we want
to find a location for $r$ that minimizes the sum of distances
$d(r,h_1)+\dots + d(r,h_{12})$.  This is a sum-of-squares problem
so it can be solved through standard sum-of-squares optimization,
or simply by solving it with gradient descent; since we're in only
2 dimensions I would expect any method to converge rapidly.
You can do this separately for each hub.
It would also be possible to try a random perturbation to the
location of the hub, re-apply harness optimization, see whether
this reduced the overall cost, and if so accept that change to the
hub location.  This might be more expensive than it is worth.
Heliostat assignment to hubs. Hold the location of hubs fixed,
and assume we're not going to make radical changes to the harness
topology.  We can try to fine-tune the assignment of heliostats to
hubs in a number of ways.
One approach is to pick a heliostat $h$ that is currently assigned to
hub $r_1$ but is relatively close to another heliostat that's wired
to a different hub $r_2$, try swapping which hub it is associated
with, and see if this leads to any reduction in cost.  To determine
whether it leads to a reduction in cost you could re-run
harness optimization from scratch in both cases (which might be slow).
Or, as a fast heuristic, you could remove $h$ from its current harness,
find the nearest heliostat $h'$ that's currently connected to $r_2$,
attach $h$ to $h'$ (splicing it into the harness for $h'$), and see
whether this reduces the cost.
Another approach would be to do a more ambitious optimization.
Pick a pair of hubs $r_1,r_2$ and the set of heliostats currently
connected to either of them.  Now try to apply joint harness optimization
to that set of heliostats and those two hubs simultaneously,
simultaneously optimizing the routing of the heliostats and also
which hub each heliostat is connected.  The methods in that other
question can be generalized to solve this problem.  When finding the
solution we may find that one or two heliostats switch which hub they
are connected to.  Now repeat this for every pair of adjacent hubs.
However, my guess is that heliostat assignment might not yield many
gains (compared to the naive method of assigning each heliostat to
the hub it is closest to), so it might not be worth implementing these
more sophisticated methods for heliostat assignment.
Overall algorithm. In summary, we start from some initial
solution (e.g., selected using the method described in the question)
and then repeatedly rotate through the following three operations:

Harness optimization: fine-tune the routing of the harnesses.
Hub location optimization: fine-tune the location of the hubs.
Heliostat re-assignment: fune-tune the assignment of heliostats to hubs.

There is no guarantee that this will lead to a global optimum, i.e.,
the absolute lowest cost possible.  This could get stuck in a local
minimum.  There are various methods for dealing with that (e.g., applying
random perturbations; re-starting from multiple randomly chosen initial
solutions), which you could experiment with as well.
Hopefully this will help you find a solution that is "good enough",
or at least, get you started on some methods to try.  Good luck!

