Advanced Question Calculating optimal thruster choice...

Hlynkacg

Aspiring rocket scientist
Addon Developer
Tutorial Publisher
Donator
Joined
Dec 27, 2010
Messages
1,870
Reaction score
3
Points
0
Location
San Diego
What does "break" mean? What happens? I have some ideas if the main loop just continues to loop until it hits MaxIterations.

And continues to add thrusters a-la our DG loop conundrum

---------- Post added at 16:11 ---------- Previous post was at 16:08 ----------

Code:
bool AAPO::FCS::FindSolution ()
{
	// Get desired torque and normalize it.
	VECTOR3 tVector = DesiredTorque;		// Direction of torque
	double tMag = length (tVector);			// Magnitude of torque
	if (tMag > 0) normalise (tVector);		// Normalise vector

	// Do same for force.
	VECTOR3 fVector = DesiredForce;			// Direction of force
	double fMag = length (fVector);			// Magnitude of force
	if (fMag > 0) normalise (fVector);		// Normalise vector

	if (tMag + fMag <= 0) return false;		// Sanity Check: Is there a solution to find?

	// Debuggery
	sprintf (oapiDebugString(), "Desired Torque x %0.1f, y %0.1f, z %0.1f, Desired Force x %0.1f, y %0.1f, z %0.1f", DesiredTorque.x, DesiredTorque.y, DesiredTorque.z, DesiredForce.x, DesiredForce.y, DesiredForce.z);

	// Initialize variables
	VECTOR3 ResultantTorque = NULL_VECTOR;
	VECTOR3 ResultantForce = NULL_VECTOR;

	double tolerance = 0.1;
	int iterations = 0;
	int MaxIterations = thrusters.size ();
	bool solution = false;

	// Create a list of thrusters to be used 
	std::vector <std::pair <THRUSTER_DATA, double>> ThrusterList;

	// Get to work
	while ((solution == false) && (iterations < MaxIterations))
	{
		VECTOR3 TorqueRemaining = DesiredTorque - ResultantTorque;
		VECTOR3 ForceRemaining = DesiredForce - ResultantTorque;

		bool GoodTorque = (length (TorqueRemaining) < tolerance);
		bool GoodForce = (length (ForceRemaining) < tolerance);

		if ((GoodTorque) && (GoodForce)) solution = true;
		else // we need to iterate
		{
			VECTOR3 tr = NULL_VECTOR;
			VECTOR3 fr = NULL_VECTOR;

			if (length (TorqueRemaining) > 0)
			{
				tr = TorqueRemaining;
				normalise (tr);
			}

			if (length (ForceRemaining) > 0)
			{
				fr = ForceRemaining;
				normalise (fr);
			}

			// find all thrusters that positively contribute to desired torque
			int best = -1;
			double highscore = 0;
			std::vector<THRUSTER_DATA>::iterator i = thrusters.begin ();
			for (int i = 0; i < thrusters.size(); i++)
			{
				double score = WeighThruster (thrusters[i], DesiredTorque, DesiredForce);
				if (score > 0)
				{
					double lvl = score / length (thrusters[i].F);
					ThrusterList.push_back (std::make_pair (thrusters[i], lvl));
				}	
				i++;
			} 

			if (best != -1) 

			// Debuggery
			sprintf (oapiDebugString(), "%i Thrusters on list", ThrusterList.size ());
			//sprintf (oapiDebugString(), "Desired Torque x %0.1f, y %0.1f, z %0.1f, Desired Force x %0.1f, y %0.1f, z %0.1f, %i Thrusters on list", DesiredTorque.x, DesiredTorque.y, DesiredTorque.z, DesiredForce.x, DesiredForce.y, DesiredForce.z, ThrusterList.size ());

			VECTOR3 T = NULL_VECTOR;
			VECTOR3 F = NULL_VECTOR;

			// determine total torque and force of all thrusters on list.. 
			for (int i = 0; i < ThrusterList.size(); i++)
			{		
				//we don't care about pos, only the force/torque values
				F += ThrusterList[i].first.F * ThrusterList[i].second;
				T += ThrusterList[i].first.T * ThrusterList[i].second;
			}
	
			// Scale to magnitude of desired force and torque
			double tMod = 1;
			if (length (T) > tMag) tMod = tMag / length (T);
	
			double fMod = 1;
			if (length (F) > fMag) fMod = fMag / length (F);

			double mod = (fMag > tMag ? fMod : tMod);

			// Update list and resultants
			for (int i = 0; i < ThrusterList.size(); i++)
			{
				ThrusterList[i].second *= mod;
				ResultantTorque += ThrusterList[i].first.T * ThrusterList[i].second;
				ResultantForce += ThrusterList[i].first.F * ThrusterList[i].second;
			}

			// Debuggery
			sprintf (oapiDebugString(), "Desired Torque x %0.1f, y %0.1f, z %0.1f, Desired Force x %0.1f, y %0.1f, z %0.1f, Torque x %0.1f, y %0.1f, z %0.1f, Force x %0.1f, y %0.1f, z %0.1f, %i Thrusters on list",  DesiredTorque.x, DesiredTorque.y, DesiredTorque.z, DesiredForce.x, DesiredForce.y, DesiredForce.z, T.x, T.y, T.z, F.x, F.y, F.z, ThrusterList.size ());
		} // End iteration back to the top;

		iterations++;

		// Debuggery
		sprintf (oapiDebugString(), "%i iterations", iterations);

	} // End "while ((solution == false) && (iterations < MaxIterations))"


	// Note: neet to make some sort of scaling scheme to prevent excessivly high torques or forces from breaking the solver  

	// Use "ThrusterList" to generate firing commands
	//if ((solution) && (ThrusterList.size() > 0))
	if (ThrusterList.size() > 0)
	{
		//double simdt = oapiGetSimStep ();

		for (int i = 0; i < ThrusterList.size(); i++)
		{
			double lvl = ThrusterList[i].second;
			FCS::FireThruster (ThrusterList[i].first.th, lvl);
		}
		return true;
	}
	else return false; 	
}

Ok most recent code iteration.

Loop issue fixed but now yaw axis is non responsive and inverted control inputs are back.

Pretty sure it's just a swapped sign somewhere but not sure where.
 

meson800

Addon Developer
Addon Developer
Donator
Joined
Aug 6, 2011
Messages
405
Reaction score
2
Points
18
I'm not sure what is causing the inverted controls, but I do have a solution for the infinite loop situation.

Basically, it involves taking the desired force and torque vectors, normalizing them, and running them through the algorithm, then using that information to find the "maximum" torque/force values. Then, if the desired force/torque exceeds the maximum value, then scale the desired vectors before trying to find the "real" solution.

Here's what my idea would look like
Code:
//assume there is a function FindSolutionThrusterList
//which is most of your code in FindSolution, until it generates the ThrusterList
//but then just directly returns a std::vector <std::pair <THRUSTER_DATA, double>>
//signature:
//std::vector <std::pair <THRUSTER_DATA, double>> FindSolutionThrusterList(VECTOR3 force, VECTOR3 torque);


void FindSafeSolution(VECTOR3 desiredForce, VECTOR3 desiredTorque)
{
  //normalize the vectors and find the solution
  VECTOR3 desiredForceNorm = normalize(desiredForce);
  VECTOR3 desiredTorqueNorm = normalize(desiredTorque);
  std::vector <std::pair <THRUSTER_DATA, double>> thrusterList = FindSolutionThrusterList(desiredForceNorm, desiredTorqueNorm);
  
  //now, find the largest ratio that we can multiple the thruster list by without pushing one
  //of the thrusters beyond 100% thrust, and we can use this to calculate maximum force/torque
  double lowestRatio = 100000000; //some large number
  for(int i = 0; i < thrusterList.size(); i++)
  {
    double ratio = 1 / thrusterList[i].second; //sets ratio equal to the amount of times we can multiply
	//the thrust level by before hitting 100% thrust
	if (ratio < lowestRatio)
	  lowestRatio = ratio;
  }
  
  //find the maximum force/torque values
  double maximumForce = magnitude(totalThrusterForce(thrusterList)) * lowestRatio;
  double maximumTorque = magnitude(totalThrusterTorque(thrusterList)) * lowestRatio;
  
  VECTOR3 maximumForceVec, maximumTorqueVec;
  //NOTE: calculate the maximum force and torque VECTORS with some method
  //cap the desiredForce and desiredTorque at this value
  VECTOR3 finalDesiredForce = (magnitude(desiredForce) > magnitude(maximumForceVec)) ? maximumForceVec : desiredForce;
  VECTOR3 finalDesiredTorque = (magnitude(desiredTorque) > magnitude(maximumTorqueVec)) ? maximumTorqueVec : desiredTorque;
  FindSolution(finalDesiredForce, finalDesiredTorque)
 

Thorsten

Active member
Joined
Dec 7, 2013
Messages
785
Reaction score
56
Points
43
Do you actually need to solve the general problem?

I would assume any designer of a spacecraft would group thrusters into 'primitives' such that firing a primitive will give you either a pure force or a pure torque.

So rather than working with a list of arbitrary thrusters, you'd work with a pre-computed list of primitives characterizing your craft. And that makes the problem comparatively simple - for a force you just fire x-y-z primitives with thrust proportional to the dot product of their orientation and the desired force vector.

A lower, much simpler layer of the control code would then assign the commands to the primitives to individual thrusters.

But I think a two layer strategy will be much faster and more transparent.

Edit: As an afterthought: Not sure how much orbiter simulates this, but in reality there'd be more benefits of grouping thrusters to primitive operations.

One example is the use of propellant. Thrusters may not have a central propellant storage. Assume we need a force into z direction, and we have a thruster which provides it directly. A loop running through all thrusters will always fire the first one found. As a result, its propellant will be depleted early. A primitive z-thrust group would be defined in such a way that thrust is preferentially shared across different thrusters, leading to a more even consumption of propellant.

A similar case may be that the thrust is a non-linear function of propellant flow. In this case, firing two thrusters at 50% could be superior to firing one thruster at 100% throttle. Again, you would want to include this into a control system.
 
Last edited:

Hlynkacg

Aspiring rocket scientist
Addon Developer
Tutorial Publisher
Donator
Joined
Dec 27, 2010
Messages
1,870
Reaction score
3
Points
0
Location
San Diego
Do you actually need to solve the general problem?

Yes, at least partially, because neither the thruster configuration nor the center of gravity are fixed.

---------- Post added at 08:36 ---------- Previous post was at 08:21 ----------

Code:
  for(int i = 0; i < thrusterList.size(); i++)
  {
    double ratio = 1 / thrusterList[i].second; //sets ratio equal to the amount of times we can multiply
	//the thrust level by before hitting 100% thrust
	if (ratio < lowestRatio)
	  lowestRatio = ratio;
  }

I don't understand how this helps. If i'm reading this correctly the loop will have already iterated itself in to idiocy prior to this function taking effect. Otherwise it'd just be multiplying a bunch of zeros by zero.

Wouldn't it be better to figure out the theoretical max torque the thruster system can produce and then have something to the effect of...

Code:
if (length (Desired) > Max)
{
   normalise (desired);
   deisred *= max;
}

...at the beginning of the solver

---------- Post added at 09:27 ---------- Previous post was at 08:36 ----------

added a "scaling function" to the beginning of FindSolution()

---------- Post added at 09:51 ---------- Previous post was at 09:27 ----------

Ok progress update...

I found the transposed sign that was causing the control inversion and at this point it will successfully solve for torque inputs along the x and z axis. The y axis is hit or miss but I'm not sure why. I know the function is in serious need of commenting / clean-up, but i'm posting the complete body anyway in the interests of maintaining a visible change-log and getting more eyes on the potential issues.

Code:
bool AAPO::FCS::FindSolution ()
{
	// Scale inputs to possible (You requested 1,000,000,000,000 n/m of torque, that doesn't mean the thruster system is capable of providing it)
	for (int i = 0; i < 3; i++)
	{
		if (fabs (DesiredTorque.data[i]) > MaxTorque.data[i])					// if absolute value of desired torque is greater than maximum torque impulse...
		{
			double scale = MaxTorque.data[i] / fabs (DesiredTorque.data[i]);	// get scalar by dividing maximum torque impulse by desired torque.
			DesiredTorque.data[i] *= scale;										// multiply input value by scalar
		}

		if (fabs (DesiredForce.data[i]) > MaxTorque.data[i])					// Repeat process for force...
		{
			double scale = MaxForce.data[i] / fabs (DesiredForce.data[i]);
			DesiredForce.data[i] *= scale;
		}
	}

	// Get desired torque and normalize it.
	VECTOR3 tVector = DesiredTorque;		// Direction of torque
	double tMag = length (tVector);			// Magnitude of torque
	if (tMag > 0) normalise (tVector);		// Normalise vector

	// Do same for force.
	VECTOR3 fVector = DesiredForce;			// Direction of force
	double fMag = length (fVector);			// Magnitude of force
	if (fMag > 0) normalise (fVector);		// Normalise vector

	if (tMag + fMag <= 0) return false;		// Sanity Check: Is there a solution to find?

	// Debuggery
	sprintf (oapiDebugString(), "Desired Torque x %0.1f, y %0.1f, z %0.1f, Desired Force x %0.1f, y %0.1f, z %0.1f", DesiredTorque.x, DesiredTorque.y, DesiredTorque.z, DesiredForce.x, DesiredForce.y, DesiredForce.z);

	// Initialize variables
	VECTOR3 ResultantTorque = NULL_VECTOR;
	VECTOR3 ResultantForce = NULL_VECTOR;

	double tolerance = 0.1;
	int iterations = 0;
	int MaxIterations = thrusters.size ();
	bool solution = false;

	// Create a list of thrusters to be used 
	std::vector <std::pair <THRUSTER_DATA, double>> ThrusterList;

	// Get to work
	while ((solution == false) && (iterations < MaxIterations))
	{
		VECTOR3 TorqueRemaining = DesiredTorque - ResultantTorque;
		VECTOR3 ForceRemaining = DesiredForce - ResultantTorque;

		bool GoodTorque = (length (TorqueRemaining) < tolerance);
		bool GoodForce = (length (ForceRemaining) < tolerance);

		if ((GoodTorque) && (GoodForce)) solution = true;
		else // we need to iterate
		{
			VECTOR3 tr = NULL_VECTOR;
			VECTOR3 fr = NULL_VECTOR;

			if (length (TorqueRemaining) > 0)
			{
				tr = TorqueRemaining;
				normalise (tr);
			}

			if (length (ForceRemaining) > 0)
			{
				fr = ForceRemaining;
				normalise (fr);
			}

			// find all thrusters that positively contribute to desired torque
			int best = -1;
			double highscore = 0;
			std::vector<THRUSTER_DATA>::iterator i = thrusters.begin ();
			for (int i = 0; i < thrusters.size(); i++)
			{
				double score = WeighThruster (thrusters[i], DesiredTorque, DesiredForce);
				if (score > 0)
				{
					double lvl = score / length (thrusters[i].F);
					ThrusterList.push_back (std::make_pair (thrusters[i], lvl));
				}	
				i++;
			} 

			if (best != -1) 

			// Debuggery
			sprintf (oapiDebugString(), "%i Thrusters on list", ThrusterList.size ());
			//sprintf (oapiDebugString(), "Desired Torque x %0.1f, y %0.1f, z %0.1f, Desired Force x %0.1f, y %0.1f, z %0.1f, %i Thrusters on list", DesiredTorque.x, DesiredTorque.y, DesiredTorque.z, DesiredForce.x, DesiredForce.y, DesiredForce.z, ThrusterList.size ());

			VECTOR3 T = NULL_VECTOR;
			VECTOR3 F = NULL_VECTOR;

			// determine total torque and force of all thrusters on list.. 
			for (int i = 0; i < ThrusterList.size(); i++)
			{		
				//we don't care about pos, only the force/torque values
				F += ThrusterList[i].first.F * ThrusterList[i].second;
				T += ThrusterList[i].first.T * ThrusterList[i].second;
			}
	
			// Scale to magnitude of desired force and torque
			double tMod = 1;
			if (length (T) > tMag) tMod = tMag / length (T);
	
			double fMod = 1;
			if (length (F) > fMag) fMod = fMag / length (F);

			double mod = (fMag > tMag ? fMod : tMod);

			// Update list and resultants
			for (int i = 0; i < ThrusterList.size(); i++)
			{
				ThrusterList[i].second *= mod;
				ResultantTorque += ThrusterList[i].first.T * ThrusterList[i].second;
				ResultantForce += ThrusterList[i].first.F * ThrusterList[i].second;
			}

			// Debuggery
			sprintf (oapiDebugString(), "Desired Torque x %0.1f, y %0.1f, z %0.1f, Desired Force x %0.1f, y %0.1f, z %0.1f, Torque x %0.1f, y %0.1f, z %0.1f, Force x %0.1f, y %0.1f, z %0.1f, %i Thrusters on list",  DesiredTorque.x, DesiredTorque.y, DesiredTorque.z, DesiredForce.x, DesiredForce.y, DesiredForce.z, T.x, T.y, T.z, F.x, F.y, F.z, ThrusterList.size ());
		} // End iteration back to the top;

		iterations++;

		// Debuggery
		sprintf (oapiDebugString(), "%i iterations", iterations);

	} // End "while ((solution == false) && (iterations < MaxIterations))"


	// Note: neet to make some sort of scaling scheme to prevent excessivly high torques or forces from breaking the solver  

	// Use "ThrusterList" to generate firing commands
	//if ((solution) && (ThrusterList.size() > 0))
	if (ThrusterList.size() > 0)
	{
		//double simdt = oapiGetSimStep ();

		for (int i = 0; i < ThrusterList.size(); i++)
		{
			double lvl = ThrusterList[i].second;
			FCS::FireThruster (ThrusterList[i].first.th, lvl);
		}
		return true;
	}
	else return false; 	
}
 

Thorsten

Active member
Joined
Dec 7, 2013
Messages
785
Reaction score
56
Points
43
Yes, at least partially, because neither the thruster configuration nor the center of gravity are fixed.

Sure - but CoG doesn't change much per frame. Neither does thruster configuration. To compute slowly changing things inside a fast control logic is bad design.

Well - that sounds harsh, sorry. I come from real time rendering where 'per frame' really is a precious quantity.

But I deduce from the shortness of the response that you're not really interested in a two layer design, so I won't pursue this idea any further.
 

Hlynkacg

Aspiring rocket scientist
Addon Developer
Tutorial Publisher
Donator
Joined
Dec 27, 2010
Messages
1,870
Reaction score
3
Points
0
Location
San Diego
Sure - but CoG doesn't change much per frame. Neither does thruster configuration. To compute slowly changing things inside a fast control logic is bad design.

Well yes, which is why the list of "THRUSTER_DATA" structs that contain the relevant torque force position data etc... are saved as class variables and updated on "as needed" rather than on a "per frame" basis.

The issue here is that the torque and force needed/desired to maintain a given acceleration and attitude is some thing that changes quickly and as such needs to be updated on a per frame basis or close to it.

The purpose of FindSolution () is to review the saved thruster data and figure out how to achieve the desired torque or force if possible.

The issue I have with your "primitive" idea is that the primitives would need to be dynamically generated, otherwise I'd just be reinvention of Orbiter's standard THGROUP handler which kind of defeats the purpose of designing an integrated FCS.
 

meson800

Addon Developer
Addon Developer
Donator
Joined
Aug 6, 2011
Messages
405
Reaction score
2
Points
18
I don't understand how this helps. If i'm reading this correctly the loop will have already iterated itself in to idiocy prior to this function taking effect. Otherwise it'd just be multiplying a bunch of zeros by zero.

Wouldn't it be better to figure out the theoretical max torque the thruster system can produce and then have something to the effect of...

Code:

Code:
if (length (Desired) > Max)
{
   normalise (desired);
   deisred *= max;
}

...at the beginning of the solver
The loop normalized the desired torque/force before running it through this loop, so it only runs the solver for 1 N/m or N or torque/force. This should not make the solver freak out.

Then, it finds the maximum torque by doing the ratio thing.

A loop running through all thrusters will always fire the first one found
Actually, Hlynkacg has a weighting function which finds the best thruster to fire, not just the first one.
 

Hlynkacg

Aspiring rocket scientist
Addon Developer
Tutorial Publisher
Donator
Joined
Dec 27, 2010
Messages
1,870
Reaction score
3
Points
0
Location
San Diego
The loop normalized the desired torque/force before running it through this loop, so it only runs the solver for 1 N/m or N or torque/force. This should not make the solver freak out.

Then, it finds the maximum torque by doing the ratio thing.

Actually, Hlynkacg has a weighting function which finds the best thruster to fire, not just the first one.

Ok I get it now,

That actually has me thinking...

So one of the biggest issues is the throttle loop. solving for the normalized case and then scaling upwards might be the best solution. The problem is that I need a way of determining how much of a given T or F vector a specific thruster constitutes.

Something along the lines of the combined T of theses three thrusters is x,y,z. thruster 2 accounts for 2% of x, 18% if y and 10% of z.
 

Thorsten

Active member
Joined
Dec 7, 2013
Messages
785
Reaction score
56
Points
43
Well yes, which is why the list of "THRUSTER_DATA" structs that contain the relevant torque force position data etc... are saved as class variables and updated on "as needed" rather than on a "per frame" basis.

Well, that's not what I was talking about. You're running an optimization routine every time you have an updated target force from the control system. And that's completely unneeded. You need to run it only when the thruster data updates.



The issue here is that the torque and force needed/desired to maintain a given acceleration and attitude is some thing that changes quickly and as such needs to be updated on a per frame basis or close to it.

Yes, so if you know primitives for Fx,Fy and Fz and Mx, My and Mz, it's down to computing dot products, which is much much faster than doing an iterative loop as you plan to do.

You need to determine the primitives of course - but you can determine them every time the thruster table changes, not every tiime the required force updates.


The purpose of FindSolution () is to review the saved thruster data and figure out how to achieve the desired torque or force if possible.

Thanks, I can read C++ code, I know what it does. What it does primitive finding and doing the thruster commands all mixed up in an unintuitive way.

Your DeltaGlider loop conundrum exemplifies that. You're trying to solve two problems (find a desired force - find a desired torque) at the same time with a third problem, the individual thrust settings, and so you need to invest quite some time to avoid control logic instability - because the system tries to do several tasks at once, and that's what leads to the oscillation in the first place.

You have to realize that once you found primitives, the thrust control part is trivial. And you have no realize that finding primitives for Fx, Fy and Fz is a much easier problem than finding a solution for a general F, because Fx, Fy and Fz can use a coordinate system which is aligned with the axes of the thrusters as stored in your table (they don't need to be really x,y,z, since we're later doing dot products anyway - they just have to span a non-degenerate coordinate system).

The issue I have with your "primitive" idea is that the primitives would need to be dynamically generated, otherwise I'd just be reinvention of Orbiter's standard THGROUP handler which kind of defeats the purpose of designing an integrated FCS.

Running a chi^2 optimization loop inside the control logic is better than dynamically generating the primitives once the thruster table updates? Come on...
 

Hlynkacg

Aspiring rocket scientist
Addon Developer
Tutorial Publisher
Donator
Joined
Dec 27, 2010
Messages
1,870
Reaction score
3
Points
0
Location
San Diego
I'm sorry Thorsten, it seems I misinterpreted what you were suggesting.

On the first pass I essentially thought you were asking why I wasn't using statically defined thruster groups or the default THGROUP architecture and replied accordingly.

---------- Post added at 18:27 ---------- Previous post was at 15:34 ----------

All right, I'm thinking about this primitive idea and if anything breaking the throttle and thruster selection processes into thier own threads does simplify things.

So what are the base components of each primitive?

We've got a list of thrusters obviously. We need throttle level for each thruster, or rather a modifier to account for asymmetrical loads. And optionally we need the magnitude of torque or force applied for use in the actual control loop. For example if we know that the "Pitch_Neg" primitive produces 2000 n/m of torque at full power we know to scale the saved throttle levels by 0.2 to get a desired torque of "-400, 0, 0"

ETA:

This might be a good use / work-around for Meson's earlier algorithm...

because all the vectors are normalized before calculating the "remaining" vectors, the solution will always induce one Nm of torque per N of force (except in the case where one of the desired vectors is zero). Although the function that generates the burn list based on the solution can calculate different magnitudes, the burn list will always have a 1:1 torque to force ratio. I'm not sure how to solve this. :shrug:
 
Last edited:

Thorsten

Active member
Joined
Dec 7, 2013
Messages
785
Reaction score
56
Points
43
Well, I'm refining the idea as I go along, but I think I've got an algorithm sketched.

Basically, when you design thrusters on a spacecraft, you wouldn't add an arbitrary collection pointing into arbitrary directions, the normal case would be to arrange thrusters orthogonally (because if you don't, you waste propellant) and in opposite pairs. So the algorithm should make use of these properties and only contain fallback in case that's not doable.

Suppose you start out with a list of thrusters T1, ...Tn which contains max. F (as vector) and max M (as vector).

* First task is to go through the list and merge thrusters into subgroups where thrust axes are identical - so you normalize all thrust vectors, go through the list, open a different group if the normalized vector is different to one you already have or add to a group if the vector points the same direction.

-> For a normal thruster system, you should end up with six groups, and we can define these as +-x, +-y and +-z directions in the internal coordinate system

* second task - find out if these groups are orthonormal and opposing (i.e. whether the scalar product of any two unit vectors is either 0 or -1) - keep a list which groups thrust into opposite directions

(* third task - if you want to deal with thruster arrangements that are bad design, i.e. not symmetrically opposing, not orthonormal, you need to start from scratch.

So go through the list, find the largest (unnormalized thrust vector). That's the x-axis candidate. Normalize it. Go through the list again and find the thruster which minimizes the dot product. That's the y-axis candidate. Then go through the list again and find the vector v that minimizes t*x and t*y - that's the z-axis candidate. That also defines -x, -y and -z. Then do a clustering - take any vector remaining in the list - compute v*x, v*y and v*z and merge it to the direction which has the larges scalar product. Say that is the x-direction, so adding it will alter the x-direction a bit. Use that new x-direction to go through the list again, merge again.

-> You end up with a non-orthonormal coordinate system that groups thrusters into what points 'most' into some coordinate direction and you get the 'best' separated thrust axes.

-> orthonormalize it (that's mathematically well-defined, and done by taking one vector and solving the rest for vx * vy = 0 and vz = vx cross vy - you get linear combinations of vectors. Thrust from such a system will always be sub-optimal and waasteful because thrusters partially conter their effects in a non-ON system. That's inherent in the design.

Note that this really shouldn't be necessary unless a designer randomly attaches thrusters to a spacecraft without thinking what they're for. Point being - it can be solved, but this general case is really what makes the math difficult to deal with, and at the same time it wouldn't actually occur in reality).

* next task - find primitives. The beauty of being in an orthonormalized (ON) system is that primitives in to x direction can only involve thrusters in the x-group, primitives for Mx can only come from x and -x. So we never need to search the whole list.

For primitive force finding, add all thrusters in the +x group. This will in general give you some net torque Mx. If it is positive, find the thruster with the largest individual positive Mx, reduce its thrust till the net Mx is 0 or its thrust is zero. Work through the list until Mx is zero - that's the maximal thrust arrangement you can fire without creating a torque.

For torque primitive finding, add all +Mx generating thrusters in the +x and -x group. You will in general end up with a net force. Take the thruster with the smallest Mx and reduce thrust till the net force goes to zero.

Proceed through the other directions.


You end up with primitives Px, Py, Pz, ... for all directions, defined in terms of lists of thrusters to fire along with throttle settings for these. All this can be done 'as needed'.


* inside the control logic, just do dot products of primitives with your desired force-torque vectors, and you do it in the xyz coordinate system your thrusters are arranged in. So for any control input CF or CM, you do CF * P = CFx * Px + CFy * Py + CFz * Pz - which gives you a series of coefficients, e.g. CFx * Px = cx.

What you need to fire is cx * Px + cy * PY + cz * Pz.

This will tell you that e.g. you need to fire 0.3 * Px + 0.4 * Py. Pass that information to a lower layer - for instance Px might resolve as 1 * T1 + 1* T4 + 0.5 * T6 (where Ti are individual thrusters). Just multiply the 0.3 inside, so 0.3 Px resolves into 0.3 * T1 + 0.3 * T4 + 0.5 * 0.3 * T6.

-> It might happen that a thruster is part of several primitives and gets a pre-factor (throttle-setting) ct >1 - catch this case, return insufficient thrust into to the control system and normalize all throttle settings by 1/ct in order to get the maximum possible thrust.

-> We're done.

(* Well, there's the case of thrusters which might gimbal - that's potentially more complicated and I haven't worked that out yet...)
 

meson800

Addon Developer
Addon Developer
Donator
Joined
Aug 6, 2011
Messages
405
Reaction score
2
Points
18
Thorsten's solution is much more efficient than the current solution.

So go through the list, find the largest (unnormalized thrust vector). That's the x-axis candidate. Normalize it. Go through the list again and find the thruster which minimizes the dot product. That's the y-axis candidate. Then go through the list again and find the vector v that minimizes t*x and t*y - that's the z-axis candidate. That also defines -x, -y and -z. Then do a clustering - take any vector remaining in the list - compute v*x, v*y and v*z and merge it to the direction which has the larges scalar product. Say that is the x-direction, so adding it will alter the x-direction a bit. Use that new x-direction to go through the list again, merge again.

-> You end up with a non-orthonormal coordinate system that groups thrusters into what points 'most' into some coordinate direction and you get the 'best' separated thrust axes.
Hlynkacg's FindSolution already does this, so much of this code is already done. The weighting function already internally finds the thruster that will minimize the dot product.
 

Hlynkacg

Aspiring rocket scientist
Addon Developer
Tutorial Publisher
Donator
Joined
Dec 27, 2010
Messages
1,870
Reaction score
3
Points
0
Location
San Diego
Note that this really shouldn't be necessary unless a designer randomly attaches thrusters to a spacecraft without thinking what they're for. Point being - it can be solved, but this general case is really what makes the math difficult to deal with, and at the same time it wouldn't actually occur in reality).

Funny you should mention that, the need to solve this exact problem is why I'm going through the trouble of programing an FCS system in the first place rather than using orbiter's built-in thruster control groups.

picture.php


The CSM's designers weren't "randomly attaching thrusters" but additional constraints such as gas impingement and the routing of wiring and fuel lines resulted in a RCS configuration that was non-orthagonal. Add in the challenges of shifting CoG and things quickly get a bit complicated.

picture.php


picture.php


Likewise for the CM's asymmetrical configuration.

The LM's thruster system is the only one close to being "optimized" but that's only after the descent stage has been separated and ascent engine propellant expended. Otherwise we run into the same sort of problems caused by the SM's shifting CoG.
 

Thorsten

Active member
Joined
Dec 7, 2013
Messages
785
Reaction score
56
Points
43
Visually they look pretty orthogonal to me - they may not be 100%, but there's a lot of symmetry in there which can be exploited - for instance they do come in opposing pairs.

CoG shifts are not an issue - that's easy to account for once you have an orthogonal system. Getting orthogonal is the only tricky part.

Hlynkacg's FindSolution already does this, so much of this code is already done. The weighting function already internally finds the thruster that will minimize the dot product.

Yes, you're quite right - the algorithm can easily be used for the orthogonalization.

Funny you should mention that, the need to solve this exact problem is why I'm going through the trouble of programing an FCS system in the first place rather than using orbiter's built-in thruster control groups.

You have to forgive me, I know nothing what Orbiter already has or provides (nor am I interested in the actual application). I'm interested in algorithms, I found the problem you posted an interesting puzzle which I thought I should be able to work out, so I invested some time. So if I say something that's completely obvious to anyone coding for Orbiter, it's because I really don't know any better.
 
Last edited:

Hlynkacg

Aspiring rocket scientist
Addon Developer
Tutorial Publisher
Donator
Joined
Dec 27, 2010
Messages
1,870
Reaction score
3
Points
0
Location
San Diego
You have to forgive me, I know nothing what Orbiter already has or provides (nor am I interested in the actual application). I'm interested in algorithms, I found the problem you posted an interesting puzzle which I thought I should be able to work out, so I invested some time. So if I say something that's completely obvious to anyone coding for Orbiter, it's because I really don't know any better.

No worries, and I do appreciate the help. :tiphat:

Currently sketching out some code and will post when I have something worth noting.

---------- Post added at 13:42 ---------- Previous post was at 11:58 ----------

I have to confess that I've been messing with the code for 2 hours now and I'm still having trouble visualizing this orthagonalization process.

My first impulse is to writ a funtion that solves for a given normal and then do something along the lines of

Code:
x_pos = BuildThrusterPrimitive (Torque, X_AXIS);
x_neg = BuildThrusterPrimitive (Torque,-X_AXIS);
y_pos = BuildThrusterPrimitive (Torque, Y_AXIS);
etc...

but that's not what you're describing is it?
 

Thorsten

Active member
Joined
Dec 7, 2013
Messages
785
Reaction score
56
Points
43
There are several ways this could be done.

My first idea was to use a data-driven approach and let the thrust axes as specified in the thruster file decide what the coordinate system should be. That would work in the case you're interested in, because there is a high degree of symmetry.

For instance, you might end up with thrusters pointing along the primed directions

x' = x
y' = 0.087 * x + 0.996 * y
z' = z

i.e. x and z, y and z are orthogonal, but x and y are not - there's a 5 degree angle from the normal.

In this case, you'd orthogonalize by introducing the group

y' - 0.087 * x'

which fires always y' and x' thrusters together when a pure y thrust is asked for, store that into somewhere on a low level and just manage the pure y on the higher control levels. That would have the advantage that the results you get look 'simple', i.e. you can read off from the tables you get that there's 5 degree offset in the y-axis.

Now, since you already have an optimization routine, you don't necessarily need to care about finesse, so I think you can just impose an orthogonal coordinate system by hand.

Then you ask your optimizer:

* this is my x direction - what is the maximal thrust with zero torque I can get for it?
-> if your optimizer works as advertized, it'll return the primitive for your x-axis

* repeat for orthogonally imposed y and z

-> since the coordinate system you specified may not be the coordinate system inherent in the thruster axes, the lists of thrusters will in general be complicated. So even assuming your original spacecraft design had an orthogonal x', y' and z' system, in general the coordinate system you impose can be different by two angles, so you get

v' = R v where R is a general rotation matrix. So you couldn't necessarily interpret the results you find that way easily, but they would be correct nevertheless.

I guess given that you have the optimization loop already, this is the way I would go, so you would indeed do as your code fragment suggests, just specify any ON system by hand and ask the optimizer for the solution for that axis.

Oh, and I don't know if you noticed, but in CF * P = CFx * Px + CFy * Py + CFz * Pz I silently assumed that the control input is normalized to the maximally available thrust of the primitive - otherwise the physical units don't drop out and you don't end up with a number between 0 and 1 characterizing the throttle setting needed.
 

Hlynkacg

Aspiring rocket scientist
Addon Developer
Tutorial Publisher
Donator
Joined
Dec 27, 2010
Messages
1,870
Reaction score
3
Points
0
Location
San Diego
I'm still interested in Thorsten's proposal, I'm just not sure how to go about implementing it.

In the meantime I've just finished a major revision of the current code...

Code:
bool AAPO::FCS::FindSolution ()
{
	// Generate primitives
	THRUSTER_PRIMITIVE x, y, z;

	if (DesiredTorque.x < 0) x = GeneratePrimitive (TORQUE,-X_AXIS);
	else if (DesiredTorque.x > 0)  x = GeneratePrimitive (TORQUE, X_AXIS);

	if (DesiredTorque.y < 0) y = GeneratePrimitive (TORQUE,-Y_AXIS);
	else if (DesiredTorque.y > 0)  y = GeneratePrimitive (TORQUE, Y_AXIS);

	if (DesiredTorque.z < 0) z = GeneratePrimitive (TORQUE,-Z_AXIS);
	else if (DesiredTorque.z > 0)  z = GeneratePrimitive (TORQUE, Z_AXIS);

	// Calculate throttle levels
	double x_lvl = fabs (DesiredTorque.x) / x.mag;
	double y_lvl = fabs (DesiredTorque.y) / y.mag;
	double z_lvl = fabs (DesiredTorque.z) / z.mag;

	// Generate firing commands
	double deadzone = 0.01;
	if (x_lvl > deadzone) FireThrusterGroup (x, x_lvl);
	if (y_lvl > deadzone) FireThrusterGroup (y, y_lvl);
	if (z_lvl > deadzone) FireThrusterGroup (z, z_lvl);
	
	return true; // Place holder till I can build a checker ensure that induced force and torque match desired)
}

Obviously the current function only solves for torque but the important thing is that it works with a minimal amount of over & undershoot.

each "Primitive" is a struct containing one of our "std::vector <std::pair <THRUSTER_DATA, double>>" ThrusterLists and the magnitude of force or torque that it produces, ".mag".

The thruster selection process itself has been moved to "GeneratePrimitive ()" and still needs refinement. Eventually I'll make the primitives class variables so that they don't have to be generated each frame, for the moment this is just a proof of concept.

here is the current body of "GeneratePrimitive ()"...

Code:
THRUSTER_PRIMITIVE AAPO::FCS::GeneratePrimitive (solve4 mode, VECTOR3 vector)
{
	// Initialize variables
	std::vector <std::pair <THRUSTER_DATA, double>> ThrusterList;

	THRUSTER_PRIMITIVE	output;
	VECTOR3				ResultantTorque = NULL_VECTOR;
	VECTOR3				ResultantForce = NULL_VECTOR;
	VECTOR3				tVector = NULL_VECTOR;
	VECTOR3				fVector = NULL_VECTOR;
	
	int iterations = 0;
	int maxiterations = thrusters.size ();
	bool solution = false;

	// Normalise Vector
	if (length (vector) > 0) normalise (vector);

	// Get to work
	switch (mode)
	{
	case TORQUE: // Solve for torque
		tVector = vector; // Direction of torque

		// Build a list of the thrusters that positively contribite to our desired torque and force vectors
		for (int i = 0; i < thrusters.size(); i++)
		{
			double score = WeighThruster (thrusters[i], tVector, fVector);
			if (score > 0)
			{
				double lvl = score / length (thrusters[i].F);					// Normalize score to get throttle level
				ResultantTorque += thrusters[i].T * lvl;						// Add thruster's torque to result 
				ResultantForce += thrusters[i].F * lvl;							// Add thruster's force to result
				ThrusterList.push_back (std::make_pair (thrusters[i], lvl));	// Add thruster to list
			}	
		} 

		// Refine solution
		while ((solution == false) && (iterations < maxiterations))
		{
			iterations++;
		}
		break;

	case FORCE: // Solve for force
		fVector = vector;								// Direction of force

		// *COMING SOON*

		break;
	}

	// Update list
	for (int i = 0; i < ThrusterList.size(); i++) output.ThrusterList.push_back (std::make_pair (ThrusterList[i].first, ThrusterList[i].second));

	if (mode == TORQUE) output.mag = length (ResultantTorque);
	else output.mag = length (ResultantForce);

	output.type = mode;
	output.nThrust = ThrusterList.size();

	return output;
}
 

Thorsten

Active member
Joined
Dec 7, 2013
Messages
785
Reaction score
56
Points
43
I'm still interested in Thorsten's proposal, I'm just not sure how to go about implementing it.

Can you outline where you're stuck?
 

Hlynkacg

Aspiring rocket scientist
Addon Developer
Tutorial Publisher
Donator
Joined
Dec 27, 2010
Messages
1,870
Reaction score
3
Points
0
Location
San Diego
Apologies for letting this sit but once again obligations in meat-space have been keeping me busy.

Anyway, the area that I'm having trouble with is the orthaginalization and procedural assignment of the thruster groups. My current solution (see code above) works after a fashion but also introduces certain inefficiencies as illustrated in the following image...

picture.php


Whats happening in this picture is that a Pitch input is resulting in +Y thrust command to the aft RCS quads and a -Y command for the forward quads. At the same time a roll input is resulting in a +Y thrust command to the starboard RCS quads and a -Y command to the port.

The result is that the two opposing quads, port aft, and starboard fwd are trying to go in both directions at once.

I understand the concept behind thorsten's algorithm but the structure of the required code and procedure for implementing it eludes me.
 
Last edited:

Hlynkacg

Aspiring rocket scientist
Addon Developer
Tutorial Publisher
Donator
Joined
Dec 27, 2010
Messages
1,870
Reaction score
3
Points
0
Location
San Diego
Ok bad news.

I added the descent module back onto the LM and the result is that the thruster grid is asymmetrical relative to the CoG and the function no-longer functional.

well it functions but does not complete it's objectives.

...so back to the drawing board.

I'm re-assessing Darren's Simplex algorithm but I could use some extra eyes on it. For one thing he has this "RCS_ENUMERATOR" class that I do not understand the exact purpose of as the only enumerators I am familiar with in C++ are the "enum XXX" variety.

I could also use some help parsing his tableau.cpp. I understand that he is essentially writing a procedural version of
but I'm having difficulty following what goes where on the code end, and what parts I actually need.
 
Top