# Moon lander project: thrust vectoring with a PID controller

1. Moon lander project: thrust vectoring with a PID controller
2. Moon lander project: PID-based vs TGO-based guidance
3. Moon lander project: preliminary study of the guidance software

The goal of my Moon lander guidance software series is to learn and gain some experience with an end-to-end aerospace software project, especially on the topic of GNC. I have already written a few articles sharing my findings: for instance, this one about guidance.

Today I want to talk about control. I will use the engine gimbal controller of my Moon lander1 project to explain how to design, tune and validate a PID controller. I will first explain how the control algorithm is designed as a two-steps process with first a translation and then a rotation controller. I will then detail how I designed the rotation controller, which is responsible for the engine gimbal. Finally, I will show the implementation and explain how the performances were validated.

## Translation vs rotation control

The control of the spacecraft is implemented as a two-steps process: first, a translation controller, and then, a rotation controller. This naming is probably not industry standard at all, so let’s explain each one in details.

The translation controller takes as input the guidance acceleration commands x and y, and outputs a thrust and an angle. The idea is to determine the ideal angle the rocket should be pointed at, as well as its thrust, to fulfill the commanded guidance. This assumes that the rocket points magically and instantaneously to the commanded angle: it does not takes into account engine gimbal. This translation control is easy: there are no PID controller, no overshoot, no steady state error. It is only some basic trigonometry. The code, implemented here is just under 30 lines. The controlled thrust is a simple function of the commanded acceleration and the spacecraft’s mass: `ctr_thrust = gui_acc * sc_mass`. The controlled angle is the output of the `arctan` function: `ctr_angle = atan2(gui_acc.y, gui_acc.x)`.

With a commanded rocket thrust and rocket angle, the second step of the control can take place. It determines the actual engine gimbal that will bring the rocket to the previously determined rocket angle. This rotational control is a bit more complicated: first of all, we have a commanded rocket angle, and (obviously) a current rocket angle. From these, we can derive an angle error, which we want to correct (bring to zero). We will need a controller, as there are no analytically perfect solution. We also have many other sources of complexity, such as the inertia of the spacecraft, or the engine gimbal limits (angular position and velocity).

This article focuses on this second step of the control. We will see together how to design, implement and tune, and finally validate the controller.

## Controller design

Unlike the example used in this article (control), the guidance is time-bounded and needs to take into account both position and velocity, which means it needs something smarter than a PID. In our case, for control, the engine gimbal is not time-bounded (it must follow a moving setpoint the best it can), and there is only one variable to control (the gimbal angle). Therefore, a PID fits very well this problem.

The PID’s input (setpoint) will be the commanded thrust angle. The PID’s process variable (the measured value) will be the current spacecraft angle. This seems to be a standard, easy PID controller: compute the error (process variable minus setpoint), and compute the PID’s output (commanded value).

This output control value is not directly the engine gimbal angle, but merely a spacecraft-related angular control. Indeed, there are a two constraints to keep in mind which require this generalization. First, the engine gimbal has a maximum angular position and angular velocity2. This is rather easy to manage and will be handled with saturation of the control output. Second, the mass (and therefore the moment of inertia) of the spacecraft changes. This needs to be taken into account somewhere in the control, to avoid having a PID performing well for only a particular mass and badly for other masses.

## Controller implementation and tunning

Now that the controller is designed, let’s implement it in Rust.

After the implementation, we will need to tune it to determine the gains `Kp` and `Kd`. The goal will be to find a couple (`Kp`, `Kd`) which gives a fast response with minimal overshoot.

Implementation

The PID’s implementation is rather simple, so I will not go too much into the details. Just two details. First, the I term is kept at zero, because in space (no atmosphere) there can not be steady state disturbances on the spacecraft. Second, instead of using the derivative of the error, `derr`, for `Kd` , I directly use the spacecraft’s angular velocity `sc_ang_vel` (which is the derivative of `sc_ang_pos`). This makes sure that `Kd` applies only on the spacecraft’s behavior, without the disturbances from the commanded angular position `ctr_ang_pos`.

The code is in Rust, with type hints to make the units more clear. It also uses the `uom` crate, to avoid mixing units.

 ``````1 2 3 4 5 6 7 `````` `````` // ang acc PID let err: Angle = ctr_ang_pos - sc_ang_pos; let derr: AngularVelocity = sc_ang_vel; let control: Angle = (kp*err + kd*derr).into(); let ctr_ang_acc: AngularAcceleration = (control * control_transfer_function).into(); ``````

The `control_transfer_function` is a simple factor to scale the PID’s output into something which makes more sense and can be used to actually control the spacecraft. In my case, it is a factor of `1/sec**2` (in fact, only the unit is important), to convert an `Angle` into an `AngularAcceleration`.

The angular acceleration can then be converted into a gimbal angle, taking into account the dynamics of the spacecraft (mass, thrust, etc). The last step is to check the engine physical limits (angular velocity and angle).

 `````` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 `````` `````` // compute torque for correction // Inertia: 1/2*m*r**2 = kg.m**2 // Torque: N*m = kg*m**2 * rad/sec**2 let sc_moment_of_inertia = 0.5 * sc_mass * squared!(conf.sc_width/2.0); let ctr_torque: Torque = (ctr_ang_acc * sc_moment_of_inertia).into(); // compute engine gimbal // Torque = L*F*sin(gimbal) let sin_gimbal: Ratio = ctr_torque/(conf.sc_height/2.0*cur_thrust); let mut ctr_eng_gimbal: Angle = sin_gimbal.asin(); // apply engine gimbal limits (angular velocity and angle) // [...] ctr_eng_gimbal = saturate( ctr_eng_gimbal, -conf.ctr_eng_gimbal_pos_max, conf.ctr_eng_gimbal_pos_max, ); ``````

Tunning

To find the gains `Kp` and `Kd` that give the best performances, I simulate the above PID controller with a standard step response. In my case, the measure of the performance is the settling time `tau`, which denotes the time it takes for the system to reach and stay within 90-110% of the setpoint. The lower this value is, the better.

The 90-110% value is the amount of overshoot which I consider acceptable. I assume that the exact value does not matter too much: the best (`Kp`, `Kd`) couple for 90-110% of the setpoint will be the same for 95-105% or 80-120%. I’m mainly interested in the relation between `tau`, `Kp` and `Kd`. I also assume that a sufficiently small `tau` will not give any oscillations, hence no (or minimal) overshoot: each oscillation outside 90-110% of the setpoint would increase the value of `tau`, making it not sufficiently small. In any case, after finding potentially good gain values, I will run some landing simulations, asses how stable the spacecraft is, and adjust them accordingly.

Here are the step responses for two differently tunned PIDs; both with a setpoint of 2 degrees. The first one, in blue, uses some default and low values for `Kp` and `Kd` (2 and -0.5, respectively). The second one, in orange, uses optimally tunned values for a `tau` of 0.5 sec: 12.1 and -3.7. The most interesting subplot is the first one, showing the angular position of the engine gimbal. Vertical bars denote the settling time `tau` for both controllers.

This first PID has a time constant `tau` of 9.0 seconds. It shows a large overshoot (~50% at 2 sec) and many oscillations before finally settling. The second PID has a time constant of 0.5 second. As I assumed, it settles very quickly and shows almost no oscillations. Clearly, the first one is not optimized, while the second one is.

One thing to keep in mind, which I did not took into account in the step response here, is the maximum angular acceleration that the engine gimbal can achieve: the first controller commands up to ~4 deg/sec2, while the second controller commands between 20 and 25. The engine gimbal can achieve a maximum angular acceleration of somewhere between 3 and 9 deg/sec2. This depends on the engine thrust and spacecraft mass (and a few other parameters which are constant during the flight, such as the inertia or the maximum gimbal angle). To avoid an unpredictable controller, the controller should be tunned to avoid too strong commands that are outside of the engine’s physical limitations.

To obtain a small `tau`, the PID needs a large `Kp`, to have a strong and quick correction. But this comes at the expense of overshoot and oscillations. `Kd` can be decreased (below zero, hence increase its absolute value) to avoid them. `Kd` takes into account the angular velocity, and the larger it is, the more it reduces the correction. This makes the PID braking in anticipation when moving too fast toward the setpoint.

In short: a small `tau` (big `Kp`) gives a fast correction, with large gimbal commands, and sometimes some overshoot and oscillations (depending on `Kd`, up to a certain amount). This tunning is good for a short flight, such as a 150 m hop (30 sec). Big `tau` (small `Kp`) gives a slow, smooth correction, with small gimbal commands and small to no overshoot. This is suitable for long flights, such as a landing (700 sec).

From `Kp` and `Kd` to `tau`

I’m lazy, and I want to specify only one parameter, not two (`Kp` and `Kd`). This one parameter is obviously `tau`, from which we can estimate the two gains.

Therefore, I need to determine two functions `Kp = f(tau)` and `Kd = g(tau)`. I will test all (`Kp`, `Kd`) couples, find their corresponding `tau` with a step response, and finally compute a regression of `Kp` and `Kd` as a function of `tau`. I used `scipy.optimize.curve_fit` to compute the regression, which proved easier than anticipated. The regression uses a power function `f(tau) = k_scale * tau**k_exponent` to fit the data. A few other regression functions were tested (linear, exponential) but did not fit as well.

The result is shown below. The first plot shows `tau`, in function of `Kp` and `Kd`. The red line highlights the lowest `tau`, for a given `Kp`. The second plot shows `Kp` and `Kd` as a function of `tau`, with the calculated regression.

The result of the regression is:

``````kp_scale = 4.224639
kp_exponent = -1.524106
kd_scale = -2.203941
kd_exponent = -0.759514
``````

For example, for `tau=0.5`, `Kp` will be `4.224639*0.5**-1.524106 = 12.150` and `Kd` will be `-2.203941*0.5**-0.759514 = -3.731` (these are the values I used in the step response in the previous section).

## Conclusion

This was a fun little subproject. I touched many different things: control theory (PID), Geogebra (for the animated gimbal), Dia (for the block diagram), matplotlib (for the other graphs), programing (Rust, Python), maths (regression), etc. I also learned and practiced some GNC. Even if it is not state of the art, it is still interesting, and good enough for many applications, including my little toy project. All the source files (code and diagrams) can be found in this archive.

I improved my GNC software step by step, learning new things and writing blog articles along the way. First, I assumed that the engine gimbal and the rocket are rigid and point immediately and perfectly the setpoint. Then, I introduced the engine gimbal and the rocket inertia, which created the need for a controller. Later I will add some noise (sensor noise, engine gimbal misalignment, etc) to perfect the simulation.

1. Actually, the plan changed a bit: it is not limited to the Moon, or landing, anymore. It can totally work for, say, an ascent from Earth. ↩︎

2. For example, Apollo had a maximum gimbal deflection of 6 deg, and maximum gimbal rate of 0.2 deg / sec. ↩︎