PID problem...

[plcking]

Have you tried kicking it? Not too hard but a delicate kick.

 

[Pandiani]

Since there were not many replies to this thread for a long time, I hope we can conclude few things.
First of all, can you agree with nonlinear model and differential equations in my previous post?
Second, do you agree with the way I linearized the model? I applied Taylor approximation consistently.
Should we proceed further to the LQR controller?

 

[Peter Nachtwey]

Go ahead

Since there were not many replies to this thread for a long time, I hope we can conclude few things.
First of all, can you agree with nonlinear model and differential equations in my previous post?
Second, do you agree with the way I linearized the model? I applied Taylor approximation consistently.
Should we proceed further to the LQR controller?

The B,C and D arrays are good. I haven't had time to look at the A array.
Your plots look good so you can't be too far off. I notice you ignored the temperature loss. I will continue with using the percentage trick mentioned above.

I have been gone for a while with out my Mathcad. I should get back this afternoon and then I can contribute more. I was at the OTC ( Off shore technology conference ) in Houston, visiting and training customers.

 

[Pandiani]

Hello guys,
it's time to move on.
Here is SCILAB code that calculates LQR gains:

clear
clc
T1 = 65;
T2= 20;
T0 = 40;

F10 = 0.0042*0.4;
F20 = 0.0042*0.5;
L0 = 1;
Area = 1;
Ac = [-(F10+F20)/(Area*L0), (F10*(T0-T1)+F20*(T0-T2))/(Area*L0^2);
	0 0];

Bc = [(T1-T0)/(Area*L0), (T2-T0)/(Area*L0); 
	1/Area, 1/Area];
Cc = [1 0; 0 1];
Dc = [0, 0; 0, 0];


T = 0.001; //Sampling time
Time = 1;

//Discrete state-space representation
I = eye (2,2); //Identity matrix
Ad = I;
for i = 1:5
	Ad = Ad + ((Ac*T)^i)/prod(1:i);
end

B = I;
for i = 1:5
	B = B + ((Ac*T)^i)/prod(1:i+1);
end

Bd = B*Bc*T;

//Cost matrix
Q = [1/T,0; 0, T];
R = [0.1, 0; 0, 0.1];

N1 = floor(Time/T);
P = Q;
 
for i = 1:N1
  K = inv( Bd'*P*Bd+R)*Bd'*P*Ad;  
  P = Ad'*P*(Ad-Bd*K)+Q;  
end

This yields to the following result:
K =
22.385663 0.0087231
- 17.908618 0.0109037

Control is basically calculated: u = -K*x, and x is state vector. Now before we continue with response and simulation I'd like to discuss how matrices Q and R can be chosen.
I tried to execute default Matlab code for (lqr function)
and found out that both matrices Q and R must be positive definitive, so it means there must be costs (positive number) involved in diagonal numbers. Not really sure how to intepret this.
I asked Peter about this and I hope he'll reply soon...

 

[Peter Nachtwey]

I think your Q array is messed up

//Cost matrix
Q = [1/T,0; 0, T];
R = [0.1, 0; 0, 0.1];

You shouldn't have T in the Q array. I now I did on sci.engr.control but that was because I was minimizing error for position, velocity and acceleration and the relationship between them is 1/T. In this case there is no time relationship between temperature and level.

I would
Q = [1,0; 0,1];
R = [0.1, 0; 0, 0.0];
The upper left entry of Q is arbitrary but all other values must be scaled from that. In this case I made the cost of level error very small. An 1 degree error in temperature has the same cost as a 1 meter error in level. I think the it would be better if
Q = [1,0; 0,100];
Now a 1 degree error costs as much a 0.1 meter error. Remember the error is squared.
The R arrays is the cost of flow. I eliminated the cost of cold water flow. I think better would be
R = [0.01, 0; 0, 0.0];
So the cost of hot water doesn't affect the control too much.
Try that.
Also, Matlab has a C2D ( continuous to discrete ) function. Use it. I made an equivalent function in Scilab. I will verify your calculations tonight.

I have been busy lately and this is my first weekend I have had in a while. I haven't been home to work on this problem but I should have time tonight. Meanwhile, it is time to enjoy the sun.

 

[Pandiani]

I would
Q = [1,0; 0,1];
R = [0.1, 0; 0, 0.0];

I have tried before C2d function and use dlqr (discrete lqr), they return more or less same results. It is interesting that Matlab refuse to calculate K matrix when one eigenvalue is 0, therefore R = [0.1, 0; 0, 0.0] will not work in Matlab. I really wonder why...

 

[Lancie1]

How to Make a Simple Problem so Complicated you Must Use Experienced Guy.

So make an indirect system.
Use the hot water direct to the dorm and mix it locally with simple mixing valves. This way everybody can set their own temp.

Another method is to simply put a water mixing valve in the design in place of where the mixing tank is now. The mixing valve should have a temperature feedback from the Hot Output Pipe to Dormintory. If that temperature is too hot, it automatically reduces the flow of hot input water. If the output temperature is too cold, it increases the flow of hot input water, keeping the dorm hot water supply at a constant temperature. Isn't hot water at a constant temperature the goal here?

I think the solution is easy, but it would eliminate all the fun that Peter and Pandiani are having.

On the other hand, it has been a month now, and the students must be getting tired of bathing in that cold water...

30 years ago I installed a solar hot water system in my home. It had a simple temperature controller with two temperature sensors, S1 one on the solar collectors, and S2 on the hot water tank. If S1 > S2+10, then the collector circulation pump ran. If S1 < S2+10, the pump did not run. It worked well. K.I.S.S.

 

[Peter Nachtwey]

You have missed the point.

Another method is to simply put a water mixing valve in the design in place of where the mixing tank is now.

Yes, yes, that has been noted a couple of times above.

This is a simple problem to gain experience controlling a MIMO system. I could have fudged something together in a fraction of the time that would work in this simple example.

I think the solution is easy, but it would eliminate all the fun that Peter and Pandiani are having.

Yes, so you do understand.

Pandiani wants to use the LQR method.
I will get around to using the method Alaric suggested.
We need a fuzzy logic solution to compare.
Why don't you try to solve this problem using fuzzy logic?
Then we can compare methods.

Pandiani and I were going to do this anyway. taz3m just happen to post the problem we were going to solve.

 

[taz3m]

Hi every body, long time no see.......

By the way, it make me please to see that atleast few people are happy with my prob......

Lancie, the fact that you mentioned about the circulation of water through panels.. the Idea with 2 temperature sensor is GREAT...
I had other ideas about dong this task which was more complicated than yours...

Thanks for the solution, its great...

taz......

 

[Pandiani]

Scilab code provided above will return result as 2x2 matrix. This means that control to the system is calculated in the following way:
ut = K11*Temp + K12*Level
ul = K21*Temp + K22*Level

Since there is no integrator in the model we can conclude that there will always be error in response and that we must use PI control. First thing that comes to my mind is to shape control in the following way:
ut = Kit*inegral(errort) + K11*Temp + K12*Level
ul = Kil*inegral(errorl) + K21*Temp + K22*Level

where errort is difference between actual and setpoint temperature while errorl is same but for the level.

In order to use code above we need to add new poles (integrators in transfer function) which will increase dimension of matrices A and B.
Now the real question is how to add new rows in A and B in such way that Scilab script returns useful results, especially because this is MIMO system. Any ideas?

 

[Pandiani]

There is one workaround that can use LQR to have zero steady state.
Since error is e = output - setpoint and sepoint is constant, we can use introduce new state variable.
err = out-ref
derr/dt = dout/dt = C*dx/dt
now we can write:
dXnew/dt = [derr/dt; dz/dt] = [0, C; 0, A]*Xnew + [0;B]*v
Lqr will give the following solution:
v = -K1*err-K2*z
or
u = -K1*integral(err)-K2*x

In this case I made this:

  Q = [ 1,0,0,0; 
  	0,1,0,0; 
  	0,0,1, 0;
  	0,0,0,1];
  
  R=[0.1,0;0,0.001];
  
  
  Aw = [[0,0;0,0;0,0;0,0],[Cc;Ac]];
  Bw = [[0,0; 0,0]; Bc]

and I used Matlab built in function: lqr
[K1,S,E] = lqr(Aw,Bw,Q,R,0);
which returned:
K1 =
0.5426 3.1154 0.5634 3.5321
-31.1539 5.4256 -31.1776 5.9467

This provides me the following model structure:
http://www.plctalk.net/qanda/uploads/Model_lqr.jpg

The following responses are obtained:
http://www.plctalk.net/qanda/uploads/Temp_resp.jpg

http://www.plctalk.net/qanda/uploads/Level_resp.jpg

I'm not satisfied with the response at all, but at least it's stable

 

[Peter Nachtwey]

Scilab code provided above will return result as 2x2 matrix. This means that control to the system is calculated in the following way:
ut = K11*Temp + K12*Level
ul = K21*Temp + K22*Level

Because I used the fractional method I don't have four gains to figure out. Just two but I think we can equate the LQR and the fractional method.

Since there is no integrator in the model we can conclude that there will always be error in response and that we must use PI control. First thing that comes to my mind is to shape control in the following way:
ut = Kit*inegral(errort) + K11*Temp + K12*Level
ul = Kil*inegral(errorl) + K21*Temp + K22*Level

where errort is difference between actual and setpoint temperature while errorl is same but for the level.

In order to use code above we need to add new poles (integrators in transfer function) which will increase dimension of matrices A and B.
Now the real question is how to add new rows in A and B in such way that Scilab script returns useful results, especially because this is MIMO system. Any ideas?

I agree, you end up with a 4x4 array. This greatly complicates things for so little gain. I think the proportional gains can be increased until the desired steady state error is achieved.

Here is my steady state analysis for the fractional gain method.
ftp://ftp.deltacompsys.com/public/NG/Mathcad%20-%20Temp%26Level.pdf

The good stuff is on pages 5/10 and 6/10. At the bottom of page 6 is the two equations for calculating the control output to the hot and cold pumps. It is really very simple. Pages 7-10/10 are equations that didn't wrap to stay in the margins. I duplicated the equations in-line so you need need to bother with pages 7-10.
I also cheated. I used Mathcad's computational abilities to solve for the steady state PVt and PVl. I didn't show the symbolic solution because it stretched across many pagess because no attempt is made to simplify or linearize the process.

I will work on a dynamic response but first I wanted to verify I finally got my model right.

I should also study Rytko's .pdf. I have alway relied on the computational power I have at hand. Without it I would need to resort to the linearization methods.

My conclusion is that a temperature gain of 1 or 0.5 and a flow/level gain of 10 to 20 will work nicely.

I'm not satisfied with the response at all, but at least it's stable

It doesn't look right. You should get nioce smooth first order responses to changes in demand.

 

[Pandiani]

Peter, I'll have a look. We're taking different approaches. First off all do you agree with the way I linearized the system around working points 40°C and 1m (what about matrix A)? I think model is right because of compared responses.

Now, do you agree with the structure of regulator that is shown in last picture (two intgrators and two gains for each input)? I have used Matlab's built in function in order to avoid possible wrong calculations of gains.
I'm not sure why I have that oscillations in the begining. Did you notice that I tried to simulate control valve as real number in interval (0-1)? Is this correct? Output from the controller is limited between 0 and 1 and it practically represent Coh (Coc) which is multiplied with constant pump flow (0.0042).

I'll try to play a little to determine what causes oscillations.
I'll need to fix my own model before I fully focus on yours.

 

[Peter Nachtwey]

Peter, I'll have a look. We're taking different approaches. First off all do you agree with the way I linearized the system around working points 40°C and 1m (what about matrix A)? I think model is right because of compared responses.

I will be able to look at it now that I have my simulation done. I updated my .pdf
ftp://ftp.deltacompsys.com/public/NG/Mathcad%20-%20Temp%26Level.pdf
It works the way I would expect it to work. I did not try to optimize the gains. I don't think it is necessary to put to much effort in the gains as they can be adjusted. within limits, until the desired error is achieved. My simulation changes the demand as a function of time. I made the demand go higher than the available flow for a few minutes to see how the system would respond. You can see the level drops. The temperature increases because when both pumps are fully on the average supplied temperature is (20+65)/2=42.5 degrees. Also, if the demand goes 0 the mixing tank temperature decreases because of heat loss. There needs to be a minimun flow to maintain 40 degrees C. I still have a few minor tweaks to make but I am done with doing it my way.

I'll need to fix my own model before I fully focus on yours.

That is what I have been thinking too which is why I haven't been much help. Now that I have mine done I will try controlling the hot and cold pumps independently like you.

First I will look at the steady state values for your model. They should be similar to mine.

I do have a problem with the block diagrams. I can't see what is in the blocks.

 

[Pandiani]

OK Peter,
In order to make model easier to follow I place mixing tank model inside subsystem.
Here, in attachment you can find model which is related to the following .pdf file http://www.plctalk.net/qanda/attachment.php?attachmentid=8141

I think I have made model correctly. Only calculation of gains and modelling control valves can be problem.
Model is based on the equation presented in .pdf (link above).