Advanced Control - System Identification

Hello guys,
we're continuing our series of threads named "Advanced control" whose purpose is to present formal methods and techniques of system identification and control. In previous threads we could see how to desing controllers for heat exchanger and tank level control. Now it is time to learn how to do system identification i.e. how to obtain (mathematical) model of process in order to gain better control. I asked our Peter Nachtwey to teach us techniques he prefers and luckily he agreed.
In next posts different identification methods for heat exchanger will be explained. To refresh memory this link:
www.controlguru.com/postpics/heatexchanger.doc
will be helpful.
Since I'm not familiar with identification methods I cannot make any detail introduction.
Peter, it's your turn...

Hold on to your seats, this will be a wild ride.

If you read through all the post on www.controlguru.com you will learn a simple way of estimating to gain, time constant and dead time. For many applications the visual estimating method is more than good enough. Pandiani want to learn a more mathematical approach that is similar to the methods that Expertune and Bestune use. This knowledge is worth $$$$ but it will take some effort or home work. Anybody that has many PIDs to tune should pay attention and struggle through this.

I am going to make two major steps.
1. manually adjust the gain, time constant and dead time to minimize the SSE or sum of squared error.
2. Using an minimization routine to minimize the SSE.

I just commited to do this so I don't have anything prepared yet. This may take a day or two. All examples will be using Scilab which is available free at www.scilab.org.

Attached is the heatexchanger data that will be required to do the system identification. It is basically a trend file in text.

Manually Minimizing the Sum of Squared Error or SSE.

Attached is a Scilab program that allows me to find the Gain, time constant and dead time of the heat exchanger. I also included a .gif that shows my results and the original heat exchanger data.

The goal is to minimize the error between the model's EV and the PV from the actual data. The EV is the model's estimated value of the temperature as a function of the actual control output data. The sum of squared error is the sum of all the errors between PV and EV squared. One can manually change Kp, Tp, and the dead time and run the simulation. If you changed the plant parameters in the correct direction the SSE will drop. If you made a bad change then the SSE will increase.

Note the SSE for the plant parameters given by the Control Guru is 2.127. I was able to reduce SSE to less that 1.3. It took a few minutes. You will see the Control Guru's numbers can be improved.

I apologize for not getting the slider bars to work. I am learning Scilab as I write these programms. You must edited and load the plant parameters manually. I will still work getting the slider bars to work but most will get excellent answers in about 15 minutes.

There are a few on the forum that remember doing this to tune systems using excel a couple years ago.

Peter,

The sliders really would help. However, after looking at the demo program that has them (TCL - Color Demo) I can see why you haven't gotten that far.

Do you have the book referenced in the comments or is that comment from something else. It looks useful.

Here's an Amazon link
http://www.amazon.com/Modeling-Simulation-Scilab-Stephen-Campbell/dp/0387278028

I believe you left your optimized tuning parameters in the file so not much to do but tweak them here.

Can't wait for the minimization routine....

OOPS!!! My mistake.

Yes, I left the changed parameters in the file which ruins the fun. My bad. Can you improve upon my plant paramters? I didn't try that hard. However, you can see that this works with minimal effort and it isn't magic. I could have done this on an Excel spreadsheet. However, this example only has three variables. Five is a little harder but companies like Boeing have systems where they have 40 parameters to find. 40 is far too many to be adjusting by hand. Therefore a more methodical or mathemagical way is required. That will be the second part.

I have the book.

The second part will use the lsqrsolve function. I choose to use the lsqrsolve because it can find values for parameters even if the system is non-linear. I will not try to explain how the lsqrsolve function works except to say that it has an almost optimal way of adjusting all the parameters at the same time. It takes only about 10 iterations for a simple 3 parameter system. A five parameter system would take about 20-30 iterations. The real trick, or problem, is the dead time because it can't be adjusted in a continous way.

Thanks for that Scilab script Peter.
I recall we discuss these issues before and one procedure would be to write loops that start with initial parameter's values and changing one at a time until SSE cannot be further minimized. Then same procedure is repeated for second and third parameter. Once whole procedure is over, obtained parameters are initial parameters for new iteration.
I wonder if this way is used in some software solutions or maybe more advanced things are applied. such as that famous lsqrsolve function...

Extra Credit

Pandiani said:

I wonder if this way is used in some software solutions

Click to expand...

No, there are more better methods you suggest below. I just wanted you to see what we are trying to accomplish.

or maybe more advanced things are applied. such as that famous lsqrsolve function...

Click to expand...

Yes, there are many methods. Some are very simple. Read this for extra credit.
http://seal.web.cern.ch/seal/documents/minuit/mntutorial.pdf
Just try to get the general feel for what they are trying to do. It will keep you 'entertained' while I figure out how to use the lsqrsolve function and sliders bars. It has nothing to do with PLCs. It is pure math. However, an engineer is always trying to find optimal solutions. Optimizing is applied to many other fields too like economics.

More good info on this topic.
http://cow.physics.wisc.edu/~craigm/idl/fitting.html

FOPDT system identification

It took a while. I definitely took longer than it should. I must have made every mistake. It didn't take long to get the lsqrsolve to work but to make it find the dead time work took a while. I should have just used a brute force method.

What I found out is this. The deadtime is not continuous so it can't be optimized using the lsqrsolve. lsqrsolve needs to find the derivatives at each point and that can't happen with a dead time that changes in steps.

The good part is that the program works well and it got the SSE down to less than 1. That is a big improvement over the manual SSE of about 1.3. Scilab can also calculate the system identification quickly. Now one can use this program along with the formulas from the control guru site to tune your FOPDT PIDs.

Note there is also a Scilab 4.1 now. The debugging tools are better.

Now for the nitty gritty. The FOPDT routine is the heart of the evaluation. The lsqrsolve thinks it is just optimzing Kp, Tp, and PVss but interally the coefficient are generated and a complete simulation is done. The lsqrsolve tries minimize the returned errors. If your system is a SOPDT ( second order plus dead time ) you would need to replace the FOPDT function with a SOPDT. The parameter list that is contains the plant parameters can be easily changed and extended. Some aircraft companies have models with 40 or more parameters. One can even make the plant gains change as a function of temperature like it really dies on a heat exchanger.

Upgrade and some bug fixes

Attached is a .zip file that has the latest verstion of the AutoTuneFOPDT.sce. I fixed a bug where the program would not find the correctly solution the first time it was executed.


MinSSE is the minumum sum of square errors. The lsqrsolve tries many combinations of Kp, Ti, DT and PVss in order to find the combination that matches the data the best. The values that match the best will have the minumum sum of square errors.
MinDT is the minimum dead time or the deat time that yield the minimum sum of squared errors.
Minp or minumum parameters is the Kp, Tp and PVss that provided the minimum sum of square errors.

One can see by the graph that the first order identificatin gets close. It is more than good enough to use in calculating the control K, Ti and Tp. This program takes only a about 10 seconds to run.

I think I will add a gain calculator function next and then a second order plus dead time system identification with a gain calculator . That will be easy. I bought a book on Tcl/Tk to help with the HMI part but I didn't have time to read it yet since I was out in the field playing with a rotary shear the week before Christmas. HMIs are the hard part.

Thank you for your time Peter. I went through your scilab code several times, and played with lsqrsolve function myself until I generally understood what is going on. This L-M method minimize sum of squares and I expected that results would be very close to those obtained in controlguru articles. To refresh our memory, parameters were: Kp=-0.53, Tp=1.3min, DT=0.8min.

This scripts gives following results:

Kp=-0.528, Tp=56s=0.93min, DT=66s=1.1min.

I tried to change initial values of parameters, but every time I got same results, which means that procedure converge to optimal parameters no matter how initial values are far from optimal.

I suppose that you used same test data as used in controlguru examples.

There, I have found your comments, I'll quote one:

1. „Peter Nachtwey Says:
   March 15th, 2006 at 12:19 am

I guess we will never agree. I don’t even agree with myself. I have computed the plant parameters two ways. The first method uses the least squares system identification.

Kp = -.551
Tp = 1.45
deadtime = .5333 min
Steady state temperature = 161.4 degrees
Mean square error = .0059

The second method uses quasi Newton minimization function to minimize the square error between the logged data and the estimated data.

Kp = -.5721
Tp = 1.56
deadtime = .8 min
Steady state temperature = 162.3
means square error = .007772

I all cases the first order tuning is very forgiving of modeling errors when simulated against a first order model. Next I will try the first order tuning against the more accurate second order model. So far, the second order model fits the data much better. The mean squared error is .001331“



I don't know exactly how L-M algorithm works internally, but I know it use derivation (approximation) and I know basics (from wikipedia). I don't know how to explain the fact that time constant and dead time are different from what one could expect from the articles.



Especially I'm interested in dead time, since it's killer of good control. I'm asking this because I was thinking to approximate deadtime using Taylor series. It is known fact that virtually any function who has derivation can be represented using Taylor series.

For example, exp(x)=sum(x^n/n!), when function is approximated near zero. I have been thinking maybe to approximate dead time with finite members of Taylor series, but since

exp(-sT)=1/sum(x^n/n!), this will result in higher order system, and it's discrete equivalent would be very hard to find.

Pandiani said:

This scripts gives following results:

Kp=-0.528, Tp=56s=0.93min, DT=66s=1.1min.

I tried to change initial values of parameters, but every time I got same results, which means that procedure converge to optimal parameters no matter how initial values are far from optimal.

Click to expand...

Yes, and that is a good thing. This just show the LM algorithm is very robust and an algorithm you can trust. Notice that the sum of square errors is much lower than what I had found before and what the control guru had found.

I suppose that you used same test data as used in controlguru examples.

Click to expand...

Yes, but I wasn't using the LM algorithm before. I was using the minimize function of my Mathcad. This Scilab program is the only the second time I written a program that uses LM algorithm.

I don't know exactly how L-M algorithm works internally, but I know it use derivation (approximation) and I know basics (from wikipedia).

Click to expand...

Basically it finds the partial derivatives in all dimensions so all the variables can be adjusted at once instead of doing one at a time using the manual method.

I don't know how to explain the fact that time constant and dead time are different from what one could expect from the articles.

Click to expand...

The dead time cannot be simulated in a continous manner. In a digital simulation it is impossible to calculate the derivative the same way as the continous or analog variables such as Kp and Tp. There are ways but you must modifiy the LM algorithm itself and we can do that because it is built in function to Scilab.

I/We use the LM algorithm for tuning our motion controllers. I have the implemented the LM algorithm in C starting with some code from the source forge. It is very fast. Minimizing is an art in itself.

Especially I'm interested in dead time, since it's killer of good control. I'm asking this because I was thinking to approximate deadtime using Taylor series.

Click to expand...

That isn't the way to do it. The common technique is to use a Pade Approximant.
http://mathworld.wolfram.com/PadeApproximant.html
Look at the different ways exp(x) can be approximated. I think the version the control guru uses is (2+x)/(2-x) assuming x is small relative to the system time constants. Now the exp term is gone and you can work in the s domain easily. I will post a link when I get home. I have the proof of how Kc, Ti and Td are calculated from the system ID. The only trick is substituting the Pade Approximant for exp(x). If you want a better approximation then use more terms.

There are methods such as the Smith Predictor that can compensate for the dead time by controlling the model without dead time and not the actual system. The actual system is used to update the model so it doesn't drift from reality.

In the digital world exp(-sT) is approximated by putting time delays on the inputs.

y=A*y(n-1)+B*x(n-int(DT/T))

where DT is the dead time and T is the sample time. If you want a better estimate you must make T small.

You ask good questions. You ask the same ones that bothered me for so many years.

Peter Nachtwey said:

You ask good questions. You ask the same ones that bothered me for so many years.

Click to expand...

Haha, I heard somewhere that stupid questions don't exist, only stupid answers

Yes, there are a lot of questions that bother me today. When I was in college, I must admit, I was primary interested to pass exams as soon as possible (and with as better mark as possible) and simply many questions didn't "occured" to me. But this forum and web sites such as www.controlguru.com helped me to start "thinking" and asking questions. Also, it's easier to ask this way, since nobody will tell me "how the hell you don't know that".

When I was student I designed controllers by simply playing with Matlab and it's rltool, I have added zeros and poles to achive desired performance, but in example such as tank level control, I learned that everything has its price. I learned that I can easily design controller wich will provide desired output performance, but controller's output would not be feasible in real world. Yes, I heard about windup and similar but simply didn't have chance to fight it. Now, I can only assume how much I don't know about these things.

Nevertheless, as usual your answers Peter open more questions (sometimes I'm afraid to ask question ).

Peter Nachtwey said:

There are methods such as the Smith Predictor that can compensate for the dead time by controlling the model without dead time and not the actual system. The actual system is used to update the model so it doesn't drift from reality.

Click to expand...

Searching the web I found model with Smith predictor as shown in the attachment.

In direct branch, there is PID controller and process, however in lower branch there is FO process model and dead time separately. Signals from these two branches are subtracted to eliminate dead time. If I understood good, the whole idea is to design good PID controller using model of process without dead time. If we design PID controller for desired response of system without dead time, then that PID is also optimal for system with dead time. I must admit I don't understand how this would be done in reality. If you know something more, please explain.

Thanks in advance.

Pandiani said:

Haha, I heard somewhere that stupid questions don't exist, only stupid answers

Click to expand...

You don't really believe that do you? There are good questions, stupid question, better questions and questions that shouldn't be asked depending on the situation.

Yes, there are a lot of questions that bother me today.

Click to expand...

I think that is good.

When I was in college, I must admit, I was primary interested to pass exams as soon as possible (and with as better mark as possible) and simply many questions didn't "occured" to me.

Click to expand...

Same here.

But this forum and web sites such as www.controlguru.com helped me to start "thinking" and asking questions. Also, it's easier to ask this way, since nobody will tell me "how the hell you don't know that".

Click to expand...

Why don't you go to sci.engr.control
http://groups.google.com/group/sci.engr.control?lnk=oa&hl=en
and ask how one does the calculations for finding the system identification for a first order or second order system just to see the answer you will get. You will find that if you understand the AutoTuneFOPDT Scilab program you are in a small minority that have seen how it is done. Do it!

There is nothing to prevent you from asking a question you already know the answer to just to find out what you know relative to everybody else.

When I was student ...

Click to expand...

Keep at it. I am much older. I had to learn on my own with no one to ask. It takes time. Some algorithms were only discovered in the early 90s. ( see Ackermann's equation ). There wasn't an internet with www.controlguru.com then where I could ask anybody. Even so, the control I do is much different because it is motion control. The theory is the same but the practical aspects are much different. I don't get involved with heat exchangers or oven except as a mental exercise.

Searching the web I found model with Smith predictor as shown in the attachment. ....

Click to expand...

The diagram is very good. Suppose we take the model of the heat exchanger as an example. What if we used our calculated PI parameters to control our FO model without the dead time and also use the output to control the real system? What do you think would happen? Notice that if we assume the dead time is zero then the controller gain will be much higher. Think about. Start another Advanced Control thread for Smith Predictor.

I think we should continue with system identification. Do you think you could write the SOPDT routine that would be called from the function fct? Assume the second order model is Kp*exp(-s*DT)/((Tp1*s+1)(Tp2*s+1)) and the digital form is
y = A1*y(n-1)+A2*y(n-2)+Kp*(B1*x(n-1-nDT)+B2*x(n-2-nDT)+C
I will help if you get stuck. It will require some thought and there are a couple of approaches to the problem.

Once we what SOPDT written we can compare the SSE of both the FOPDT model and the SOPDT model to find which is the best.

Then we will be ready for the Smith Predictor. BTW, you should know that I have never written a smith predictor but it should be easy. What I propose it that we use the SOPDT model to simulate the real plant and then use our cruder FOPDT model and the resulting PI gains to control the FOPDT model. It should be interesting.

Peter Nachtwey said:

In the digital world exp(-sT) is approximated by putting time delays on the inputs.

y=A*y(n-1)+B*x(n-int(DT/T))

where DT is the dead time and T is the sample time. If you want a better estimate you must make T small.

Click to expand...

Peter and Pandiani, thanks for starting this very useful thread.

In above algorithm, how to determine A and B? Do I need to store DT/T number of samples? If I don't have enough memory to store DT/T samples but still simulate DT then what shall I do?

I want to implement Deadtime in PLC for simulation.

SVN

svn said:

Peter and Pandiani, thanks for starting this very useful thread.
In above algorithm, how to determine A and B?

Click to expand...

You must download the AutoTuneFOPDT.zip file to see how A and B are calculated. You should do a search for Advanced Control to find the threads that we have started and read them all.

Do I need to store DT/T number of samples?

Click to expand...

Yes if you are going to implement the Smith Predictor. If your dead time isn't too big then use the formulas posted at www.controlguru.com.

If I don't have enough memory to store DT/T samples but still simulate DT then what shall I do?

Click to expand...

Get more memory if you are going to implement the Smith Predictor. If you don't then use the formulas posted on the www.controlguru.com blog.

I want to implement Deadtime in PLC for simulation.

Click to expand...

That is easy enough. You must only download Scilab from www.scilab.org and run the programs that I have posted.