Advanced control

Pandiani · Sep 11, 2006

I use Matlab student version, but I'll ask friend with ADSL to download SciLab
and it will take some time for me to learn it.
I hope user interface is similar to Simulink's.
I made Matlab model of heat exchanger http://www.controlguru.com/?p=45 using
first order + dead time model.
So I recommend to all interested in this discussion to read related articles from
www.controlguru.com
By following procedure described on web site provided above:
heat exchanger can be approximated with First order + dead time system.
Since Simulink use system model expressed in therms of transfer functions
(Laplace transform) I'll make a little introductionary first.
Beacuse I cannot type equations here, I placed text in pdf file in attachment.
Thank you for your time.

P.S. Since english is not my native language, I appologize for all spelling and semantic errors. Please feel free to correct me where is that case, I'll benefit from that.

Peter Nachtwey · Sep 11, 2006

First we must get the simple things right.

Pandiani, there is a problem with fig4. The error should come from between the summing junction and the gain. What you call the output is not right either. That should come from where your error comes from now.

What is value in figure 5? If you are copying the Control Guru example then the left side of the graph should be label temperature SP and PV.

Fig 7 doesn't is wrong for much of the same reasons as fig4.

I would plot the output the setpoint and and PV.
See my Control Guru article http://www.controlguru.com/?p=80#more-80
Below is the file I sent to the Control Guru before editing.
I will say more later. I must go to work.
Notice that I use Mathcad. Your PI graph should look like mine.

Pandiani · Sep 11, 2006

You're right Peter, I stand corrected.
Here's my response (with corrections)...

Peter Nachtwey · Sep 11, 2006

A Scilab Example

Those that want to play along should copy and paste the program below into their Scilab editor and load it.

Code:

// SIMPLE HEAT EXCHANGER CALCULATION FROM CONTROL GURU
 // [url="http://www.controlguru.com"]http://www.controlguru.com[/url]
 
Kp=-0.533		// PLANT GAIN, UNITLESS
Tp=1.3*60		// PLANT TIME CONSNTAT IN SECONDS
DT=round(0.8*60) // PLANT DEAD TIME IN SECONDS
T=1			  // UPDATE TIME IN SECONDS
PVss=140-39*Kp   // STEADY STATE TEMPERATURE

// CALCULATE THE MODEL COEFFICIENTS

A=exp(-T/Tp)	 // DISCRETE TIME TRANSITION MATRIX WITH ONLY ONE ELEMENT
B=Kp*(1-A)	   // DISCRETE TIME INPUT COUPLING MATRIX
C=(1-A)*PVss	 // OFFSET OR SYSTEM BIAS SO SYSTEM GOES TO PVss WITH NO CONTROL

// CALCULATE THE PI CONTROLLER GAINS

Tc=max(Tp,8*DT)  // DESIRED CLOSED LOOP TIME CONSTANT
Kc=Tp/(Kp*(Tc+DT)) // CONTROLLER PROPORTIONAL GAIN
Ti=Tp			// CONTROLLER INTEGRATOR TIME CONSTANT

for n=1:3600	 // SIMULATE 60 MINUTES OF TIME.  n = SECONDS

  // INITIAL CONDITIONS, NOTE SCILAB ARRAYS START WITH 1
  
  if n == 1 then
	SP(1)=140;		// INITIAL SET POINT
	PV(1)=140;		// PROCESS VARIABLE, TEMPERATURE = 140
	CV(1)= 39;		// CONTROL VARIABLE, OUTPUT % = 39%
	CVi(1)=CV(1);	 // THE INTEGRATOR CONTRIBUTION TO THE OUTPUT.  39%
	CVp(1)=0;		 // PROPORTIONAL CONTRIBBUTION TO THE OUTPUT
					  // ASSUME INITALY THERE IS NOT ERROR
  end

  // GENERATE THE SP AS A FUNCTION OF TIME
  
  if n < 25.5*60 then
	SP(n)=140;
  else
	SP(n)=138.4;
  end

  // GET THE DELAYED VALUE FOR THE OUTPUT OR CVd
  
  if n <= DT+1 then
	 CVd(n)=CV(1);		// IF THE DELAY IS PRE-HISTORY THEN USE THE INITIAL VALUE
  else
	 CVd(n)=CV(n-DT);	 // DELAY THE OUTPUT BY THE DEAD TIME
  end

  // UPDATE THE PV USING THE SIMPLE FORPDT MODEL
  
  PV(n+1)=A*PV(n)+B*CVd(n)+C;
 
  // CALCULATE THE PI CONTROLLER OUTPUT
  
  Error(n) = SP(n)-PV(n);
  CVp(n+1)=Kc*Error(n);				// CALCULATE THE PROPORTIONAL CONTRIBUTION TO THE OUTPUT
  CVi(n+1)=CVi(n)+(Kc*T/Ti)*Error(n);  // CALCULATE THE INTEGRATOR CONTRIBUTION TO THE OUTPUT
  CV(n+1)=CVi(n+1)+CVp(n+1);		 // SUM THE P AND THE I OUTPUT TERMS

end // for i

PV($)=[];	// REMOVE THEN n+1 INDEX SO THE ARRAY SIZE AND INDEX MATCH SO TH ARRAY CAN BE PLOTTED
CV($)=[];	// REMOVE THEN n+1 INDEX SO THE ARRAY SIZE AND INDEX MATCH SO TH ARRAY CAN BE PLOTTED

// PLOT THE SIMULATION
// CREATE TO PLOTS USING THE SUBPLOTS FUNCTION
// THE TOP PLOT SHOWS THE SP AND PV. THE BOTTOM PLOT SHOWS THE CV 

subplot(2,1,1);
plot([1:n],[SP PV]);
xtitle('Heat Exchanger Simulation', 'Time In Seconds', 'Temperature');
legend("SP","PV");
subplot(2,1,2);
plot([1:n],[CV]);
xtitle('Heat Exchanger Simulation', 'Time In Seconds', 'Control Output%');
legend("CV")
// END

This is what the output looks like

I am just learning Scilab too so I don't have the pretty formating down yet.

Pandiani · Sep 13, 2006

It would be good that anyone more experienced confirms if my thoughts about Siemens FB41 and it's DISV parameter are correct.
Also, every suggestion is more than welcome.

Peter Nachtwey · Sep 13, 2006

First we learn to crawl.

It isn't clear what you want to know about the DISV. It is a bias that is added on to the output. First you need to understand the control guru example. Then I can modify the program do use the bias. Have you downloaded Scilab yet?

I should have had lots of questions about the .pdf and Scilab program by now.

When we use the DISV, bias or feed forward we will not use the Ki or the integrator gain.

Pandiani · Sep 14, 2006

Hello, Peter

I have doewnloaded it, but it will take some time to learn it so I can participate. Since we want to crawl first, I guess I'll need to take some time.
When you say that DISV is not used if Ki is used, that is similar to what I wrote earlier. Beacuse integral part is continuous sum, controller output will have constant value on it's output after some time. However that control guru example use CObias + K*e(t)+K/Ti * integral(e(t)). Why Cobias is used here when there is integral part. I need to understand this first.
I guess I would be able to make script files in SCILAB to simulate outputs and system response since matlab only shows dynamic. I'll wait weekend to make my first crawl in scilab.

Peter Nachtwey · Sep 14, 2006

Pandiani, I use feed forwards, bias or DISV all the time in motion control. In temperature control the SP is usually changed in steps. This lowers the usefulness of the bias.

There are some times when a bias could be used. If you have an oven you can increase the bias as a function of cold material going in and heated material going out. What ever is being heated is absorbing energy so the faster. The bias can be calculated on-the-fly to account for the energy removed by the product. That will come later.

You be able to cut and paste that example in to the scilab editor and load it. I tried it so I know it works. I expected about a half dozen questions at least from you and others by now.

Pandiani · Sep 14, 2006

It will be dozen of questions, but not before late Saturday...
P.S. I'm really willing to learn, I start to work in one firm and if I understand and learn these issues it will be exactly what I want to do in my life.
Thanks Peter!

Gerry · Sep 16, 2006

playing along...

I downloaded scilab and pasted in Peter's script and got the same plots, but I don't feel enlightened. So I started tinkering with ControlLogix and made this to model the heat exchanger.

This is the response for a step change in CO from 39% to 42%. I've inserted the bias term Inlet_Temperature to get the Outlet_Temperature into the right range (presumes a linear Kp).

Later, I'll try to progress through the ControlGuru tutorials.

Peter Nachtwey · Sep 16, 2006

Getting the same result is only the beginning.

Try simulating the PID too.

Aren't you curious? What happens if:
1. You change the PID update time?
2. What are the minimum and maximum limits?
3. What happens if the dead time is increased and decreased.
4. What happens if the Tc ( closed loop time constant ) is increased or decreased?
5. Can you control the system with just a P gain?
6. Can you control the system with a bias and a P gain?
7. Do you believe the formulas on the control guru site work if you know the plant gain, time constant and dead time?
8. How do you find the plant gain, time constant and dead time?
9. What good is the derivative gain?
10. The simulation uses double precision floatiing point. What if the temperature feedback and the output have limited resolution so that the input and output are using 12 bit AtoD and DtoA converters?

Lets assume Kp is linear for now.

Pandiani · Sep 16, 2006

Okay, I have read your script but I must say I'm pretty much confused.
First you decide to use discrete model of this control system. It's logical, because only discrete systems can be simulated the way you did it, with arrays. I used not to do anything that I don't understand, and this is what I don't uderstand.
Differential equation is:
Tp*DPV/dt+PV = CO(t-theta)
and tansfer function is: Kp*exp(-s+theta)/(Tps+1). How did you get your model:
PV(n+1)=A*PV

+B*CVd

+C;
where A = exp(-T/Tp), B = (1-A) and C = (1-A)*PVss?
I have searched my theory books and tried to obtain Z transform, but that didn't help.
Also I tried to use state space model of system with dead time and I got difference equation that was similar but not exacltly like yours. So I'm pretty much stuck at the begining. I'm afraid we're getting far from what we want. I wanted to make connections between theory and practice. If I just play around with model which I don't understand I'll learn nothing.
Please can you post link or explain how you obtained this equation?
Thanks

Pandiani · Sep 16, 2006

Here in two pdfs you can find mathematical procedure that I tried to get your equation. This is similar procedure I ave fond in my book and also result is the same to those I found on the Internet. The main problem is that C parameter. You said that is tricky part and you need it because of PVss. However I cannot find mathematical justification for this parameter. I got your A and B parameters by applying procedure you can find in pdfs.
If I understood correctly editor of controlguru is professor with phD and it would be nice to hear his opinion if we can include his somehow.
The problem and confusion is pretty much same I have with that CObias, because whole math theory in books (Modern Control Systems by Dorf and Bishop and in many control books on my language simply don't include it).
In these pdfs I marked PV as 'y' and CO as 'u' to control myself easier.
Method I used and which is used in many books involve initial condition, but result difference equation doesn't have that C = 1-A parameter. It's easy for me to understand why it is needed when somebody gives you equation. The main problem is how to obtain that equation. I tried many methods and every failed. This state space approach use initial condition and equation is not the same you got it. How did you get that famous PV(n+1) equation???
Please, let's solve this together, I'm getting a headache...

Peter Nachtwey · Sep 17, 2006

Down the rabbit hole

I have attached a .pdf file the shows TWO step by step methods of calculating the difference equations for the heat exchanger.

Note, I am self taught. I had to learn control theory because I must compete with the best in the world. Don't ever think the waht the professors taught you is all there is.

Matlab and Simulink are good for getting answers if you know how to ask the questions. Mathcad make you find the answers the hard way but then you know how the answers are calculated.

Pandiani · Sep 17, 2006

Let's go one step at a time

Thanks Peter for your time and effort.

What you have done is exactly the same to what I have done. You use Z-transform approach to get A and B constants in your PV equation using ZOH method or method of step invariance how it is also called. (It's called step invariance vecause you matching step response of continuos system to response of discrete system in sampling points). I used state space approach and we got exactly the same equation (update time is 1 sec rather then 1 minute). That's good, very good. Also in my work you can find proof for Ac and Bc using state space.

You have done just like me, first you have found equation by assuming there is no time delay.

First thing that is little unclear is CO when introducing delay. In your article from controlguru, you wrote CO(n-thetap) (in the begining, but later use n*thethap), and in this pdf you wrote thetap/T. Since it is discrete equation I agree there should be CO(n-trunc(thetap/T)). I do understand that this is only approximation with trunc function that has no big influence because theta is much larger than T. Theory in book about discrete system with dead time is full of advanced math, but we can leave it this way.

So, until part with headline „Steady State Compensation. THIS IS NOT IN THE TEXT BOOKS!!!!!!“ everything is more or less same like I did in my pdfs. You described it and it makes perfect sense why is needed to have that C part, so I'm accepting this. Unfortunately many of real time things aren't described in textbooks, I'm living proof, I have very good grades, but it seems I'm going too deep into theory trying to find explanations with constraints I found in textbooks and simply have no sense what some formula means in real life.

I'm becoming frustrated because, althoug I make good controller (PID, or state feedback) in Matlab, when specifications are stated, for example, design controller for system described with G(s) so that output (step response) has steady state error zero, overshoot <10%, settling time <3s, when I face real time problem, I'm not able do solve it on my own.

However I need formal method to apply to wide range of processes and I found that www.controlguru.com website, where is shown how to obtain model (FOPTD) from real process and how to calculate PID gains.

Now I need to fully understand everything that is going on.

PVss=140-39*Kp // STEADY STATE TEMPERATURE

You wrote this:

From the heat exchanger data plot, we can see that the heat exchanger temperature (PV) is 140 °C degrees when the controller output (CO) is 39%. If these values are plugged into the formula y=mx+b, we can solve for the steady state temperature, PVss:

140 = 39Kp + PVss

I don't understand how. If Kp is process gain, then you assumed linear relationship:

PV = Kp*CO+PVss and since 39% of CO will give 140 deg you got this relationship:

140 = 39Kp + PVss. Since Kp is process gain and it's Kp = -0.53, this means that PVss will be 160.7 deg???? What is conclusion of this? Temperature is system's output, so steady state value means that system will have this value after some specific amount of time. So one can conclude that when valve is opened 39%, after some time, output temperature will have have value of 160 deg. But that doesn't happen. Steady state temperature is 140 not 160 Temperature 160 will be when CO goes to 0%. Maybe steady state is not the right therm, it's better to call it offset or something.
Also, you assumed linear relationship y = mx+b? I understand this is an approximation and because of nature of system it's likely that PVss will be little different than 160 deg when CO is 0%. Relationship is more exponential than linear but I assume thi gives good results in practice, right?
Please, explain this!

Advanced control

Pandiani

Lifetime Supporting Member

Peter Nachtwey

Member

Pandiani

Lifetime Supporting Member

Peter Nachtwey

Member

Pandiani

Lifetime Supporting Member

Peter Nachtwey

Member

Pandiani

Lifetime Supporting Member

Peter Nachtwey

Member

Pandiani

Lifetime Supporting Member

Gerry

Member

Peter Nachtwey

Member

Pandiani

Lifetime Supporting Member

Pandiani

Lifetime Supporting Member

Peter Nachtwey

Member

Pandiani

Lifetime Supporting Member

Similar Topics