Smarter PID, PID autotune etc

It's all about the code!
Post Reply
User avatar
Site Admin
Posts: 9852
Joined: Wed Aug 28, 2013 1:28 am
Location: Jersey City
Soldering skill: yes
Coding skill?: yes

Smarter PID, PID autotune etc

Post by AndreyB » Thu Oct 24, 2019 2:04 am

@andreika is leading our PID modernization. We are trying to get a solution for PID autotune, so far our first use-case is alternator PID on Miata NB2.

For industrial PID proposal, we are thinking derivativeFilterLoss and antiwindupFreq.

Some research is done in a separate repository at

Keywords: ... _algorithm

Code: Select all

C:\stuff\autopid.rusefi\examples\2003_mazda_miata_alternator>..\..\pid_from_msl.exe "miata_alt_log.msl" 66.737 71.662 54 10 13.8 
PID_FROM_MSL - find PID controller coefficients based on a measured step response in a rusEFI log file.

Version 0.3 (c) andreika, 2019

Reading file miata_alt_log.msl...

Project Settings: FLT_EVAL_METHOD=0 sizeof(float_t)=8 sizeof(double_t)=8

Measuring Settings: targetValue=13.8 minValue=20 maxValue=30 stepPoint=178 maxPoint=449 numPoints=462 timeScale=91.1675

1) Trying method AUTO1 on "Regulator" PID model:

* modelBias=7.88118 avgMin=13.4117 avgMax=16.1769

* Params0: K=0.276525 T1=0.440947 T2=0 L=0.285518

* Solving...

[0] (0.276313,0.598212,0.234637,0) l=0.001 merit 34.0112->27.4315, dm=-6.57969

[1] (0.276887,0.664806,0.191864,0) l=0.0001 merit 27.4315->26.6702, dm=-0.761283

[2] (0.276883,0.684747,0.171325,0) l=1e-05 merit 26.6702->26.4857, dm=-0.18457

[3] (0.276877,0.698491,0.156908,0) l=1e-06 merit 26.4857->26.4043, dm=-0.0813515

[4] (0.276871,0.705527,0.149414,0) l=1e-07 merit 26.4043->26.3831, dm=-0.0212164

[96] (0.276862,0.708473,0.145976,0) l=1e-11 merit 26.3801->26.3801, dm=2.4869e-14

[97] (0.276862,0.708473,0.145976,0) l=1e-10 merit 26.3801->26.3801, dm=3.44613e-13

[98] (0.276862,0.708473,0.145976,0) l=1e-11 merit 26.3801->26.3801, dm=-3.2685e-13

[99] (0.276862,0.708473,0.145976,0) l=1e-12 merit 26.3801->26.3801, dm=-4.51195e-13

* The solver finished in 100 iterations! (Merit: 34.0112 -> 26.3801)

* Solving for better coefs:

* The solver aborted after 0 iterations! (Merit: 2783.17 -> 2783.17)

2) Trying method AUTO2 on "Regulator" PID model:

* modelBias=7.88118 avgMin=13.4117 avgMax=16.1769

* Params0: K=0.276525 T1=5.06304 T2=0.246866 L=-0.968021

* Solving...

[0] (0.2776,1.99779,0.1,0.01) l=0.001 merit 704.978->191.579, dm=-513.399

[1] (0.275877,0.53173,0.1,0.105225) l=0.0001 merit 191.579->27.8961, dm=-163.683

[2] (0.275877,0.53173,0.1,0.105225) l=0.001 merit 27.8961->48.379, dm=20.483

[3] (0.276545,0.513536,0.1,0.0688601) l=0.0001 merit 48.379->33.1391, dm=-15.2399

[4] (0.276545,0.513536,0.1,0.0688601) l=0.001 merit 33.1391->51.093, dm=17.9539

[5] (0.27759,0.884054,0.1,0.01) l=0.0001 merit 51.093->29.9555, dm=-21.1375

[6] (0.27759,0.884054,0.1,0.01) l=0.001 merit 29.9555->319.579, dm=289.624
[97] (0.276477,0.663123,0.1,0.0921907) l=0.01 merit 26.7273->26.0569, dm=-0.670377

[98] (0.276477,0.663123,0.1,0.0921907) l=0.1 merit 26.0569->39.1889, dm=13.1319

[99] (0.276638,0.643867,0.1,0.0330709) l=0.01 merit 39.1889->27.3687, dm=-11.8202

* The solver finished in 100 iterations! (Merit: 704.978 -> 27.3687)

* Solving for better coefs:

[0] (19.1186,27.5346,0.588119,0.0194127) l=1 merit 266.429->252.835, dm=-13.5936

[1] (16.5605,19.5184,0.377857,0.0194127) l=0.1 merit 252.835->185.372, dm=-67.4635

[2] (13.5241,12.8894,0.0185283,0.0194127) l=0.01 merit 185.372->134.617, dm=-50.7552

[3] (12.9896,19.017,0.567793,0.0194127) l=0.001 merit 134.617->134.409, dm=-0.207611

[4] (10.8851,12.9655,-0.356962,0.0194127) l=0.0001 merit 134.409->117.901, dm=-16.5077

[5] (10.8851,12.9655,-0.356962,0.0194127) l=0.001 merit 117.901->119.807, dm=1.90552

[6] (10.8851,12.9655,-0.356962,0.0194127) l=0.01 merit 117.901->118.346, dm=0.445136

[7] (10.8851,12.9655,-0.356962,0.0194127) l=0.1 merit 117.901->120.412, dm=2.51062

[8] (11.4456,16.8133,0.122291,0.0194127) l=0.01 merit 117.901->117.269, dm=-0.632365

[9] (11.4456,16.8133,0.122291,0.0194127) l=0.1 merit 117.269->138.488, dm=21.2192

[10] (11.4456,16.8133,0.122291,0.0194127) l=1 merit 117.269->126.189, dm=8.92042

[11] (11.4422,17.299,0.358441,0.0194127) l=0.1 merit 117.269->117.063, dm=-0.206191

[12] (11.4422,17.299,0.358441,0.0194127) l=1 merit 117.063->119.914, dm=2.85096

[13] (11.4422,17.299,0.358441,0.0194127) l=10 merit 117.063->121.764, dm=4.70095

[14] (11.4422,17.299,0.358441,0.0194127) l=100 merit 117.063->117.584, dm=0.521006

[15] (11.4422,17.299,0.358441,0.0194127) l=1000 merit 117.063->117.089, dm=0.0257627

[16] (11.4422,17.299,0.358441,0.0194127) l=10000 merit 117.063->117.066, dm=0.00360832

[17] (11.4422,17.299,0.358441,0.0194127) l=100000 merit 117.063->117.063, dm=0.000368467

[18] (11.4422,17.299,0.358441,0.0194127) l=1e+06 merit 117.063->117.063, dm=3.69496e-05

[19] (11.4422,17.299,0.358441,0.0194127) l=1e+07 merit 117.063->117.063, dm=3.69601e-06

[20] (11.4422,17.299,0.358441,0.0194127) l=1e+08 merit 117.063->117.063, dm=3.69611e-07

[21] (11.4422,17.299,0.358441,0.0194127) l=1e+09 merit 117.063->117.063, dm=3.69611e-08

* The solver finished in 22 iterations! (Merit: 266.429 -> 117.063)


Model-1 Params: K=0.276862 T1=0.708473 T2=0 L=0.145976

  PID0: P=25.78098074 I=32.99071470 D=1.70595759 offset=21.00000000 period=11.00000000ms

  Metric0 result: 2783.17 ITAE=2783.17 ISE=26945.9 Overshoot=0%

  PID1: P=25.78098074 I=32.99071470 D=1.70595759 offset=21.00000000 period=11.00000000ms

  Metric1 result: 2783.17 ITAE=2783.17 ISE=26945.9 Overshoot=0%

Model-2 Params: K=0.276638 T1=0.643867 T2=0.0330709 L=0.1

  PID0: P=19.57612646 I=28.91866148 D=0.61577173 offset=21.00000000 period=11.00000000ms

  Metric0 result: 266.429 ITAE=138.854 ISE=180.596 Overshoot=11.2949%

  PID1: P=11.44222099 I=17.29904723 D=0.35844114 offset=21.00000000 period=11.00000000ms

  Metric1 result: 117.063 ITAE=109.151 ISE=155.763 Overshoot=2.81278%

The best PID: P=11.44222099 I=17.29904723 D=0.35844114 offset=21.0 period=11.0ms
very limited telepathic abilities - please post logs & tunes where appropriate -
my skype is arro239

User avatar
Posts: 421
Joined: Mon Feb 13, 2017 2:35 pm
Location: Kiev

Re: Smarter PID, PID autotune etc

Post by andreika » Sat Oct 26, 2019 9:52 am

Let me try to explain the main idea behind this autotune method.

** Part I **

Let's assume that we want to find some working PID coefficients for a real device ("plant"). I'll use @Russian's Miata alternator as an example.
The plant has an input (the PWM signal which controls the alternator) and an output (the battery voltage which we can measure). The PID controller also has an input and an output. And finally we have our setpoint value (the desirable target voltage). Both plant and controller are connected:
pid_01.png (2.28 KiB) Viewed 279 times
If we want to determine the coefficients of the controller, we need to understand the plant's "properties" first. And the key point is that we can measure the plant's "response" to some impulse, and then determine the properties of the plant from there:
pid_02.JPG (38.87 KiB) Viewed 279 times
We change the plant's input (pwm signal is white) and measure the output (the battery voltage is red). We see that the plant doesn't react instantly, there's some inertia or delayed processing taking place. And that's exactly what we want to know about the plant!
The plant itself can be rather complex device (especially if we want to control the combustion engine), and we don't want to be concerned with its complexity. We want to find a rather simple mathematical model which responds to the signal nearly the same way as our real-world plant.
We can analyze our measured data and treat is as a transfer function:
pid_03.png (15.43 KiB) Viewed 279 times
And find an approximation of the response using one of the well-known smooth analytical curves:
pid_04.png (16.5 KiB) Viewed 279 times
The Automatic Control Theory says that most of real-world plants can be approximated with the "first order", "second order" or high-order transfer functions, which differ by curvature:
pid_05.png (19.98 KiB) Viewed 279 times
Also, the transfer function may have a time delay - this helps to better approximate some really slow-starting processes.
In our research we've confined ourselves to the 1st and 2nd order functions plus time delay - FOPTD and SOPTD.
The transfer function has the following parameters:
1) The gain - K. It affects the amplitude of the plant's output in relation to the input.
2) The time constant - T. First order function has 1 parameter, and the 2nd order has two: T1 and T2. It can be shown that the SOPDT degenerates to the FOPDT if one of the time constants becomes very high and another is very low:
pid_06.jpg (39.61 KiB) Viewed 279 times
3) The delay time - L.

So our task is to:
1) measure the step response data (it's called "manual bump test");
2) choose FOPDT or SOPDT function curve;
3) find it's 3 or 4 parameters (K, T(T1/T2), L) for the measured data.

We need some sort of regression analysis ("least squares problem"). There're several mathematical methods that would help.
We've used an optimization method called Levenberg-Marquardt (LM) to fit the curve. It combines the Gauss-Newton and "gradient descent" methods which takes less iterations to converge. And uses a linearization approximation to compute the partial derivatives for the Hessian matrix.

This method requires that we have some "reasoned" initial parameters and tries to make them better.
We use several different methods for that:
1) "2 points" for FOPDT and "3 points" for SOPDT.
See: Rangaiah G.P., Krishnaswamy P.R. Estimating Second-Order plus Dead Time Model Parameters, 1994.
2) Harriott's method for SOPDT. See: Harriott P. Process control (1964). McGraw-Hill. USA.
Also other methods are described in the book: Visioli A. Practical PID Control, pp. 169-187.

We also need some initial signal pre-processing in order to get the basic knowledge:
- when the input step starts;
- when the output saturation ends;
- what's the output min. and max. levels, the signal noise ratio etc.

This is another example of LM algorithm fitting the FOPDT function to the noisy engine RPM signal (thanks @darxfame):
pid_07.png (55.19 KiB) Viewed 279 times

*** to be continued ***

User avatar
Posts: 421
Joined: Mon Feb 13, 2017 2:35 pm
Location: Kiev

Re: Smarter PID, PID autotune etc

Post by andreika » Sat Oct 26, 2019 9:53 am

** Part II **

Now we have the analytical transfer function defined and we want to find the best PID controller coefficients for it.
But first we need to define "The best".
There are two types of controllers:
1) "Servo", the set-point is ever-changing and the controller changes the output as fast and precise as the input;
The example: ETB controller.
2) "Regulator", where the set-point is quasi-fixed, but the output is disturbed by some external influences.
The examples: Alternator controller, Idle IAC or Timing controllers.

Also, there's another trade-off: response speed vs overall controller accuracy (which can be measured differently: overshoot factor, oscillation and instability factors etc.).

So the best way to choose between different PID coefficients is to simulate each and every controller in a virtual "sandbox" on some pseudo-real data and measure it's speed/accuracy using our own "merit function".

In order to achieve this, we need:
1) Calculate the initial "estimated" coefficients for our transfer function parameters;
2) Choose a PID controller model and input data depending on it's kind (servo or regulator);
3) Prepare a merit function to compare the controllers;
4) Do a regression analysis ("least squares problem") using our favourite LM algorithm;

* * *

The initial coefficients are rather important because the residual function is not smooth anymore as it was in transfer response fitting. The solver can easily "stuck" in some local extremum point. There're many known formulas for PID estimations. Each author has tried to reach his own goals by tuning the coefficients and solving trade-offs. One of the most complete references with all methods gathered is: O'dwyer A. Handbook of pi and pid controller tuning rules. 3rd Edition, 2009.
pid_10.png (94.96 KiB) Viewed 279 times
We've implemented several formulas:
1) Chien-Hrones-Reswick PID implementation for FOPDT: generic model, set-point regulation, disturbance rejection;
2) "IMC-PID" Rivera-Morari-Zafiriou for FOPDT;
3) "IMC-Chien" (aka Rivera/Smith) for SOPDT;
4) Van der Grinten Model (1963) for SOPDT (disturbance);
5) Haalman-Pemberton model for SOPDT (for overdamped response and significant delay);
6) IMC-Maclaurin model for SOPDT;

The default are Chien-Hrones-Reswick models (CHR). See: Chien, I.L., IMC-PID controller design - An extension.:
pid_11.jpg (100.32 KiB) Viewed 279 times
* * *

We use two types of "sandbox" input data for PID controller simulation:
1) For "servo" type, we imitate a step response:
pid_12.png (18.69 KiB) Viewed 279 times
2) For "regulator" type, we imitate output disturbance with a step and periodical noise:
pid_13.jpg (40.26 KiB) Viewed 279 times
Also, we've implemented and tested two controller models:
1) Parallel controller (identical to the PID implementation in rusEFI firmware);
2) Industrial controller, with derivative filtering and integrator anti-windup;

* * *

The merit function for the LM solver may use various criteria:
1) Integral time-weighted absolute error (ITAE);
2) Integral square error (ISE);
3) Max. overshoot factor;

We've combined both ITAE and squared Max.Overshoot (in %) as a balanced factor.

* * *

** The conclusion **

We've made an utility (PID_FROM_MSL) to calculate "the theoretical best" PID coefficients from user-supplied log-file of the plant's step response.
This research has shown that it's possible to use mathematical numerical methods based on abstract models for PID controller tuning.
At the same time it's not possible to take into account all real-world factors, so the obtained results can be used as a starting point for futher manual tuning.

Thanks for reading! :)

Post Reply