I've got a few comments on this. Apologies if any of this sounds rough, it's not intentional but there is a lot of material so I'm going to speak plainly.
No problem at all - I prefer frankness. Saves everyone time
Very long post (as are mine...) I'm not going to comment on it point by point (dangit I nearly ended up doing just that...)
Just a few take aways:
20 permutations vs 24 permutations - D'oh on my part.
I of course meant 24, no idea why I kept writing 20 - It was late

Nice catch.
---------
Ignore the below. Can't seem to find a
strike tag - So using color instead

Had gotten my p1-p3 vs p1-p4 permutations mixed up.
That along with the --aside-- experiments below got the thing mixed up for me - apologies.
"the order of the 3 references used to generate the candidate points is not important."
I'm sorry - But it is.
Due to the rounding of the distances - Ie. the "delta" between the "true" value and the rounded value.
Some of those deltas are larger than the others.
The formula (wikipedia) first calculates x, then y and z.
The calculation for y depends on x, and z on both x and y.
If the calculated value for x is "very accurate" (due to the deltas being small) then y will be "less off" than if two other systems (with larger deltas) had been used.
etc.
So basically - If the rounded distance to say reference star 1 (r1 in wikipedia formulas) just so happens to have a very large delta, then the trilateration is more likely to be in-accurate.
Trilateration is very sensitive to the measured distances - hence why 4 digits precision works so much better (obviously).
So yes - The order does matter (not always -again, depending on how large/different the deltas are)
-- aside
I've done tests on this as well, putting in the deltas in the formulas - To see if one could somehow "isolate" them, or at least try and use them for accuracy prediction - based on how the coordinates are related to each other (another dead-end btw

)
The formulas get extremely unwieldy ofc - But it's sorta interesting (if you are math inclined).
But put in max delta (+/-) and it's very visible what happens to the end coordinates (you get 8 formulas - for all possible combinations of max delta (+/-))
x is fairly predictive - always having a max and min range (so you can say if I err on the side of max delta on this, and min delta on these two - then i will always get the most accurate result)
y pretty much the same, but more spread out.
z on the other hand... appears entirely chaotic as you change the coordinates of the reference systems.
As said - a dead end for any predictiveness with regards to accuracy. But was a "fun" exercise in error-calculation.
-- end aside
----------
"Trilateration requires good reference systems, as has been mentioned repeatedly in this thread. I think it is better to use an algorithm that will work well with a specified set of references or that will warn when the chosen references are bad"
"I'm not sure there is much value in testing with random reference systems because we already know that can get you into trouble."
This is likely controversial - But I do not subscribe to the idea of "fixed reference systems".
The reason being my above observations regarding deltas on distances.
I know it has kind of been accepted as "the way to go" in this thread - But I'm sorry, the math says it's irrelevant.
Any set of "fixed reference systems" is just as likely to have large deltas as any randomly picked systems.
It all depends on what system you want to pin down - and if it just so happens that the coordinates for that systems will give rise to large deltas on the distances or not.
No "fixed reference systems" that one picks can ever guarantee that there wont me max deltas (large rounding errors on the distances)
And those deltas are what matters.
A "good reference system" is one that (just so happens) to give small deltas.
And that is entirely dependent on the target systems coordinates.
I'll gladly do a math post proving this if required. Let me know.
-------
Ignore the below.
Again my mixup of my p1-p3 vs p1-p4 permutations - With just p1-p3 permuted our results agree.
"I never get more than 2 coordinates (which would indicate that the order was significant)."
This flies in the face of my experimentation - With the same setup.
When you calculate r1,r2,r3 (which would be the user inputted values) do you remember to round them to 3 digits of precision before using them in the trilateration formulas?
I once made an "ooops" forgetting this - Only reason I mention it...
Otherwise I cant explain the difference in our results.
And this really is something we should look at.
Our test setup is the same (3 random systems (well 4..) - to determine 4th (that we already know what the result should be for)).
So that we get such divergent results is troubling.
------
"Our output point is an approximation of the shape generated by the intersection of those three shells. The intersection of the three shells can be quite large (if the reference stars are poorly chosen), large enough to contain multiple points on the 1/32 Ly grid."
Again, if by "non poorly chosen", you mean 3 (or more) fixed reference systems that are always used for the calculations of p0 - Then I have to disagree.
You can not pick such systems such that the delta on the distances to p0 will always be minimal.
Move p0 a bit to the side and delta-r1 will go up, and delta-r2 will go down (or also up - who knows) etc.
I do agree that any picked reference systems should be "spread around" the target system (p0) as this is more likely to reduces the area of the intersecting shells.
----
"Which means the reference stars are degenerate in some way, e.g. they all lie in a plane."
The degenerate case is if the reference stars are
*collinear*.
The reference systems p1-p3 will always be in a plane (3 points always will be)
"Trying to use linear algebra to predict if your reference systems are a good set seems like an incredibly roundabout (and difficult) way to tackle the problem"
Very much so

.
But that was to try to determine if p4 where in a plane with p1-p3, and if that was bad for the result of trilateration formulas.
The consensus in this thread seems to have been that that would be the case - as your comment just above also shows.
So I ended up disproving that - after all too much liniar algebra
Wether p1-p4 - or indeed p0-p4 all are very close to being in the same plane or not - has no impact on the final calculations.
Again - It all boils down to the deltas on the distances - Regardless the spatial positions (minus the colliniar case)
-----
"In my opinion the best measure of accuracy is how well the distances derived from the calculated coordinates match the input distances. "
The difference being the deltas as I call them.
And yes we agree - The smaller the deltas, the more accurate the result is likely to be.
"If the distances between p0 and p1 .. p5 match the supplied distances d1 .. d5 to the precision of d1 .. d5 (3 decimal places in the galaxy map) then p0 is as accurate as we need it to be (and in fact as accurate as we can test for)."
This is one test I haven't done for some odd reason.
Comparing the p1-p3 distances after the calculation...
Now I have something to do the rest of the day instead of writing mega posts here.
I'm going to correlate this with the results of my other tests - and see if they match up or not.
My suspicion is that that still leaves room for error though...
Replacing one of the reference systems (p1-p3) with another system (and then using my draconian acceptance measures), I still think will be better.
I might prove myself wrong though - I'll let you know what I come up with.
-------
resolved - see earlier strikeouts
As a last PS in this ridiculous long post, so that it doesn't disappear in the wall of text: I really am very very interested in finding out why our tests give different results, as mentioned further up (the "maybe you didn't round the distances" suggestion - I'm sure you did but then what?)