Pidloop improvements

[ayzenquwe]

This PR contains several improvements, fixes and unit tests for PIDLoop. Related issues: #2919, #2920, #2921.

These changes should not affect existing kOS scripts and I tried my best to preserve current PIDLoop behavior (minus bugs).

Motivation

Several months ago I decided to start learning how PID controllers work. I've been using MATLAB to model simple things for KSP and that was a very fun experience: model something in MATLAB and make it work in KSP. After some time playing with MATLAB and KSP I noticed that PIDLoop in some cases works differently than PID controller in MATLAB, so I thought that it's a nice challenge for me - to understand why.

The first problem (#2919) I found was that integration in kOS PIDLoop skips the first error value. I fixed it in the first commit by using correct Forward Euler method [1] for integration. This change has almost zero impact on kOS PIDs, but still it was important for me to fix it before moving forward.

In the second commit I added simple unit tests to make sure I don't break anything during my next improvements and refactoring. I wasn't able to put unit tests before any changes due to the bug I had to fix in the first commit.

The second problem (#2921) is related to kOS PIDLoop anti-windup method. PIDLoop has a different and a bit simpler method than methods which can be found in MATLAB. This makes modeling anything complex in MATLAB near impossible. The third commit addresses this problem and adds additional anti-windup methods: clamping and back-calculation. [2] If anti-windup method is not set explicitly PIDLoop will work without any changes - using the currently implemented method.

ANTIWINDUPMODE (StringValue) is added to PIDLoop with these values:
NONE - output saturation (min/max limits) is applied, but anti-windup is not used;
CLAMPING - clamping algorithm (see [2]);
BACK-CALC - back-calculation algorithm (see [2]);
DEFAULT or any other value - currently implemented anti-windup method.

For the back-calculation algorithm a coefficient KBACKCALC (ScalarValue) is added. See [2] for details.

The diagram below shows the difference between these methods (telemetry for the diagram was gathered in KSP - vertical speed stabilization to 10, then to 0, using four J-20 "Juno" engines).

Besides new anti-windup methods I have refactored the PIDLoop Update method, removed obsolete variables and tried to make it more clean.

The third problem (#2920) is a bug in the derivative part: currently, PIDLoop calculates it based on input, but an error should be used instead. The fourth commit contains a fix for that.

If you think these fixes/changes may be helpful for kOS, I can work on documentation and add it as a separate commit to the PR (although I'm not a native speaker, so someone will need to check it :) ). If needed I can break the PR down into several PRs.

References

[1] https://www.sciencedirect.com/topics/engineering/forward-euler
[2] https://www.matec-conferences.org/articles/matecconf/pdf/2017/26/matecconf_imane2017_05013.pdf

[Dunbaratu]

Can you explain why is it supposed to be better to time-shift the I term one iteration like this so it's using the previous error value instead of the current one from this iteration? The issue #2919 only has one sentence in the description which doesn't explain why it's desired to delay it like this.

[ayzenquwe]

I have added some explanations to PR comments, but I shouldn't have changed the integration method, so I'll revert this.

[Dunbaratu]

If you do this, don't forget that it's not enough to just add documentation of the new changes, but also change the existing documentation that shows pseudocode to the reader showing what kOS PIDloop does.

Specifically I mean this: https://ksp-kos.github.io/KOS_DOC/structures/misc/pidloop.html?highlight=pidloop#pidloop-update-derivation

[Dunbaratu]

If this value is going to be checked repeatedly during Update it's probably better to store it internally as an enum type and check the enum here instead of doing string manipulation on the fly again and again each Update(). (i.e. when a script uses set to change the value of the suffix, check it then to see if the string is one of the known values and if it is then set the enum, else throw an error.) That way the string work only happens when the value is set, rather than 50 times a second.

[ayzenquwe]

You're right, I'll fix this. Thanks!

[Dunbaratu]

I have to admit the math here is a bit beyond me. (the reference links provided seemed to skirt around the issue without really describing it head-on). But I'm willing to accept the PR if people who understand it better can chime in and check it.

[ayzenquwe]

If you do this, don't forget that it's not enough to just add documentation of the new changes, but also change the existing documentation that shows pseudocode to the reader showing what kOS PIDloop does.

Specifically I mean this: https://ksp-kos.github.io/KOS_DOC/structures/misc/pidloop.html?highlight=pidloop#pidloop-update-derivation

Thanks for your quick response! I'll definitely update the pseudocode. I just wanted to get some feedback before I start doing this.

I'll answer your other questions later. Thanks!

[nuggreat]

I see 2 possible issues.

First the change in how the D term is calculated from using the input to using the error could break existing scripts that are relying on this functionality. See dunbaratu's aircraft autopilot mentioned in issue #2693.

Second the changes to how the I term is calculated don't quite look right to me for a reason. Because it looks like the I term will now lag one update call behind the P and D terms. As a result of the lag introduced into the I term this means that the applied dt value will now no longer be from the same update that the I term was calculated in and as a PIDupdate won't always be called with the same dt this could cause major issues. See example code below for what I think might happen assuming I am reading the code right.

LOCAL demoPID IS PIDLOOP(0,1,0,-100,100).

PRINT demoPID:UPDATE(0.0, 0).// prints 0
PRINT demoPID:UPDATE(1.0, 1).// prints 0
PRINT demoPID:UPDATE(1.1, 1).// prints -0.1

[mgalyean]

Neither skipping the first I term as the pre-PR code did nor time-shifting the I terms as this PR does makes sense to this math dilletante. Does Forward Euler really imply time shifting like that for the PID implementation? I've done some quick searches and reading and I don't see where the need for the time shift comes in and am open to education on it. Given that the I term has a kind of built in "lag" (the integration aspect), it just doesn't seem right to delay that signal another tick. Especially given odd behavior that can occur when dealing with time-skewed simulation signals (though I think the I term would be a bit more immune to this skew sensitivity than P or especially the D I suppose)

[ayzenquwe]

First of all, thanks for your feedback, I really appreciate it.

Regarding the I term calculation, both pre-PR and PR methods are approximations, they are very similar and also they aren't really correct compared to real integration. Pre-PR method is called Backward Euler Method, PR method is Forward Euler Method. They're both related to finding numerical approximations to the solutions of ordinary differential equations and also they can be used for numerical integration. They're used in MATLAB for discrete time integrators and as far as I know there are many other, more complex methods, that are used for continuous time simulations, but it's really beyond my knowledge. MATLAB uses Forward Euler method as a default method for discrete time integrators and for discrete PID controllers (their documentation states that it's the best choice for small time steps, but they don't explain why).

MATLAB discrete-time integrator documentation has quite simple explanations of how these methods are calculated. Unfortunately all other articles I found are quite complex with lots of math.

To illustrate how these methods approximate integration here is a very simplified example:
Backward Euler Method

Forward Euler Method

The green line is the error signal we want to integrate, the red line is approximated error signal and gray blocks are the result of integration. As you can see both methods are very similar and they're not precise: integration from 1 to 3 is 4, but Backward Euler gives 5, Forward Euler - 3. There are some subtle differences between the methods, but they are related to solving differential equations. Also the Forward Euler method can be unstable when discretizing a system that is stable in continuous time, but as far as I understand it's not applicable to our cases.

My initial intention was only to fix the skipped first value problem, but eventually I changed the integration method to Forward Euler. I have to admit that wasn't a very good idea, so I'll revert it back to Backward Euler with fixed first value. Also this will allow me to put some initial unit tests (where the first value is 0) in the first commit.

As for D term, I completely missed the part of documentation that describes D term differences from classic PID. Although this looks strange to me I understand that it works as intended.

So the only difference from pre-PR PID implementation will be the first value problem fix (I'll revert integration method and D term calculation). It's a very small fix, but still can affect some scripts in rare cases... Maybe it's better to introduce a "classicPid" mode and put it there along with classic D term calculation? What do you think?

Also, should I create a new PR for my new changes or just modify this one?

[KrazyKerbalnaut]

There is another approximation of integrals. It is called the trapezoidal sum and is relatively inexpensive for computers since we are already doing most of the work already. There are 3 additional operations to perform. This gives us a better approximation than the Riemann method because it takes slope into account. While not perfect, it gets far closer than the left or right Riemann sums.

The function is illustrated below:

This can be programmed as follows in kerboscript and can be easily translated to c#:

// initialize required variables that are static throughout the main loop
SET t0 TO time:seconds.
SET errorLast TO 0.0.
SET Iarea TO 0.0.

...
// some delay before we ever call it for the first time
// could also have 0 delta checking in the script like the original PID loop
...

// integral trapezoidal summation
FUNCTION intTrapSum {
    // get the current error state
    PARAMETER error.

    // time functions
    LOCAL t1 TO time:seconds.
    LOCAL deltaT TO t1 - t0.

    // update time for the future iteration
    SET t0 to t1.

    // perform the trapezoidal summation of the last error to the current error for the I term (Iarea).
    SET Iarea TO Iarea + deltaT * (( error + errorLast) / 2).

    // update the last error for the future iteration
    SET errorLast TO error.

   RETURN Iarea.
}

Also, I had closed the bug report for item #2938. Before any limiting function can be done to a PID loop, we must always compute the terms and then clamp to the requested min/max. We cannot perform clamping on the integral term based on that. If we do, we have a time lagging system that grows in error relative to the input error, until the system is unclamped as a whole.

Instead, to perform integral anti windup AND max PID values, we need 2 sets of limiting terms. We need the minimum/maximum allowable iterm values AND we need the minimum/maximum allowable PID output values. Then we figure out if we need to clamp the output AFTER the P, I, and D terms are computed. This way, when the error reverses, the I term doesn't depend on the P (not really the P term since it will also immediately reverse) and the D term (yes, the calculates the rate of change and applies resistance to that. So if we suddenly go negative, it will resist the change and actually increase the overall output, keeping it clamped, and the Iterm is unable to start its own reversal until the system is unclamped. At least in the current system).

If we perform anti-windup before we clamp the system as a whole, even if the derivative pulls us even higher, the I term will begin descending, actually helping us get out of the system clamp, even if I term was already at its max windup.

I know this will break some code. Even if the kOS library is able to execute the function 3 times faster, I still get better performance, and the expected performance when I use my own code in pull request #2938.

Cheers!

[KrazyKerbalnaut]

There is another approximation of integrals. It is called the trapezoidal sum and is relatively inexpensive for computers since we are already doing most of the work already. There are 3 additional opperations to perform. The function is illustrated below:
...

I made a boo boo and almost did derivation instead of integration on my last reply. I've updated the formulation

[Dunbaratu]

@ayzenquwe What is the state of this PR? Is it ready to merge? I see the most recent discussion said you were going to make another change (going from forward-euler back to backward-eurler) but I see no new commits after that comment, which is why I never merged this yet. I thought you weren't done yet.

As for this question:

Also, should I create a new PR for my new changes or just modify this one?

I would say modify this one (push a new commit to this branch). It's easier for me to merge in that way.

[ayzenquwe]

Hi @Dunbaratu. I was waiting for answers. :)

Can you please answer to this one?

So the only difference from pre-PR PID implementation will be the first value problem fix (I'll revert integration method and D term calculation). It's a very small fix, but still can affect some scripts in rare cases... Maybe it's better to introduce a "classicPid" mode and put it there along with classic D term calculation? What do you think?

[Dunbaratu]

Maybe it's better to introduce a "classicPid" mode and put it there along with classic D term calculation? What do you think?

I think the hassle of maintaining two PID codebases would be more hassle than dealing with the occasional person who's PID controller was dependent on the old behavior. I think we don't need a "classicPID". I'll just make sure the changelog is very clear.

That is, assuming the following phrase happens:

(I'll revert integration method and D term calculation)

If that's not done then there would be enough difference to need a classic version, but since you said you were going to do that, I don't think the remaining difference is big enough to matter.

[ayzenquwe]

Understood. I will need some time to implement it, since I won't be near my computer for 4 weeks. Will do it as soon as I return.