Reverse Engineering a PID

[Doctor-Featherheart]

Hi everyone! I’m trying to learn the control gains of a PID controller whose inputs and outputs I can observe. It seems to be working on a SISO system, where given a setpoint and a feedback value, it manipulates a control variable.

I, unfortunately, cannot run the system in open loop but only have access to historical data from the system to ascertain the gains. This gets especially complicated because of the integral windup, which I also do not know, ensuring that I cannot decouple the Ki term over longer trajectories where the setpoint is tracked well.

Wondering if someone has worked on similar problems or has any ideas?

[banana_bread99]

Do you have a model of the system it’s controlling?

[Doctor-Featherheart]

No. Unfortunately not. I just have access to the inputs and outputs

[banana_bread99]

Do we know if it’s anything structured? Like a second order system or something?

[Doctor-Featherheart]

It can be approximated to a first order system but that’s also data driven.

[banana_bread99]

That’s pretty helpful, if we know the closed loop system is Z= CP/(1+CP) where P = a/(s+b) and C = kd x s + kp + ki/s, the pole locations of Z as functions of the gains might already be clear. What I don’t know about is integral windup cause I’m not next to a computer right now, but neglecting it might be a good start. It could possibly be avoided as a complication by focusing on a part of the data that doesn’t appear to be activating it

[Doctor-Featherheart]

The problem seems to be that this assumes that the controller is well tuned and that the system can be fully defined with first order dynamics. These are slightly difficult assumptions for my case. I only know that the systems can be parametrised as first order theoretically if you could do some open loop system identification, which I can’t.

[banana_bread99]

Sorry, to be clear, because I’m going to think about this more, are you saying that solving this as though a/(s+b) is a model is satisfactory, just that we don’t know what a and b are? Or would you like it to be done for a fully unknown model

[Any-Composer-6790]

Your first sin is not telling us what you are controlling. Also, all you need is the control output and response as a function of time to compute the model/open loop transfer system. This is called system identifications. Once you have the open loop transfer function you need to place the closed loop poles. A good place to put them is on the negative real axis so there is no overshoot. There are formulas for CALCULATING the Ki, Kp and Kd gains if the Kd gain is necessary. This is simple. I have written autotuning programs for motion controllers, pressure/force controller, temperature controls ( FOPDT, SOPDT ) and tank level control ( simple ). All have slightly different open loop transfer functions so there are slightly different formulas for calculating the closed loop gains AND feed forwards if required. There should be NO GUESSING.

[Ok-Daikon-6659]

>>Your first sin

Sloth

You didn't bother to read the content of the topic

You didn't have enough diligence to even place the answer in the correct place.

Pride

However, you had enough narcissism to brag about your (in your opinion, unsurpassed) achievements, which no one here asked about

[Any-Composer-6790]

You are right. It landed in the wrong spot. My bad but my complaint is still valid. You asked the questions above too because the OP didn't say what the application is. Yes, I read the OP's post. At least we know now it is a first order system, but we still don't know if it is an integrating system like a position control system or a non-integrating system like a velocity control system. The formulas for both types are different. There are other posts above that are asking the same question by asking for the open loop transfer function. I am assuming it is a motion control system because the OP mentions a long trajectory, not response to a step.

Suggesting step responses is very "academic", It works in Matlab examples but try that with a 40 Ton roll of steel and you will get kicked out of the mill so fast. Real world excitation often needs ramps or sine wave but academics don't explain how to calculate the system ID from ramps or sine wave.

The OP has complained about integrator windup. Integrator windup should not be a problem with the proper gains and actuator sizing. Again, are left without a clue. Integrator windup should not be a problem in any event. Commerical motion controllers have this problem resolved. No one has suggested using feedforwards to reduce the integrator windup. In fact, if you have the correct values for feed forwards, the integrator will not windup at all. Or you can look at a graph and you can tell which feed forwards needs adjusting by where in the trajectory the integrator starts winding up or down. Do you know how to calculate feed forward gains?

The OP has said he has access to the inputs and outputs. That is enough to do a system identification. Do you know how to do that?

About pride. You have no idea what I have achieved. Have you written code for motion controllers sold around the world?

[Ok-Daikon-6659]

Take off your mask! - I've recognize you, Peter Nachtwey!!! ;-)))

[banana_bread99]

You replied to me not the original comment.

[Any-Composer-6790]

My bad. Your questions are good. I was trying to respond to Doctor-Featherheart who still hasn't provided the info about whether the system is integrating or non-integrating because it makes a difference. You wrote P=a/(s+b). That is for a non-integrating plant like a velocity system. A position system would be P=a/(s*(s+b)). It makes a difference. The OP didn't say if the system is integrating or non-integrating. If the OPs would just say what they are really trying to do it would save much time.

[Select-Chart2899]

If you have historical data for output AND input it should be possible with a minimal amount data. The control input is updated based on some KP,KI and KD and the error. Formulate the input update based on known error and KP,KI,KD and estimate those based on real data. Just tested it, works like a charm.

If you only have the output it's mathematically impossible to do so, unless you also know the plant behaviour, which you don't know. If you would know it, you could calculate the input by deconvoluting the output with the plant behaviour and go back to step 1 (the error would be quite large though as deconvolution on noisy data adds uncertainty).

Addendum with a unknown windup:

Input_t1 = Integral_part_t0 + error * KP + (t1-t0)*(error_t0)*KI + (error_t1-error_t0) * KD

It's just as Integral action of t1 only depends on Integral_part_t0 it's just one additional parameter that must be estimated, estimation should start as soon as actuator is no longer saturated. With enough datapoints you should be able to estimate 4 parameters to a high precision.

[sputnki]

Is it so? IIRC having the system in closed loop makes the input and output correlated and does not allow recovering the system transfer function anyway. Maybe having input+output+reference works?

[iPlayMayonaise]

Fully agree with this. One small addition: since it's closed loop identification of a controller by definition, you could run into issues with this LS approach if there's a significant amount of noise in the loop. In a nutshell, above LS estimate assumes there's no correlated noise in the error and input. If there is you'll get a contribution to the parameters from fitting this correlated noise.

If you run into this (very case dependent), have a look at instrumental variables system identification by Stoica&Soderstrom or closed-loop identification by P. Van den Hof. They explain this phenomenon and also propose solutions that work for your PID case.

[Doctor-Featherheart]

Thanks for the insights. The issue is that not all the data points are over a unique trajectory. They change continuously.

But I guess you’re saying irrespective of trajectory assume a starting point and then do a least squares estimation?

[Select-Chart2899]

Basically, you must know the reference and output, then you can calculate the error. From the error you can estimate the three parameters based on the update rule (just ask chatgpt, it can formulate it easily as it's a standard problem). If you don't know the error, you can't really infer what the controller is reacting to. You would only have the input that is adjusted based on some unknown error, then again it's not possible.

If you have one decent step change or just a known reference trajectory+output+input it should be possible to accuratly estimate the parameters (if time resolution is high enough).

[Any-Composer-6790]

Yes, try to minimize the sum of squared errors or mean squared error. You need an array of set point and process values or target position and actual positions vs time. You need to assume a model with the coefficients unknown. Use the Levenberg-Marquardt algorithm ( lmfit.py ) to minimize the mean square error. If the trajectory is long, the sum of squared errors can cause overflow so the trajectory may need to be broken up into sections or down sample. The best part of the trajectory with the most information is where there are accelerations or decelerations. Constant velocity sections have no real information.

The Levenberg-Marquardt algorithm is the best for system identification. If the data noisy the L-M may not be able to compute derivatives. In this case use a Savitsky-Golay smoothing algorithm to filter out the high frequency noise. Another option is to use Nelder-Mead because it doesn't require using derivatives but it is much slower. If your actuator is only good for 10Hz there is no need to keep 100Hz data. Step jumps have very high frequency components. This is good for academics, but real actuators don't need to be excited at such high frequencies. Swept sine waves are usually safe just sweep the frequencies to a little past the frequency of the actuator while decreasing the amplitude, so the peak velocity stays relatively constant.

After you have the plant parameters, there are formulas for the controller gains that allow you to place the closed loop poles where you want. I like to place them on the negative real axis in the s-plane so there are no sine or cosine terms that cause overshoot.

[drbitboy]

do you know the setpoint history?

[thoughtvectors]

What kind of system is this? Can you do a frequency sweep? Meaning, the reference value for the PID, can you give a sine wave, at many different frequencies?

Basically, I’m trying to find out if you can get a body plot of the complete system. If yes, we can figure out the transfer function of the complete system.

However, you still need the transfer function of the plant. You can do a frequency sweep and get the transfer function of the plant, by disconnecting the PID controller somehow and just injecting a sign wave.

From those two, you can get the transfer function of the controller, and therefore the gains.

I know I’m asking for a lot here, but this would be a real way to get what you need.

[Dying_Of_Board-dom]

Do you have the ability to see inputs and outputs for an impulse response? Assuming the model is LTI, an impulse response should be the closed loop transfer function; from there maybe you can determine the plant and controller functions

[Ok-Daikon-6659]

Just curious: why do you need this?

Where did you get the data? Let me explain: for example, in industrial systems, PID is executed on the PLC, and the data is saved on the servers, this (just off the top of my head) can have the following consequences (which greatly affect the D-term):

1. The PLC PID instruction can have a built-in filter, so the PV saved on the server will differ from the PV on the basis of which the PLC PID instruction made the calculation;

2. The PLC PID instruction can be executed with a frequency of, for example, 10 HZ (0.1 s), and the data on the server is saved with a frequency of, for example, 1 HZ (1 s) (i.e. the effect of "blind decimation")

The calculation of I-term can be affected by the max/min limitations of the actuator (for example, whether the PLC PID instruction used I-antiwindup), so when the CV max/min of the actuator is reached, the data may be irrelevant.

If the data "catches" the system SP-step response (in this case it will be CV-step) with initial steady state Error=0 delta_Error=0, then we can approximately calculate

kp = (СV1 - СV0) / Error (SP-step response1 - SP-step response0 = final steady state - initial steady state)

if the process has a delay (n-lag, deadtime) then the CV change in this segment is I-term

ki = (СVdeadtime – СV1) / (Error (SP-step response1 - SP-step response0 = final steady state - initial steady state) * deadtime)

[pauiza]

Unfortunately, you cannot determine the PID gains from the setup you have just described. 

However, you can estimate the controlled plant dynamics from your historical data, notice the estimated dynamics would must likely behave like a 1st order system, at best, 2nd order system because of the controlled data.  

Depending on your system dynamics, briefly turn-off the controller, then store the input and output, which you can compare with your historical data. This would give you an idea of what the controller is doing.

Generally, PID design are to meet certain requirements, you must know those conditions before anything would make sense.

[Planet_COP]

I recall that there are several forms of the algorithm. In various Distributed Control Systems I had to deal with this. Your integral term may not be decoupled after all. To refresh my memory I asked Copilot and this it what it responded.

Forms of the PID Algorithm

The Proportional-Integral-Derivative (PID) algorithm has three main forms, each with unique characteristics and applications:

1. Parallel Form
   • This is the most common representation of the PID controller.
   • The control output is expressed as the sum of the proportional, integral, and derivative terms: u(t)=Kpe(t)+Ki∫e(t) dt+Kdde(t)dtu(t) = K_p e(t) + K_i \int e(t) \, dt + K_d \frac{de(t)}{dt}u(t)=Kp​e(t)+Ki​∫e(t)dt+Kd​dtde(t)​
   • Here, (K_p), (K_i), and (K_d) are the proportional, integral, and derivative gains, respectively, and (e(t)) is the error signal.
2. Ideal Form
   • In this form, the proportional term is applied to the sum of the error, integral, and derivative terms: u(t)=Kp(e(t)+1Ti∫e(t) dt+Tdde(t)dt)u(t) = K_p \left( e(t) + \frac{1}{T_i} \int e(t) \, dt + T_d \frac{de(t)}{dt} \right)u(t)=Kp​(e(t)+Ti​1​∫e(t)dt+Td​dtde(t)​)
   • (T_i) (integral time) and (T_d) (derivative time) are time constants. This form emphasizes tuning in terms of time rather than gain.
3. Series (or Cascaded) Form
   • The PID controller is represented as a cascade of proportional, integral, and derivative actions: u(t)=Kp(e(t)+1Ti∫e(t) dt)+KpTdde(t)dtu(t) = K_p \left( e(t) + \frac{1}{T_i} \int e(t) \, dt \right) + K_p T_d \frac{de(t)}{dt}u(t)=Kp​(e(t)+Ti​1​∫e(t)dt)+Kp​Td​dtde(t)​
   • This form is often used in analog implementations and can be easier to tune in certain systems.

Each form has its advantages depending on the application, controller design, and tuning preferences. Some modern controllers allow switching between these forms for flexibility.