# Why does candidate optimal group computation produce strange results when thrusters are too close together?

I have two thruster configurations I'm testing. One is almost symmetric about its center-of-gravity (it's a little bottom-heavy, but the thrusters themselves are the same fore and aft). The other has four additional thrusters aft for more translational delta-v.

That image shows the asymmetric configuration. The symmetric one is pretty similar, but without the 15cm separation between the second set of translational thrusters. The radius is about 2 and a quarter meters, and the length is around 9 meters; CG is 4 meters from the rear.

But I get really strange results with the asymmetric configuration.

## Candidate Optimal Groups

Unfortunately, I haven't been able to find the original paper for COGs. But the algorithm itself is described in two places I've found:

• Shoemaker, D. 2013. A survey of spacecraft jet selection logic algorithms. In AAS Guidance & Control Conference, Breckenridge, CO.
• Wie Bong. 2008. Chapter 7: Rotational Maneuvers and Attitude Control. Space Vehicle Dynamics and Control. (This one seems to be in PDF form on the Google. In my copy, the algorithm is given in section 7.6.2, starting on p. 474.)

The plots below are my attempt to reproduce those given in the Shoemaker review. The rough idea is that you have a command vector $\mathbf{b}$ of size $d$ (one entry per degree of freedom), normalized to length 1, and a thruster on-time array $\mathbf{x}$ of length $N$ (one for each thruster). You're trying to solve

$$\mathbf{A}\mathbf{x} = \mathbf{b}$$ $$\textrm{subject to the constraint that }\mathbf{x} \geq \mathbf{0}$$

You can't really solve that directly, since it's under-constrained. So instead we assume that for $d$ degrees of freedom, we need $d$ thrusters to be turned on, and partition $\mathbf{x}$ into $\mathbf{\hat{x}}$ (for selected thrusters) and $\mathbf{\tilde{x}}$ (for thrusters we're definitely keeping off).

We can accordingly partition $\mathbf{A}$ into corresponding $\mathbf{\hat{A}}$ (square) and $\mathbf{\tilde{A}}$ for each $\mathbf{x}$ partition. Using that, we can solve

$$\mathbf{\hat{A}}\mathbf{\hat{x}} = \mathbf{b}\,\textrm{.}$$

That gives all the candidate groups of thrusters, only some of which are optimal. The full algorithm describes how you reduce the search space, eliminating groups that are clearly non-optimal — those where $\mathbf{\hat{A}}$ is singular, where turning on some thruster not in $\mathbf{\hat{x}}$ produces a better score, or where any value in $\mathbf{\hat{x}}$ is less than $-\epsilon$.

The final candidate optimal groups are those which maximize the score $\mathbf{b}\cdot(\mathbf{\hat{c}}\mathbf{\hat{A}}^{-1})^\mathrm{T}$, where $\mathbf{c}$ is the cost vector (basically the flow rate for each thruster when turned on) and appropriately partitioned as with $\mathbf{A}$ and $\mathbf{x}$. (Note that most people assume all costs are 1 for simplicity.)

So! That brings us to the plots.

I'm only showing 3 DOF here, so I've done the computations for rotation-only or translation-only. Thus, commands $\mathbf{b}$ are of size 3. So, for pure $\hat{x}$ translation, $\mathbf{b}_\textrm{translation} = \left[ 1\,0\,0 \right]$; for pure $\hat{y}$ translation, $\mathbf{b}_\textrm{translation} = \left[0\,1\,0\right]$; etc. I do the same for rotation.

So I generate a unit sphere mesh of commands. For each COG, I find the optimal set (or sets) of thrusters to use, and compute the score and the total on-times (the sum of $\mathbf{\hat{x}}$). The on-times figure, then, shows the unit sphere scaled by the sum of computed on-times. The authority figure is supposed to be, according to my reading, the unit sphere points scaled by the reciprocal score. The color map is computed identically.

## Results

The translation results are identical, as expected:

I've been struggling with this for a few days. I suspect the issue relates to the LU decomposition I'm using for inverting the candidate authority matrices. Why have I come to this conclusion? In my translational authority matrix, all of my cells have only two magnitudes (0 and a hundred something). But in my rotational authority matrix, I get a couple of magnitudes which are small (~8) and a bunch that are bigger (200 to 600).

Has anyone else seen anything like this? I've now tried two different LU decomp algorithms, and they produce the same output. (Both use partial pivoting, and they should be stable; one is from GSL. I'm also checking the reciprocal condition number of the result and discarding candidate groups where rcond < eps.)

When I first wrote the code in Python, I used SVD instead of LU decomposition, and it didn't seem to suffer from this particular problem. But now that I'm writing in a compiled language, SVD's relative algorithmic complexity is of concern, so I'd rather use LU, if possible.

• Can we see your code? – Schlusstein Mar 23 '17 at 0:09
• @uhoh I have added a whole bunch more information to answer your questions. I hope this helps and isn't overwhelming. – Doctor Mohawk Mar 23 '17 at 15:28
• @uhoh Just added the word "score" for clarification where it's computed. It's the thing we're trying to maximize. – Doctor Mohawk Mar 23 '17 at 15:38
• Oh, $\mathbf{b}\cdot(\mathbf{\hat{c}}\mathbf{\hat{A}}^{-1})^\mathrm{T}$. OK got it. Thanks! – uhoh Mar 23 '17 at 15:45

It turns out I made a dumb math error (a transposition, specifically), related to this math.stackexchange question. In short, I was computing $\mathbf{\lambda} = \mathbf{\hat{A}}^{-1}\mathbf{\hat{c}}$ instead of $\mathbf{\lambda}^\textrm{T} = \mathbf{\hat{c}}^\textrm{T}\mathbf{\hat{A}}^{-1}$.

It still produced what appeared to be valid answers in many cases — probably because my $\mathbf{\hat{A}}$ matrices were often symmetric, particularly for translation.

I do want to note that the values one uses for the different epsilons can make a large difference in some edge cases. So while it didn't solve my problem, it may come in handy for others. Use machine epsilon for your matrix inversion (or don't invert, if you can help it, since forward and back-substitution is more efficient). Use larger negative epsilon values to determine whether $\mathbf{x}$ is non-negative and for the delta-cost computation.

I hope this comes in handy for someone.

• +1 because it's always good practice to post your own solution to an unanswered question if you find it - it may come in handy to someone some day, and it can potentially avoid someone wasting time working on a solution that you don't need anymore. Even if it's a supplemental answer to another solution, you never know when someone else could find the additional insight helpful. – uhoh Apr 4 '17 at 8:18

It looks like the problem stems from this layout being ill-conditioned for numerical calculation.

As you wrote "almost symmetric about its center-of-gravity", the word "almost" is the red flag in any computational models.

In particular, the calculation probably attempts to determine pitch and yaw caused by the minuscule offset of the aft thrusters (instead of simply assuming they cancel each other, it calculates torque caused by first from the pair, the second from the pair, and only then sums these). This introduces almost zero torque; almost - but not quite. The value will be within floating point error range, these errors will accumulate and can completely trash your results.

Either make the configuration completely symmetric, or strongly asymmetric. Anything "almost" is asking for trouble.

• The length is 9, but the cg is at 5. How asymmetric does it need to be? – Doctor Mohawk Mar 23 '17 at 16:49
• @Dr.JohnnyMohawk: is it asymmetric in all axis, or just length-wise? Because it seems the rear thrusters are offset in the Z axis (vertical) by a really token amount. It must be either completely symmetric or considerably asymmetric in all axis, not just one. – SF. Mar 23 '17 at 17:34
• Also - I don't know if it's the root cause. I don't know your program or what and how it actually does, but I know an ill-conditioned problem when I see one. Still, sometimes the ill conditions may be canceled out, smoothed out or rendered unimportant by other means; so it might not be the cause. But if you say "15cm apart is okay, almost in line is not", that really looks like the case. Make the thrusters offset by 4cm and see if the problem persists. – SF. Mar 23 '17 at 18:23
• So, you were not wrong. The matrices were ill-conditioned, because of a math error I'd made (a transposition). Thanks for the tip. – Doctor Mohawk Apr 3 '17 at 16:31