General Question Pointing a ship towards a specific position

evilfer

New member
Joined
Dec 9, 2010
Messages
53
Reaction score
0
Points
0
Hi,

I'm trying to build an MFD that fires RCS thrusters in order to point a ship towards a specific position.

I have something what works, but my current implementation overshoots the desired direction. Since the Orbiter autopilots do a nice work at this, I'd like to ask about a proper algorithm to set the RCS thruster levels in order to reach the vessel's attutude without overshooting and waisting fuel.

Here's what I currently do:

Code:
VECTOR3 wantedDirection;
VECTOR3 currentDirection; // (in cartesian [x, y, z] with origin [0, 0, 0] at the ship position, calculated from AROT from vessel information.
VECTOR3 angularVelocity; // = AVEL from vessel information.

VECTOR3 cross = currentDirection x wantedDirection;
VECTOR3 neededAcceleration = 0.1 * cross - angularVelocity; // so that it's not too fast

vessel->SetAttitudeRotLevel(neededAcceleration);

I'm a little lazy to do the maths, and actually don't feel comfortable working with rotations. Any idea/pointer to a better algorithm than this?

Thank you!

BTW, I love the orbitersdk apis :)
 

Wishbone

Clueless developer
Addon Developer
Joined
Sep 12, 2010
Messages
2,421
Reaction score
1
Points
0
Location
Moscow
While I'm about to leave for the weekend...
Read the following source codes: Attitude MFD, Space shuttle ultra, Attitude RCS (Lua script) and your humble servant's NTR core propulsion stage.
Read up on quaternions, Euler matrices and frame transformations. SPICE has nice background docs, so it is closer than you think.

You'd also need some reading on optimal control theory. Try rummaging through NTRS, there are Apollo and Space Shuttle guidance docs.

Wertz and Larson is the book to read as well.

P.S. Take what I write with a huge grain of salt. It is not for nothing that I put the word "clueless" in my title :)
 
Last edited:

evilfer

New member
Joined
Dec 9, 2010
Messages
53
Reaction score
0
Points
0
While I'm about to leave for the weekend...
Read the following source codes: Attitude MFD, Space shuttle ultra, Attitude RCS (Lua script) and your humble servant's NTR core propulsion stage.
Read up on quaternions, Euler matrices and frame transformations. SPICE has nice background docs, so it is closer than you think.

You'd also need some reading on optimal control theory. Try rummaging through NTRS, there are Apollo and Space Shuttle guidance docs.

Wertz and Larson is the book to read as well.

P.S. Take what I write with a huge grain of salt. It is not for nothing that I put the word "clueless" in my title :)

I wasn't sure AttitudeMFD implemented that functionality after quickly reading its description, I'll check it, as well as the other addons you indicate.

Your answer actually temtped me to look into this :D. Other physic issues motivate me more, but rotations, angular velocities and momentums are really not my favourite. Anyway I'll see if I get some time to read a bit on this...

Thanks!

PD: Still, I have the feeling that there is a simple algorithm that does the job! :D
 

Enjo

Mostly harmless
Addon Developer
Tutorial Publisher
Donator
Joined
Nov 25, 2007
Messages
1,665
Reaction score
13
Points
38
Location
Germany
Website
www.enderspace.de
Preferred Pronouns
Can't you smell my T levels?
I don't quite understand. You said that you have the desired rotation properly calculated, right? Does the same apply to the adjustment rotation? Because what bothers me, is that I'd rather do:

e = length(vTarget - vCurrent);

So now assuming that you want to lessen the error signal you'd do something like:

acc = e * Kp

where Kp = 0.1, giving you a P controller.
Adding module of velocity towards the target (future position anticipation) to the equation would give you:

acc = e * Kp + v * Kd

where v = length(velocity) and Kd = 0.05, or whatever, giving you a PD controller.
You could also add past position data as an integral of e through time of the controller operation, giving you a PID controller, but Artlav, facing the same problems, told that he uses PD only

Having this required module of angular acceleration, you need to point it somehow to the target. I'm not sure what I'd do, but maybe I'd do the following:

ve = (vTarget - vCurrent) / length(vTarget - vCurrent);

which gives you a normalized target vector - pure direction, waiting to be multiplied by some meaningful value, which is acc, so:

vacc2BUsed = ve * acc;

Maybe you could try this PID autopilot:
[ame="http://www.orbithangar.com/searchid.php?ID=754"]PIDPilot[/ame]
And for some after hours read:
[ame="http://en.wikipedia.org/wiki/PID_controller"]PID controller - Wikipedia, the free encyclopedia[/ame]

Note that the velocity in my equation assumes that you accelerate towards the target, which may not necessarily be true. Maybe you could do the following:

Code:
VECTOR3 diffPos = vTarget - vCurrent;
static VECTOR3 prevRelPos = diffPos;
VECTOR3 currRelPos = diffPos;
v = length(currRelPos - prevRelPos) / oapiGetSimStep();
VECTOR3 prevRelPos = diffPos;

oapiGetSimStep() is the frame passed by this frame, not the MFD refresh time! I assume that you execute your code from plugin's context, not MFD?
You also have to take care of the sign of v somehow. Maybe you can try it without it first.

Let me know how it suits you. I want to play around with it someday too.
 
Last edited:

evilfer

New member
Joined
Dec 9, 2010
Messages
53
Reaction score
0
Points
0
Hi Enjo,

I do need some after hours read to understand your post :D, but it seems very interesting. I have to check the PD and PID controllers.

Currently I don't user rotation matrices, only cartesians vectors that indicate the direction. This means that I ignore roll. You can see how I get the current direction vector in this format in the code below (lines 6-8). This is the same format as the desired direction vector, which in the code is: r_current->getThrustVector().

For your method, do these vectors work?

Yesterday night I dedicated the energy that I should have kept for today's work to implement another method. Here's the code:

Code:
static const double maxvel = .1, angleThreshold = .5, thrustK = 20,
		angleerrormargin = .0000000001;

static VECTOR3 arot, cross, avel, need;
pV->GetGlobalOrientation(arot);
pV->GetAngularVel(avel);

m_currentDirection.x = sin(arot.x) * sin(arot.z) + cos(arot.x) * sin(arot.y) * cos(arot.z);
m_currentDirection.y = -sin(arot.x) * cos(arot.z) + cos(arot.x) * sin(arot.y) * sin(arot.z);
m_currentDirection.z = cos(arot.x) * cos(arot.y);

cross = crossp(m_currentDirection, r_current->getThrustVector());
double crossMod = length(cross);
m_currentManeuverAngleDiff = asin(length(crossMod));

if (m_currentManeuverAngleDiff > angleerrormargin) {

	double desiredAngularVelocity = 
		m_currentManeuverAngleDiff < angleThreshold ? 
			maxvel * (m_currentManeuverAngleDiff / angleThreshold) : 
			maxvel;

	m_currentAngularVelocity = length(avel);

	need = cross * (desiredAngularVelocity / crossMod) - avel;
	need *= thrustK;
} else {
	need.x = need.y = need.z = 0;
}

pV->SetAttitudeRotLevel(need);

Now this works reasonably well I'd dare to say. Starting from a 45º angle error, it reduces it under the constant angleerrormargin in about 100 seconds, without overshooting.

This constant angleerrormargin is ~5e-9degrees. When the angle error is below that value, the MFD turns the RCS off. One problem is: if I lower this value, the operation becomes unstable: after reaching a very low angle error, the MFD appartently overshoots the desired direction and the error starts growing exponentially until very large values (> 10degrees). The other constants should be decided according to ship properties (RCS max angular acc and moment of inertia at least), but I cannot do this now.

About oapi, I derived oapi::Module and overwrote its clbkPreStep(double simt, double simdt, double mjd) method. This method calls the MFD code to let it update itself :). All the calculations are done by the MFD object, while the Module object acts only as a connector between Orbiter events and the MFD. Is this a good practice for MFDs?

Thanks,
evilfer

- Edit: about the stability issue that I mentioned before, I found that it's not because of the algorithm, but due to a wrongly configured testing scenario file and a SUBorbital trajectory :facepalm:. The algorithm appears to be stable as long as the vessel remains in space...
 
Last edited:

Enjo

Mostly harmless
Addon Developer
Tutorial Publisher
Donator
Joined
Nov 25, 2007
Messages
1,665
Reaction score
13
Points
38
Location
Germany
Website
www.enderspace.de
Preferred Pronouns
Can't you smell my T levels?
Currently I don't use rotation matrices, only cartesians vectors that indicate the direction. This means that I ignore roll. You can see how I get the current direction vector in this format in the code below (lines 6-8). This is the same format as the desired direction vector, which in the code is: r_current->getThrustVector().

For your method, do these vectors work?

Sure they will. For roll, I'd wait until the ship is in proper position, or even better, roll first, and then yaw / pitch in one axis only - saves some calculations for a good start. For the second method, first try to get the current relative bank (roll) to the target, by printing it on MFD. In Launch MFD I do the following (see inputBank):

Code:
VECTOR3 velRot;
data->GetVessel()->GetAngularVel( velRot );
double inputBank = atan2( targetVector2D.x , targetVector2D.y );

fuzzyControllerAPBankAtm.SetInputPosition( targetVector2D.len() );
fuzzyControllerAPBankAtm.SetInputVelocity( -velRot.z );
fuzzyControllerAPBankAtm.SetInputBank( inputBank );
fuzzyControllerAPBankAtm.Evaluate();
double bankAdjust = fuzzyControllerAPBankAtm.GetOutput();

data->GetVessel()->SetAttitudeRotLevel( 2, -bankAdjust );

The fuzzyController is the thing which you'd replace with PD controller. As you can see, the hardest thing is to get the targetVetor (actually it's the adjustmentVector). Perhaps you could get it by doing:

adjustmentVector = target - current ?

Please explain why you are crossing these vectors? After all crossing vectors gives a perpendicular vector, which is by no way towards the target :)

About oapi, I derived oapi::Module and overwrote its clbkPreStep(double simt, double simdt, double mjd) method. This method calls the MFD code to let it update itself :). All the calculations are done by the MFD object, while the Module object acts only as a connector between Orbiter events and the MFD. Is this a good practice for MFDs?

As long as you're testing one thing, it's OK, but not in the long run. Perhaps you could use this template, which would give you additional multiple vessels support? By the way, I'm glad that your oapi::Module callbacks work, because my don't. Maybe you could have a look into that template, after you finish of course, and tell me what's wrong there, or paste code of your module interface and methods implementations?

And my last word is that I'm glad that something works for you, but by using a proper PD controller, you wouldn't need to check any margins. When you're close to the target, the adjustments would automagically degrade to values close to 0.
 
Last edited:

evilfer

New member
Joined
Dec 9, 2010
Messages
53
Reaction score
0
Points
0
Please explain why you are crossing these vectors? After all crossing vectors gives a perpendicular vector, which is by no way towards the target :)

I cross the vectors in order to get the axis along which the vessel should rotate. For instance, let's imagine my ship is facing X (1,0,0), and I want tha have it facing Y (0,1,0). The cross product results in (0,0,1), or Z, which is the rotation axis, that is, the axis around which I must rotate the vessel to change its orientation. I'm sorry if this is not well expressed, I hope it's understandable though.

The nice thing is that the Orbiter API seems to like this perpendicular vector. If I pass the vector to the method void VESSEL::SetAttitudeRotLevel ( const VECTOR3 & th ) const , then the RCS thrusters will provoke a rotation exactly around the th axis.

In my code, before passing the cross vector to this method I first compare it with the current angular velocity (which also seems to represents a rotation axis), in order to adjust for a proper angular velocity (closing to 0 as the ship gets to the desired attitude).

As long as you're testing one thing, it's OK, but not in the long run. Perhaps you could use this template, which would give you additional multiple vessels support?

I have to look into that. Multivessel is not currently in my todo list, but I'll consider it.

By the way, I'm glad that your oapi::Module callbacks work, because my don't. Maybe you could have a look into that template, after you finish of course, and tell me what's wrong there, or paste code of your module interface and methods implementations?

Umm, I can only notice that in addition to oapiRegisterMFDMode(...), I also use the method void oapiRegisterModule ( oapi::Module * module ):

Code:
DLLCLBK void InitModule(HINSTANCE hDLL) {
	module = new TrajectoryModule(hDLL);
	oapiRegisterModule(g_module);
	...
	mode = oapiRegisterMFDMode(spec);
}

Might be it :).


And my last word is that I'm glad that something works for you, but by using a proper PD controller, you wouldn't need to check any margins. When you're close to the target, the adjustments would automagically degrade to values close to 0.

This philosophy could be applied to all my development of this MFD and the Trajectories simulator :). I guess it may be considered stupid, but since this is a hobby, my approach is 1. Try something; 2. Realize it doesn't work. 3. Study the thing and (maybe) 4. Improve it.

Now I have to add controllers to the study queue!
 

Enjo

Mostly harmless
Addon Developer
Tutorial Publisher
Donator
Joined
Nov 25, 2007
Messages
1,665
Reaction score
13
Points
38
Location
Germany
Website
www.enderspace.de
Preferred Pronouns
Can't you smell my T levels?
I cross the vectors in order to get the axis along which the vessel should rotate. For instance, let's imagine my ship is facing X (1,0,0), and I want tha have it facing Y (0,1,0). The cross product results in (0,0,1), or Z, which is the rotation axis, that is, the axis around which I must rotate the vessel to change its orientation. I'm sorry if this is not well expressed, I hope it's understandable though.

The nice thing is that the Orbiter API seems to like this perpendicular vector. If I pass the vector to the method void VESSEL::SetAttitudeRotLevel ( const VECTOR3 & th ) const , then the RCS thrusters will provoke a rotation exactly around the th axis.
That's fantastic, indeed! Having this, you can include the roll independently, because at each frame, OrbiterCore will know where to point the calculated velocity.


I have to look into that. Multivessel is not currently in my todo list, but I'll consider it.
Multivessel is a pro-bono feature, resulting from a proper code design, whose other purpose (optional) is to allow for communication between oapi::Module and the MFD with as little global / visible variables as possible. There's nothing "TO DO" there, just derive 2 subclasses :)

Umm, I can only notice that in addition to oapiRegisterMFDMode(...), I also use the method void oapiRegisterModule ( oapi::Module * module ):
Yep, this would be it. Thanks! :cheers:

Now I have to add controllers to the study queue!

But do you realize, that you don't even need to read that long Wikipedia article to make it work? In fact, what you're doing is implementing a controller. It's just one that I've never seen in my life :D. Making it a PD would be very elegant. Before I finally leave it all up to you, I'd like to propose to combine what you and I've wrote:

Code:
VECTOR3 cross = crossp(m_currentDirection, r_current->getThrustVector());
double error = length(r_current->getThrustVector() - m_currentDirection);
double Kp = 0.1;
double Kd = 0.1;
double needMod = error * Kp + m_currentAngularVelocity * Kd;
VECTOR3 need = needMod * cross;
pV->SetAttitudeRotLevel(need);

All you need to do now is play around with the K's, probably make them dependent on some ship's characteristics. I've also found a real PID goodie by Face. Attached. Whenever you want to try out the PID, let me know if it works. I'm intrigued :)
Good luck.
 

Attachments

  • PidMFD_0.1.zip
    47.7 KB · Views: 6

evilfer

New member
Joined
Dec 9, 2010
Messages
53
Reaction score
0
Points
0
Multivessel is a pro-bono feature, resulting from a proper code design, whose other purpose (optional) is to allow for communication between oapi::Module and the MFD with as little global / visible variables as possible. There's nothing "TO DO" there, just derive 2 subclasses :)

Very convincing, specially the second purpose (for my case). Where can I get the code, with the callback fix? :)

But do you realize, that you don't even need to read that long Wikipedia article to make it work? In fact, what you're doing is implementing a controller. It's just one that I've never seen in my life :D. Making it a PD would be very elegant. Before I finally leave it all up to you, I'd like to propose to combine what you and I've wrote:

Again, very convincing. I'll make it a PD controller, using your code. Later some day I'll check optimal control theory :D.


I've also found a real PID goodie by Face. Attached. Whenever you want to try out the PID, let me know if it works. I'm intrigued :)
Good luck.

I checked the PidMFD code. Pid.cpp shows that two old error values are kept in memmory (e1 = error[n-1] and e2 = error[n-2]), but they are not used I think. Rather, the integral value is calculated simply as error*timeStepDuration (not really an integration right?).

I don't like the idea of using error values from previous iterations, but I can test the integral value error*timeStepDuration. Does this value have any physical meaning?

The code for the (pseudo) PID would be:

Code:
VECTOR3 cross = crossp(m_currentDirection, r_current->getThrustVector());
double error = length(r_current->getThrustVector() - m_currentDirection);
double Kp = 0.1;
double Kd = 0.1;
double Ki = something;
double needMod =      Kp * error + 
                      Kd * m_currentAngularVelocity +
                      Ki * timeStepDuration * error;
VECTOR3 need = needMod * cross;
pV->SetAttitudeRotLevel(need);

And, finally, a question that must be understood as a product of the late hour: why not a PID2D (2D = second derivate) controller?

--- Edit ----

Ops, another thing. In addition to the error calculation function, I can think of a second:

Code:
double error = length(r_current->getThrustVector() - m_currentDirection);
double error = asin(length(cross));

Would the second line be against the PID controller design? The second error value is actually the angle difference between r_current->getThrustVector() and m_currentDirection, so I think it's more adequate.
 
Last edited:

Enjo

Mostly harmless
Addon Developer
Tutorial Publisher
Donator
Joined
Nov 25, 2007
Messages
1,665
Reaction score
13
Points
38
Location
Germany
Website
www.enderspace.de
Preferred Pronouns
Can't you smell my T levels?
Very convincing, specially the second purpose (for my case). Where can I get the code, with the callback fix? :)

Always here: http://sourceforge.net/projects/enjomitchsorbit/
You can go to Code -> SVN Browse and from there you can Download GNU Tarball, if you don't use Tortoise SVN that is. On the main page, there's also a downloadable html documantation. Beware - it's LGPLed. so you either share your code, or link the lib dynamically. I'll be preparing a dynamic lib later this week. Until then (whenever it happens), you can experiment with it as you like. The reason why I'm licensing it this way is that I want people to at least mention that they use my and Computerex' work while having multiple vessels support, so others can have it too.

Again, very convincing. I'll make it a PD controller, using your code. Later some day I'll check optimal control theory :D.

I forgot to mention, that you need to limit the control signal somehow, that the resulting needMod is always in range <-1, 1> by dividing it by some maxAccel, obtained as a ship's characteristic, through oapi. Notice that we're not steering by velocity here, but acceleration.
I don't remember how I calculated it for my controllers, but I've found this in FuzzyControllerAP:
Code:
double FuzzyControllerAP::m_statDeltaGliderRefRotAcc = 5000.0 / 11000.0; // Yaw thrust / empty mass
So you need to get thrust of your rotation group. You could use:

Code:
DWORD GetGroupThrusterCount (THGROUP_TYPE thgt) const
THRUSTER_HANDLE GetGroupThruster (THGROUP_TYPE thgt, DWORD idx) const
I checked the PidMFD code. Pid.cpp shows that two old error values are kept in memory (e1 = error[n-1] and e2 = error[n-2]), but they are not used I think. Rather, the integral value is calculated simply as error*timeStepDuration (not really an integration right?).

I don't like the idea of using error values from previous iterations, but I can test the integral value error*timeStepDuration. Does this value have any physical meaning?

I agree. It doesn't look like an integration. Integration would be:
Code:
IntegratorDiscrete::IntegratorDiscrete()
{
   Reset();
}

void IntegratorDiscrete::Reset()
{
   m_integ = m_ePrev = 0;
}

double IntegratorDiscrete::GetIntegral()
{
    return m_integ;
}

double Integrator::IntegrateTrapezoidal(double e, double dt) {
   m_integ += (e + m_ePrev) * dt / 2.0
   m_ePrev = e;
   return m_integ;
}
... based on equation for area of trapezoid: http://pl.wikipedia.org/wiki/Trapez#Trapez_prostok.C4.85tny
Rotate it 90* left, and it will resemble your error in time. So actually the code would look like:

Code:
double Ki = something;
double needMod =      Kp * error + 
                      Kd * m_currentAngularVelocity +
                      Ki * m_integrator.IntegrateTrapezoidal(error, dt);
[EDIT] It's possible, that your Kd should be negative, because the faster you are proceeding, the lesser the power should be. Noticed any similarities to your code? :) You assumed there that Kd = -1. Above you have a more formal controller's description. Even more formally, the Kd should stay positive when you apply an assumption, that positive velocity means moving away from target and negative means towards target.

And, finally, a question that must be understood as a product of the late hour: why not a PID2D (2D = second derivate) controller?

I don't know. Keep it simple. PIDs are widely used. I've never heard of PID2D.

Ops, another thing. In addition to the error calculation function, I can think of a second:

Code:
double error = length(r_current->getThrustVector() - m_currentDirection);
double error = asin(length(cross));
Would the second line be against the PID controller design? The second error value is actually the angle difference between r_current->getThrustVector() and m_currentDirection, so I think it's more adequate.

No, it would definitely not be against the PID design, if it works the way it's supposed to work. What's needed is a linear (I think) decrease of the error value towards 0, as you reach your desired orientation. In fact, if the code works, it gives you an opportunity to limit the maximal error up to +/- 180*. It's an interesting idea, since you already operate on one axis, thanks to SetAttitudeRotLevel(VECTOR3 &). One problem would be that you use asin, which is limited to (-PI/2, PI/2). You want (-PI, PI) to properly operate on all quarts. Therefore atan2 is usually preferred.

To get an angle between 2 vectors, in Launch MFD I do the following:

Code:
double VectorMath::dot(const Vect3 & a, const Vect3 & b )
{
  return a.x*b.x + b.y*a.y + a.z*b.z;
}

Vect3 VectorMath::cross(const Vect3 & a, const Vect3 & b )
{
  Vect3 v;
  v.x = a.y*b.z - b.y*a.z;
  v.y = a.z*b.x - b.z*a.x;
  v.z = a.x*b.y - b.x*a.y;
  return v;
}

double VectorMath::angle(const Vect3 & a, const Vect3 & b )
{
    Vect3 aNorm = a.norm();
    Vect3 bNorm = b.norm();

    return atan2( cross(aNorm, bNorm).len(), dot(aNorm, bNorm) );
}

double Vect3::len() const
{
    return sqrt(x * x + y * y + z * z );
}

Vect3 Vect3::norm() const
{
    return Vect3(*this) /= len();
}
The (reusable) classes for that are in Math and Systems folders. You'd only need to make a converter from VECTOR3 to my Vect3 and vice versa. You can mimic the SystemsConverter in Systems for that.

PS. Here's another thread I've found about PID.
 
Last edited:

evilfer

New member
Joined
Dec 9, 2010
Messages
53
Reaction score
0
Points
0
Always here: http://sourceforge.net/projects/enjomitchsorbit/
You can go to Code -> SVN Browse and from there you can Download GNU Tarball, if you don't use Tortoise SVN that is. On the main page, there's also a downloadable html documantation. Beware - it's LGPLed. so you either share your code, or link the lib dynamically. I'll be preparing a dynamic lib later this week. Until then (whenever it happens), you can experiment with it as you like. The reason why I'm licensing it this way is that I want people to at least mention that they use my and Computerex' work while having multiple vessels support, so others can have it too.

I published my code also in SF under GPL :). I think it's easier to add your source code files to my own. In this case, if I upload the code including your files to SF, would it be ok? (copyright and licences would be of course kept).


About the controller, it's now a PID. Currently the error is a VECTOR3, calculated in this way:

Code:
	m_error = crossp(m_currentDirection, m_targetDirection);

	double sin_angle = length(m_error);
	double cos_angle = dotp(m_currentDirection, m_targetDirection);
	m_currentangle = atan2(sin_angle, cos_angle);

	m_error *= m_currentangle / sin_angle;

That is, it's a vector whose direction is the rotation axis, and its magnitude is the angle difference between current and target attitudes. BTW, than you for the tip on asin and atan2!!


I forgot to mention, that you need to limit the control signal somehow, that the resulting needMod is always in range <-1, 1> by dividing it by some maxAccel, obtained as a ship's characteristic, through oapi. Notice that we're not steering by velocity here, but acceleration.
I don't remember how I calculated it for my controllers, but I've found this in FuzzyControllerAP:

I'm not limiting the signal and it seems to work. In addition to Kp and Kd, the key seems to be a multiplier that converts the result of P+I+D into the actual vector passed as the vessel RCS levels. If this multiplier is too low, the operation is slow; if it is too large, it overshoots. This multiplier will be related to the maxAccel :).


I agree. It doesn't look like an integration. Integration would be:

I tried an integrator that keeps N old values and timesteps. However, I'm not happy with it. This integral value depends on the time ellapsed along the N iterations; therefore, the algorithm depends on the simulator time steps. I don't like this.

I reimplemented it as an integration between t = 0 (when the controller starts to operate) and T. According to wikipedia, "The contribution from the integral term is proportional to both the magnitude of the error and the duration of the error." In my case, I don't think I need the duration of the error. Thus, so far I'm using Ki = 0.

It's possible, that your Kd should be negative, because the faster you are proceeding, the lesser the power should be. Noticed any similarities to your code? :) You assumed there that Kd = -1. Above you have a more formal controller's description. Even more formally, the Kd should stay positive when you apply an assumption, that positive velocity means moving away from target and negative means towards target.

Right, since error = target - current, then d(error)/dt = -d(current)/dt. I'm using the following, with Kd > 0:

Code:
double needMod =      Kp * error + 
                      Kd * (-m_currentAngularVelocity) +
                      Ki * integrationValue;

Next thing I'll try to use your MultiVessel template.
 

Enjo

Mostly harmless
Addon Developer
Tutorial Publisher
Donator
Joined
Nov 25, 2007
Messages
1,665
Reaction score
13
Points
38
Location
Germany
Website
www.enderspace.de
Preferred Pronouns
Can't you smell my T levels?
I published my code also in SF under GPL :). I think it's easier to add your source code files to my own. In this case, if I upload the code including your files to SF, would it be ok? (copyright and licences would be of course kept).

Sure. It goes along the license(s) anyway

I'm not limiting the signal and it seems to work. In addition to Kp and Kd, the key seems to be a multiplier that converts the result of P+I+D into the actual vector passed as the vessel RCS levels. If this multiplier is too low, the operation is slow; if it is too large, it overshoots. This multiplier will be related to the maxAccel :).
For the best experience, you must ultimately somehow limit the vector's components passed to the RCS. If you don't do it, it will work I think, because Orbiter understands values <-1, 1>. If they're beyond that, it probably limits them up to 1 or -1. If it's a vector, it may limit it so that its length is = 1, but although it will work, you're very much limiting the possible output resolution of your controller. For example, if your controller is able to smoothly output values from 0 to 5, Orbiter would use only 1/5th of this range - <0, 1> out of <0, 5>, meaning that it won't matter for your system that you're a bit off the target and you need to accelerate more, but not to the max (for example output = 4.75). Your system will accelerate at full speed, giving you a risk of oscillations that you can't remove with K's. On the other hand, when you're quite close, the controller could say, that you need to accelerate "only" with 0.75 out of 5, while for <0, 1> range, 0.75 is still quite a lot. A good exercise would be to output the resulting RCS power by using oapi, to see what actually goes on.

I tried an integrator that keeps N old values and timesteps. However, I'm not happy with it. This integral value depends on the time ellapsed along the N iterations; therefore, the algorithm depends on the simulator time steps. I don't like this.

I reimplemented it as an integration between t = 0 (when the controller starts to operate) and T. According to wikipedia, "The contribution from the integral term is proportional to both the magnitude of the error and the duration of the error." In my case, I don't think I need the duration of the error. Thus, so far I'm using Ki = 0.

First off - you don't necessarily need integration, unless you want a quicker reaction.
Secondly, my integration also keeps track of the previous steps - notice how I accumulate calculations:
Code:
m_integ += (e + m_ePrev) * dt / 2.0
which means that I keep track of the trapezoids from all previous steps until I call Reset(). This way you don't need to remember any N previous steps. Just one and the sum of all previous!
Lastly, you can't really complain about the integration part being dependent on variable timestep without blaming all the other parts for the same :) After all, all you get is an error from the current frame, velocity from the current frame, and after some time step, before which you're completely blind, you get another error and velocity. It's all the same for all parts, but generally, yeah. Don't integrate until you get the P and D parts right, with some oscillation maybe. The I part will remove these oscillations up to some degree (AFAIR)

Right, since error = target - current, then d(error)/dt = -d(current)/dt. I'm using the following, with Kd > 0:
No. It's because if error is defined as target - current and if target > current, then error > 0, and if target < current, then error < 0. Atan2 will handle it.

[EDIT]
I don't know why you've done it with the error. Pure (atan2(v1, v2) / PI) would be enough, where obviously (to me at least), the atan2 is current error, and PI is the maximal error, which gives you an error in range <-1, 1>. It seems to me that you're overcomplicating it to some degree of unreadability :)
 
Last edited:

evilfer

New member
Joined
Dec 9, 2010
Messages
53
Reaction score
0
Points
0
For the best experience, you must ultimately somehow limit the vector's components passed to the RCS. If you don't do it, it will work I think, because Orbiter understands values <-1, 1>. If they're beyond that, it probably limits them up to 1 or -1. If it's a vector, it may limit it so that its length is = 1, but although it will work, you're very much limiting the possible output resolution of your controller. For example, if your controller is able to smoothly output values from 0 to 5, Orbiter would use only 1/5th of this range - <0, 1> out of <0, 5>, meaning that it won't matter for your system that a bit off the target and you need to accelerate more, but not to the max (for example output = 4.75). Your system will accelerate at full speed, giving you a risk of oscillations that you can't remove with K's. On the other hand, when you're quite close, the controller could say, that you need to accelerate "only" with 0.75 out of 5, while for <0, 1> range, 0.75 is still quite a lot. A good exercise would be to output the resulting RCS power by using oapi, to see what actually goes on.

You're right!! Orbiter is limiting the RCS vector components to [-1,1]. I think it limits each component independently. This means that the limited vector does not keep the same direction. I've been checking the vector length, and it exceedes 1 sometimes. The effect is that the ship rotation does a curve (due to the change in direction in the limited vector), and actually the ship rolls a little (which is unnecessary).

Besides, the limitation in acceleration is very relevant, and provokes the system to overshoot :S.


First off - you don't necessarily need integration, unless you want a quicker reaction.
Secondly, my integration also keeps track of the previous steps - notice how I accumulate calculations:
Code:
m_integ += (e + m_ePrev) * dt / 2.0
which means that I keep track of the trapezoids from all previous steps until I call Reset(). This way you don't need to remember any N previous steps. Just one and the sum of all previous!
Lastly, you can't really complain about the integration part being dependent on variable timestep without blaming all the other parts for the same :) After all, all you get is an error from the current frame, velocity from the current frame, and after some time step, before which you're completely blind, you get another error and velocity. It's all the same for all parts, but generally, yeah. Don't integrate until you get the P and D parts right, with some oscillation maybe. The I part will remove these oscillations up to some degree (AFAIR)

Sorry, I missread your code. My complaint about timestep dependence does not apply to your integration, only to one that considers only a subset of previous values, but this does not make sense, it was my mistake.


No. It's because if error is defined as target - current and if target > current, then error > 0, and if target < current, then error < 0.

I don't know why you've done it with the error. Pure (atan2(v1, v2) / PI) would be enough, where obviously (to me at least), the atan2 is current error, and PI is the maximal error, which gives you an error in range <-1, 1>. It looks to me that you're overcomplicating it to some degree of unreadability :)

:S. That was because of the D component of P+I+D. D is the angular velocity, which it's a vector.

I've changed it again. error is scalar, d(error)/dt is calculated by comparing it with the previous error so it's also scalar. This is the current implementation:

Code:
	double pid = m_error * m_kP +
		m_integral * m_kI + 
		errorDt * m_kD;

	m_desiredVelocity = cross * pid / sin_angle;

	if (length(m_desiredVelocity) > m_maxAV) {
		m_desiredVelocity *= m_maxAV / length(m_desiredVelocity);
	}

	m_thrust= (m_desiredVelocity - avel) * 100;

	double max = 0;
	if (fabs(m_thrust.x) > max) max = fabs(m_thrust.x);
	if (fabs(m_thrust.y) > max) max = fabs(m_thrust.y);
	if (fabs(m_thrust.z) > max) max = fabs(m_thrust.z);
	
	if (max > 1.0) {
		m_thrust/= max;
	}

	r_vessel->SetAttitudeRotLevel(m_thrust);

I'm going swimming a little, need a rest :D.
 

Enjo

Mostly harmless
Addon Developer
Tutorial Publisher
Donator
Joined
Nov 25, 2007
Messages
1,665
Reaction score
13
Points
38
Location
Germany
Website
www.enderspace.de
Preferred Pronouns
Can't you smell my T levels?
:S. That was because of the D component of P+I+D. D is the angular velocity, which it's a vector.
Some misunderstanding here. What I meant is that I have no idea why you've done this:
Code:
m_error *= m_currentangle / sin_angle;
Why / sin_angle? I'm still with m_error = atan2(v1, v2) / PI, which stands for - current / maximal, thus naturally limiting the error to <-1, 1>

I've changed it again. error is scalar, d(error)/dt is calculated by comparing it with the previous error so it's also scalar.
So errorDt = (m_error - m_errorPrevious) / dt ? If so, than it's a great idea, because the sign of errorDt works for you now :)

Referring to your code - I think that you're doing a correct post-processing, but even before that, you should make sure that the following line:

Code:
if (length(m_desiredVelocity) > m_maxAV)
evaluates to true only on rare occasions, effectively reserving it for situations like having a given error and a huge velocity away from target. Then naturally your RCS should work on full power. Otherwise - less than full power. It should be achieved only by manipulating Kp and Kd. After all, the error and error/s are already in range of <-1, 1>. Because of that, again I frown at / sin_angle for m_desiredVelocity. What is it for? It's something I might simply not understand, but as far as I can see, this just gives some unneeded unpredictability, and let me rephrase - you should be able to limit the pid variable by only manipulating the K's. That's predictable and guaranteed to be very smooth in operation, because this way you'd be using full resolution of your PID.

The pid * cross is ok, as long as cross is normalized. See my Vect3 code for this. OrbiterSDK.h may have such function already, if you insist on a less object oriented approach :)

I'm going swimming a little, need a rest :D.
Heh, I've just came back from a swimming pool :)
 
Last edited:

evilfer

New member
Joined
Dec 9, 2010
Messages
53
Reaction score
0
Points
0
Some misunderstanding here. What I meant is that I have no idea why you've done this:
Code:
m_error *= m_currentangle / sin_angle;
Why / sin_angle? I'm still with m_error = atan2(v1, v2) / PI, which stands for - current / maximal, thus naturally limiting the error to <-1, 1>

...

Because of that, again I frown at / sin_angle for m_desiredVelocity. What is it for? It's something I might simply not understand, but as far as I can see, this just gives some unneeded unpredictability, and let me rephrase - you should be able to limit the pid variable by only manipulating the K's. That's predictable and guaranteed to be very smooth in operation, because this way you'd be using full resolution of your PID.

The pid * cross is ok, as long as cross is normalized. See my Vect3 code for this. OrbiterSDK.h may have such function already, if you insist on a less object oriented approach :)

Readibility issues :S. sin_angle = length(cross). Therefore, cross / sin_angle is a normalized vector (pointing to the rotation axis).


Referring to your code - I think that you're doing a correct post-processing, but even before that, you should make sure that the following line:

Code:
if (length(m_desiredVelocity) > m_maxAV)
evaluates to true only on rare occasions, effectively reserving it for situations like having a given error and a huge velocity away from target. Then naturally your RCS should work on full power. Otherwise - less than full power. It should be achieved only by manipulating Kp and Kd. After all, the error and error/s are already in range of <-1, 1>.

I have to play with the values, and also with Ki, see if I can achieve that :). Wikipedia has some info on tuning PID controllers... For manual tuning, I can also include buttons in the MFD to modify Ki, Kp and Kd, an have it record the time it takes to reach a given attitude.

Anyway, the controller does its job now :D.


----

By the way, after calculating pid, this value is used to calculate a desiredAngularVelocity. This is compared with the current angularVelocity to get a thrust value. Is this system then a cascade PID controller? (the first being a PID for attitude->angularVelocity, and the second being a P controller for angularVelocity->thrust)?
 
Last edited:

Enjo

Mostly harmless
Addon Developer
Tutorial Publisher
Donator
Joined
Nov 25, 2007
Messages
1,665
Reaction score
13
Points
38
Location
Germany
Website
www.enderspace.de
Preferred Pronouns
Can't you smell my T levels?
Readibility issues :S. sin_angle = length(cross). Therefore, cross / sin_angle is a normalized vector (pointing to the rotation axis).
See? This is why I have Vect3 with norm() method. This way everybody knows that I'm calling a normalization only of some specific vector, and nothing else.

By the way, after calculating pid, this value is used to calculate a desiredAngularVelocity. This is compared with the current angularVelocity to get a thrust value. Is this system then a cascade PID controller? (the first being a PID for attitude->angularVelocity, and the second being a P controller for angularVelocity->thrust)?
Not sure. Maybe, but I'd use the pid signal as the thrust right away, of course after limiting it. Anyway, I don't want to interfere anymore. Our philosophies are a bit different, so I'd rather do my own stuff with this, just to prove that it's possible. What's your sf project page?
 

evilfer

New member
Joined
Dec 9, 2010
Messages
53
Reaction score
0
Points
0
See? This is why I have Vect3 with norm() method. This way everybody knows that I'm calling a normalization only of some specific vector, and nothing else.

Definitely I agree. Anyway since I'm loading orbitersdk apis, I think I'll use the functions for VECTOR3. I'll change the name of that variable though :).


Not sure. Maybe, but I'd use the pid signal as the thrust right away, of course after limiting it. Anyway, I don't want to interfere anymore. Our philosophies are a bit different, so I'd rather do my own stuff with this, just to prove that it's possible. What's your sf project page?

What I like of the current implementation is that I have control on both the first and second derivates of the parameter I want to set (attitude); since the system can affect only the second derivate, I feel comfortable this way. Anyway, if you improve the algorithm, please let me know. Anyway, I could get to the current algorithm only with your help!!! :)

My project page: http://sourceforge.net/projects/trajectories/
 

Enjo

Mostly harmless
Addon Developer
Tutorial Publisher
Donator
Joined
Nov 25, 2007
Messages
1,665
Reaction score
13
Points
38
Location
Germany
Website
www.enderspace.de
Preferred Pronouns
Can't you smell my T levels?
Definitely I agree. Anyway since I'm loading orbitersdk apis, I think I'll use the functions for VECTOR3. I'll change the name of that variable though :).
That's because your apps are purely Windows based and you can include OrbiterSDK because you always have windows.h on your system. I use wxWidgets which is multiplatform and also supports OpenGL and SDL. I build my apps mainly on Linux. Better yet, even some classes containing visual logic are completely independent of any GUI toolkit.
Have a look:
https://sourceforge.net/projects/enjomitchsorbit/files/1.4.5/wxApps/

I can't run your app. zlib1.dll is missing


What I like of the current implementation is that I have control on both the first and second derivates of the parameter I want to set (attitude); since the system can affect only the second derivate, I feel comfortable this way. Anyway, if you improve the algorithm, please let me know. Anyway, I could get to the current algorithm only with your help!!! :)
It was fun refreshing my memory. Firstly I would seperate as many small things as possible. PID may be a seperate class for example. I'll experiment some and show you the results. What you do with it is up to you.
 
Last edited:

evilfer

New member
Joined
Dec 9, 2010
Messages
53
Reaction score
0
Points
0
That's because your apps are purely Windows based and you can include OrbiterSDK because you always have windows.h on your system. I use wxWidgets which is multiplatform and also supports OpenGL and SDL. I build my apps mainly on Linux. Better yet, even some classes containing visual logic are completely independent of any GUI toolkit.
Have a look:
https://sourceforge.net/projects/enjomitchsorbit/files/1.4.5/wxApps/

I can't run your app. zlib1.dll is missing

Only the MFD is windows based, the simulator uses glut and sometime ago it did compile and run on Linux. Then I included the FTGL and cspice libraries and I stopped trying to compile it on linux :S.

For GUI toolkit, in my sim I decided to implement my own widget system on top of OpenGL, it was an interesting exercise.

About the library... I don't use any compression function so I have no idea why it asks for it. This is the first time for me to try to distribute an c++ application, java was easier.


It was fun refreshing my memory. Firstly I would seperate as many small things as possible. PID may be a seperate class for example. I'll experiment some and show you the results. What you do with it is up to you.

:). How about a template <class errorType, class outputType>. Please keep me informed. I think I'll move on to the next step: performing the maneuver main thrust.
 

Enjo

Mostly harmless
Addon Developer
Tutorial Publisher
Donator
Joined
Nov 25, 2007
Messages
1,665
Reaction score
13
Points
38
Location
Germany
Website
www.enderspace.de
Preferred Pronouns
Can't you smell my T levels?
Wow. I've replaced space autopilot with PID and it's so damn smooth! I'll get rid of all this fuzzy logic overdose :)

For GUI toolkit, in my sim I decided to implement my own widget system on top of OpenGL, it was an interesting exercise.
Yeah. I've seen it. A lot of work.

About the library... I don't use any compression function so I have no idea why it asks for it. This is the first time for me to try to distribute an c++ application, java was easier.
That's probably for texture support. For distributing windows apps, you'll need "Dependency walker" which shows the dlls used and missing.

:). How about a template <class errorType, class outputType>. Please keep me informed. I think I'll move on to the next step: performing the maneuver main thrust.
Nah, that would be too much. Perhaps I'll finish this today and I'll show you.
 
Top