PID tuning - modified Skogestad

As mentioned in the previous post, in some cases, particularly for servo positioning systems, it's a good idea to dive a bit deeper into the theories of tuning. One method I have found useful is the Skogestad method.

All the details can be found the linked paper. This is a modified form that I have found useful for hydro-turbines in combination with the trial and error method. The method consists of:

• A step in input, plotting the output and the input.
• Finding characteristic values from that plot.
• Using a single parameter called tau_c as tuning parameter.

The reason that this method is sometimes useful for hydro turbines, is that these systems can be very dynamic and oscillating. It can therefore be very difficult to identify the 15% and 30% overshoot from the trial and error method with any accuracy at all (although not entirely impossible). The Skogestad method (modified) will give usable values in these cases, at least for Ti.

The modifications to the method are:

• Disregard everything about time delay. Although one can often find something looking like a time delay, and it can be 2-3 seconds long, this is NOT a time delay in the normal sense, and should not be treated as one. Treating it as a time delay will in most cases produce unworkable PID parameters (instability or way too much overshoot).
• Use a tau_c of approximately 5 as a starting point (subject to tuning, see further down)
  

The disadvantage of the method, is it only works for PI. It does not work for PID unless approximating the system as a second order system (which a hydro power plant definitely is not, and a rather impossible thing to do in the field anyway). A good point of the method, is that Ti always comes out with a good value. It's really only tau_c (Kp in practice) that would eventually be adjusted to fit performance specs. A tau_c of 5 is found by me as a good starting point on idealized systems. It could however be way off on real systems.

Well behaved system

Looking first at a "well behaved" system. In hydro-power plant circumstances, this will mean a very simple system, an idealized plant.

The system is set in open loop, and a step in "kappa" (guide vanes) of 10% (0.1 pu) is done. This is seen as the green line, where we also can see the speed of the servo. The red line is the turbine speed in "per units", starting at 1.0. Here we can clearly see the speed going in the opposite direction during the first 2 seconds. This is fairly typical, but it's not really a time delay. The pink line shows the 63% rise in speed and the corresponding time for that rise (approximately 5.8 seconds).

This plot is all the information we need to tune according to Skogestad, which is rather neat. We find:

• Time constant tau_1 = 5.8 seconds
• Steady state gain k = delta(n)/delta(kappa) = 0.063/0.1 = 0.63

According to the rules of Skogestad, pretending this is some kind of first/second order system, disregarding time delay and setting tau_c = 5, we get:

• Kp = (1/k)(tau_1/tau_c) = (1/0.63)(5.8/5) = 1.84
• Ti = tau_1 = 5.8 s

That's it, super simple and super fast. Tuning this plant using the trial and error approach we get Kp = 1.9 and Ti = 6.3. From a practical point of view, this is more or less identical. Both will give fairly equal performance, and both well withing specs. It's not coincidental that a well behaved system has a generator inertia constant of approximately Ta = 6, and this corresponds well to the "optimal" Ti in this case.  

A realistic system

This is a real system. It's certainly on the "fuzzy" and oscillating side, but far from unusual. Doing the same as for the idealized system, we get the plot:

We find:

• Time constant tau_1 = 9.5 seconds
• Steady state gain k = delta(n)/delta(kappa) = 0.022/0.1 = 0.22

Then with tau_c = 5 as before:

• Kp = (1/k)(tau_1/tau_c) = (1/0.22)(9.5/5) = 8.63
• Ti = tau_1 = 9.5 s

This time Ti seems reasonable, and is in line with what the trial and error method gives. Kp on the other hand seems way off, and it is. It will produce a system behaving like this:

The system will oscillate forever, the PID is unusable. Note also the slow servo speed, needed for safety reasons on this plant to prevent water hammer. Since Ti is OK, we can simply use step 5 from the trial and error approach to get a working Kp. Doing this, I find that Kp = 1.4 is as "good as it gets", and still be within performance specifications. A step response in set-point of the closed system, certainly does not look good, but the physics of the system does not allow for anything better.

In the last plot below, I have introduced anti-windup. It doesn't improve the closed loop response, but makes a smoother servo run (less wear and tear) and will have a much better performance in terms of handling disturbances because Ti is reduced considerably. This is done by trial and error.

For real systems, some trial and error usually has to be done no matter what, but Skogestad method is very helpful in obtaining Ti for systems with difficult behavior. An open loop step plot will also reveal useful information of the system. In the next post, a combined approach is shown.

PID tuning the simple way

Tuning that works "every" time

Throughout the years I have tried every method I could find, from text book classical analytical s-plane methods to trial and error types. It's mostly turbine governors I have been working with, both in theory and in practice. The only method I have found that works almost every time, in theory and on site, is a trial and error method. Sadly I cannot find the copy of the original text, nor can I remember where I found it. I have written down the method though, ages ago.

The method requires that you are able to set the set-point arbitrarily and can read off the process value. A graphical representation is of much help rather than just numbers, although numbers works as well. The method is also described in the manual of LVTrans. The method requires no knowledge of the system, although such knowledge certainly helps.

The first thing to do is to stabilize the system. This can usually be done by setting Kp to 1 or less, and Ti to 10 or more and Td to 0. Then set the set-point (SP) to a desired value of which you are going to tune around, the operational point OP, and let the whole system settle. It's useful to have some integral value here to get the SP and the process value (PV) equal, but not strictly necessary. When this is done, and looks OK and steady state, the tuning can start.

1. Set Ti to a really large number, or deactivate further integration if possible. In parallel form, Ki can be set to zero. Keep Td = 0. Let the system stabilize if it isn't already.
2. Do steps in SP of approximately 10% up and down from the OP and observe the PV.
3. Increase Kp bit by bit until an overshoot of the PV of about 15% is observed. That is 15% of the size of the step in SP, not 15% of the total PV.
4. With the Kp from 3, decrease Ti bit by bit until the overshoot increases to around 30%.
5. For a PI governor: Decrease Kp bit by bit until the overshoot is 5-20% or looks OK - finished.
6. For a PID: Increase Td bit by bit until the overshoot is 5-20% or looks OK - finished.

That's it. An "optimal" PID will have some overshoot, but for some applications an overshoot may not be wanted. Then simply adjust Kp or Td in step 5/6 until little or no overshoot is seen. By coincidence these steps looks very much like the graphical animation in the Wikipedia article.

Note. This method makes it possible to tune most PID controllers in most systems, but doesn't necessarily create the "best possible parameters" whatever that may be. For turbines the specs for the governor are in frequency domain (a separate analysis) and the ability to follow changes in SP and changes in frequency. Any tuning that makes the system behave according to spec is good enough. If no formal specs exists, then this method results in no better or worse parameters than any other method. It does however create good and stable results almost every time, and with little knowledge of the system.

On a second note, the tuning of systems where servos saturate can be tricky. By saturation I mean the output is restricted in amplitude and speed. Some anti-windup algorithm should then be included for best behavior of the PID controller. This should be in effect for step 4 and 5/6 if it isn't already. Most commercial PID's have anti-windup implemented.

On a third note, a pure 1st order system may not give any noticeable overshoot in step 3. If that happens, then you know it's a first order system. A typical example is PID servo positioning. For such systems you can set Kp to almost "whatever you want", and adjust Ti to achieve a response with no, or desired overshoot. If performance is specified, you simply adjust to achieve the wanted performance. In this case, it is however a good idea to dive a bit further into the theories of PID tuning, starting for instance here.

Background

Tuning a PID governor is the process of determining the proportional gain, the integral gain and the derivative gain so the process behave the way you want it to. A key element is "behave the way you want it to".

If the stability is important, then there could be industry standards specifying the stability gains. Typically there are gain margins and phase margins that are standardized. These are obtained in frequency domain analysis, either directly or through FFT analysis. An industrial standard could also specify certain time parameters, how fast the PID should be, how accurate, the overshoot and so on. If no specifications are given, then the usual way is to get the PID so "optimal" as possible while keeping it stable.

There are therefore no universal descriptions of what a proper tuned PID actually is. It will be dependent on the application, the industry and the very process it is governing. The only common denominator is stability of the system as a whole. Theory and practice are usually different when working with PIDs, especially when it comes to derivative. This is sometimes referred to as realizable forms of the PID vs mathematical/theoretical forms. 

Two popular versions of the PID exists. The standard form and the parallel form. From Wikipedia the standard form is:

and the parallel form is:

where

As long as the relations between the constants (above) are remembered and taken into account, these two are equivalent mathematically. The tuning procedure shown above uses the standard form. The parallel form could also be used as long as the relations above are taken into account. Both these forms above happen to be mathematical forms rather than realizable forms, but that is of little relevance here.

Neglecting the convective terms

General

A thing that was done when deducing the governing equations was to neglect all convective terms. The details of this is rather elaborate, and involves lots of equations (different forms of Navier-Stokes equations). A full treatment is not done here, only the key points are included. The full analysis can be found in Appendix A of my thesis.

The basic principle is to make the Navier-Stokes equations from Part one non-dimensional in the same manner as was done in the Long wavelength approximation in the previous post.

Why would we neglect the convective terms in the first place? The reason is that they add considerable complexity in solving these equations. At the same time, if they happen to be unimportant for the solution nonetheless, it's only silly to include them.

We know that liquids are compressible. It's what makes transient analysis of liquid filled pipe interesting and important. However, we also know that liquids are only slightly compressible. Gas on the other hand is highly compressible. Can this compressibility be related in some way to some characteristic property, like for instance fluid velocity? Yes it can.

The key property describing the compressibility in relation to transient piping flow analysis is the Mach number. The Mach number M is defined as:

M = V/a

Where V is the fluid transport velocity and a is the wave propagation velocity (commonly known as the speed of sound). The speed of sound in water is typically around 1400 m/s, while the speed of sound in air is about 340 m/s at room temperature and atmospheric pressure. The (transport) velocity of the fluid in a pipe is for liquids around 1-10 m/s at max, while for gas it can be 100 m/s or more in some cases depending on application. For choked gas flow the Mach number is reached. It's safe to say that for liquids in pipes, the Mach number is very small, while this is certainly not the case in general for gas flow (although it could be the case in some specific circumstances as shown later).

Continuity

We make the NS equations non-dimensional carefully choosing appropriate characteristic properties. The continuity equation (in 2D cylindrical coordinates) pops out to become:

Here we see that the order of smallness for all convective terms is related to the Mach number squared. This is beneficial. If the Mach number already is small, a typical value is 0.004 for water in a pipe, then this squared becomes 0.000016, which is even smaller. A table for all the convective and space derivative terms is set up:

We see that the order of smallness for all the convective terms is very small, as long as M << 1, and they can be neglected. This means that when M << 1, then the convective terms can be neglected both for 1D and 2D equations in the continuity equation. Even if the Mach number is not super small, it has very little effect in the continuity equation. For typical flow in (large) gas lines, the Mach number is around 0.1 -ish, which means the order of smallness becomes 0.005 -ish (0.5 %) which is well inside the accuracy we can hope to achieve in a calculation.

Momentum

Doing the same for the momentum equation, the following non-dimensional equations in axial and radial direction pops out.

Now, this is unfortunate. The momentum equation is not related to Mach number at all. No matter what Mach number we use, the order of smallness of the convective terms becomes O(1) and cannot simply be neglected. There's nothing we can do about it, it's just the way things are. 

But, there are a few other things we can do. One thing that can be done, is in combination with the continuity equation to force the convective terms to be small. We say that the transport velocity is small and negligible. This is called the acoustic approximation, and is used for acoustic analysis (typically in 3D). In 1D this will simply produce the ordinary water hammer equations, but without the friction terms. A kind of side effect is that the ordinary water hammer equations can be used to calculate acoustics in gas pipes, and very accurately so, as long as the long wavelength approximation is taken into account. This was actually the first application of LVTrans, calculation of acoustics in gas risers.

What is done in my thesis is to apply the long wavelength approximation, L/D >> 1 to make the whole thing 1D. The continuity and momentum equations then become:

A.13 and A.14 is combined to obtain:

Now the Mach number is back, and we can compare the time derivative terms for smallness:

This shows that the in the combined set of equations for continuity and momentum, the residing convective term becomes insignificant when M << 1. This means that M << 1 does not in general remove convective terms in the momentum equation at all. It's only when applying the long wavelength approximation (1D flow) that the combined set combined will show the convective terms to be insignificant, and only when M << 1. the final non-dimensional momentum equation becomes:

Equations for (slightly compressible) gas flow

For the sake of completeness, it is worth considering the more general form of the non-dimensional equations where only the long wavelength approximation have been introduced (1D flow). We simply remove all the radial terms to obtain (removing all the subscripts for clarity):

These are two equations governed by two constants: The Mach number and the Reynolds number. It must be noted that the variation of density with respect to pressure is constant for these equations, the equation of state is over-simplified. The following non-dimensional variables are used in this post.

The long wavelength approximation

The long wavelength approximation is not as much a "mathematical trick" as it is a fundamental part of nature itself. The meaning of that is a the liquid in a pipe of some length compared to the diameter cannot only be treated as a 1D phenomenon, it essentially IS a 1D physical entity. Treating it in full 3D certainly can be done, but the 1D nature will completely overshadow any added wanted increase in accuracy from the 3D approach. The only valuable output from a 3D approach, at least for engineering purposes, will be the 3D approach showing it is essentially a 1D phenomenon.

This has all to do with inertia and the characteristic length and diameter of the pipe. A constraint is also that any externally forced frequency must be lower than first radial resonance frequency, which normally starts in the kHz range.

Why this is so is explained next. We look at equation 3.6 and 3.7 from part 1.

Equation 3.5 is not needed, because the constraint is frequencies below the first radial filled mode. In other words, what we are looking for is when we can use that constrain, and which parameters decides it. This is done using non-dimensional analysis, so we can see the relative magnitude of radial modes vs axial modes. The following non-dimensional variables are used:

These are inserted into B.1 and B.2 and we get:

These two equations are very similar. In fact, they are exactly equal, both mathematically and in magnitude (physically) when:

This means that if we have a pipe, and the L/D ratio is 1/2 (it's half as long as it's diameter), then radial and axial inertia will have equal effect on the dynamics of the system. Is a pipe where the length is only half as long as the diameter really a pipe? Hardly. So we introduce the constraint that L/D >> 1, or conversely that D/L << 1. We then see that the first term in Equation B.6 becomes very small (much less than 1), but nothing happens to the other terms in B.5 or B.6. The larger the ratio L/D becomes, the smaller the inertia term in B.6 becomes.

With L/D = 5, the effect or the radial inertia is already an order of magnitude less than the effect of axial inertia. Another way to see this is to derivate Equation B.5 with respect to r, and B.6 with respect to x, and combine them to get:

This equations shows the magnitude of the change of velocity in radial direction with respect to change in x and the magnitude of change of velocity in axial direction with respect to change in r. When inserting L/D  >> 1 (D/L << 1), the first term simply vanishes. Equation B.7 then becomes:

This equations say that the change of axial velocity with respect to r is zero, thus the axial velocity is constant across the cross section.

A key element of this analysis is of course non-dimensionality, which is a very powerful, and often necessary tool in all of engineering. Non-dimensionality enables us to extract the physical essence that we otherwise would not be able to see.

One question remains. What exactly is L/D in terms of a numerical value so that L/D >> 1? Is it 5? 100? 8.67? There is no simple answer. For instance, a complex piping system may have a few very long pipes, a few rather short ones and lots of pipes with lengths in between. We will then be more interested in the characteristics of the total system, than in each individual pipe. For the vast majority of systems, the longest pipes will decide the total characteristics. If the average L/D for all pipes is larger than say 10, then the 1D phenomenon will be the important one. If the average L/D is smaller than say 2, then this is hardly a piping system anymore, but a complex 3D shape and should be analyzed as such. The boundary conditions are also very important in this respect, as well as the reason for doing a transient analysis in the first place.

This is the case, not only for transient analysis, but also for steady state analysis. After a sharp 10 degree bend, the flow is not considered to be fully reorganized before after 10 diameters in length and so on. This doesn't mean that 1D analysis cannot be done. In most cases it only means that we have to be aware that the accuracy may not be as good as it normally would be. Another question is why would you do a transient analysis in the first place? As long as it is conservative and roughly accurate, it would be difficult to justify the added cost (and added other uncertainties) of a full 3D analysis, unless other considerations, typically FSI, also becomes part of the picture, and failure of the system is unacceptable.

Governing equation for transient pipe flow, part 4

From the last part we ended up with these two equations:

These equations are on a form that is suitable for FSI analysis when combining them with the equations of motion for piping systems (longitudinal, rotational and bending vibrations). They could also be used as is for pure hydraulic transients, but the normal starting point is to use pressure head, H, instead of the pressure, P, and velocity, V instead of flow rate Q (even though Q usually is eventually used in the end) and the diameter, D, instead of the radius, R.

Neglecting the Poisson coupling (sturdy piping structure), the last term in equation 3.19 and reorganizing we get:

Here the wave propagation velocity, a, becomes from the outset (for thin walled pipes):

For most cases this is a good approximation, but Wylie & Streeter have elaborated this for lots of different conditions in their books. Based on experience and practicalities however, the wave propagation velocity is a very fuzzy thing to determine. It can hardly be experimentally verified for any real piping system within +- 20%, mostly due to small, and changing, amounts of gas content, gas pockets and so on. Wylie & Streeter writes a lot about that too, and in the MOC (Method of Characteristics) it is used as a parameter that is changed within limits to discretize large systems. Even though it can be calculated theoretically at very high accuracy, that accuracy simply does not exist in real systems. Having a good starting point, is a good idea nonetheless though, because the difference between polymer based pipes and steel pipes can be more than 100% for instance.

What's missing now is the pipe friction. Diving into this, one discover it is super complex material. To this day, no good general transient friction term exists as far as I know. I have written a lot about it in my thesis. For steady state, Darcy-Weisbach and the Moody diagram is what is used. In lack of a better general non steady state model, Darcy-Weisback is also used for transient simulations, even though "everybody" knows this is far from correct. The best we can say about this, is that at least the calculations becomes conservative. Calculating water hammers that are smaller than in real life, does not normally happen. The exceptions are when FSI becomes important and for fluid cavity collapse. If this is not taken into account, then the calculations are usually not conservative anymore.

The Darcy-Weisbach equation:

Here f is the Darcy-Weisbach friction factor we normally obtain from the Moody diagram. The minus sign is there because of integration when positive x is in positive flow direction. This is simply the steady state result of the momentum equation above, and we can "reverse" it back into that equation to achieve exactly that. The final equations for liquid filled transient pipe flow analysis therefore becomes:

The wave propagation velocity, a, is as shown above.

Governing equation for transient pipe flow, part 3

 To get further from part 2, Hooke's law in cylindrical coordinates is introduced.

Then by combining Equation 3.13, 3.15 and 3.16 we get:

We have already discussed that the radial inertia term is unimportant (Long wavelength approximation). We can therefore substitute the tangential stress term for the quasi-stationary relation for a thin walled cylinder found in most text books, and assuming the pipe wall is much thinner than the pipe diameter.

A bit reorganizing, and Equation 3.17 becomes:

This is now fully 1D, and the stress term can be substituted using Hooke's law in 1D:

We now get the final fluid equations:

 

This is the same equations for transient pipe flow as can be found in any text book (when re-organized a bit), except the last term in Equation 3.19. This is the so called Poisson coupling of the fluid to axial motion of the pipe. The fluid wave propagation velocity is approximated by parenthesis before the second term (when multiplied by the fluid density), and is also essentially the same as found in most text books.

An interesting point is that when deducing the Water-hammer equations from Navier-Stokes equations, the FSI coupling pops out automatically, even for fully 1D equations. The only way to remove that coupling, is to force the axial movement of the pipe to be small compared with the pressure forces. We simply say that axial displacement of the pipe wall, u, is "small". In other words, it is assumed that the piping system is in practice 100% rigid. The other way is to set the Poisson ratio to be zero, which it isn't. The Poisson ratio of steel is a constant, and about 0.3.

For normal engineering purposes we simply say that the piping system is "rigid enough" for the transient fluid simulations to be correct. This is usually the case for the vast majority of systems. The only way the piping system is normally included, is in the calculation of the pressure wave velocity, see for instance Equation 3.18 and 3.19 where it can change the wave propagation velocity substantially. This will be correct also for thick walled pipes, except in cases where the pipe wall material becomes very "soft".

For the pipe itself, the axial motion ends up being described by:

The pressure is coupled to axial displacement by the Poisson ratio. When we say that the axial displacements are small compared with pressure forces, we can flip it around and say that the pressure forces are large compared with the axial displacements. Thus, vibrational analysis of the mechanical piping system without including the fluid transients, is therefore fundamentally wrong from the start. We see that the pressure forces are proportional to R/e, the radius of the pipe over the wall thickness. By increasing the wall thickness we will also reduce FSI effects, which fits well with intuition in this case.

For a full FSI analysis, junction coupling, bending and torsion also must be included, as well as all the clamping structures. It is a bit funny however, that transient simulations of the fluid inside the pipe can be done in most cases by disregarding the piping structure elasticity, but doing a vibrational analysis of the piping structure cannot normally be done correctly without including the fluid transients.

Governing equation for transient pipe flow, part 2

In the previous part, the so called extended water hammer equations was discussed (FSI), and a reduction of the Navier-Stokes equation was started.

Imagine a valve being slam shut at the left end. Water streaming from the right. The valve causes a pressure wave, a shoch wave (water hammer) propagating to the right at about 1200 m/s. The widening of the pipe wall causes axial stress, which will shrink the diameter through the poisson coupling. This causes the precursor wave propagating at about 5000 m/s.

The equations were 2D. We want them 1D. This is done by averaging the pressure and flow across the cross section. Mathematically this is done by multiplying with 2𝜋R, integrating with respect to r from 0 to R and dividing by 𝜋R^2.

V and P is the average velocity and pressure, and expressed as:

This is just a fancy way of saying V and P are the same all over the cross section of the pipe. This is of course not 100% correct physically, particularly not for the velocity V, but this is the simplification we must do to make them 1D.

But we had one more equation from the last part:

The assumption is that the characteristic fluid length scale, L is large compared to the characteristic diameter, D. That is L/D >> 1 under the constrain that the wavelength is "long". With this assumption, it can be shown that the radial inertia is negligible compared to the axial inertia. We can simply neglect the first term in equation 3.7. This is called "the long wavelength approximation", and is found to be correct for frequencies up to 63% of the first radial liquid filled pipe frequency (which normally is much higher than the first axial frequency). Proof of this is found in my thesis. Then we are left with the last term in equation 3.7 to be equal to zero, which simply states that the pressure is constant in the cross section. Thus "the long wavelength approximation" will also "one dimensionalize" the pressure, and remove one equation from the set.

On the pipe wall, the radial fluid velocity must be equal to the radial pipe velocity:

Inserting this, and we are left with only these two equations:

We have removed one equation. We have transformed them from 2D to 1D, but we have added a radial wall velocity term. This looks like a complication, but is a key point for the end result in part 3.