How to maintain tank water level in Matlab

Can you please help me on this, my brother asked on how to solve it but I am struggling. Thanks in advance.

The question is in the image attached.

[William Rose]

Show the work you have done so far on this problem. What is your answer to question 1, the material balance equations for the system? You will need that answer before you write Matlab code. What are the dimensions for Alpha1, alpha2, alpha3, and what do those alphas represent? (Presumably hydraulic resistance or conductance, but which?) Please include the code you have written so far, to attempt this problem.

[Anass Naoushy]

close all

clear all

clc

%% Our Material Balances:

% Around Tank One: dh1/dt = (F1 - F2 - F5) / A1;

% Around Tank Two: dh2/dt = (F4 + F2 - F3) / A2;

% Around The System: dH/dt = ??

%% Defining The Variables:

syms F1 F2 F3 F4 F5 h1 h2 A1 A2 a2 a3 a5;

%% Defining The Material Balance Equations:

F2 = h1/a2;

F5 = h1/a5;

F3 = h2/a3;

f1 = (F1 - F2 - F5) / A1;

f2 = (F4 + F2 - F3) / A2;

%% Performing The Jacobians:

A = jacobian([f1 f2], [F1 F4]);

B = jacobian([f1 f2], [h1 h2]);

g1 = h1;

g2 = h2;

C = jacobian([g1 g2], [F1 F4]);

D = jacobian([g1 g2], [h1 h2]);

%% Defining The Values:

F1 = 0.35;

F4 = 0.215;

A1 = 2;

A2 = 1.25;

a2 = 0.15;

a5 = 0.22;

a3 = 0.18;

h1ss = 2.0228;

h2ss = 1.9616;

h1 = h1ss;

h2 = h2ss;

%% Evaluate The Jacobians:

AA = eval(A);

BB = eval(B);

CC = eval(C);

DD = eval(D);

%% Convert To Space System:

[q1, p1] = ss2tf(AA, BB, CC, DD, 1);

[q2, p2] = ss2tf(AA, BB, CC, DD, 2);

%% Find The Transfer Functions:

G1 = tf(q1(1, :), p1); % h1/F1

G2 = tf(q1(2, :), p1); % h1/F4

G3 = tf(q2(1, :), p2); % h2/F1

G4 = tf(q2(2, :), p2);% h2/F4

%% Plotting The Figures:

time = 0:1:50;

response = zeros(size(time));

response(10: 15) = 3;

simResponse = lsim(G1, response, time);

finalResponse = simResponse + h1ss;

figure;

plot(time, finalResponse);

xlabel('Time (s)');

ylabel('Height Tank One');

title('Rectangular Response of h1 wrt Disturbing F1');

responseTwo = zeros(size(time));

responseTwo(10: 15) = 3;

simResponseTwo = lsim(G3, responseTwo, time);

finalResponseTwo = simResponseTwo + h2ss;

figure;

plot(time, finalResponseTwo);

xlabel('Time (s)');

ylabel('Height Tank Two');

title('Rectangular Response of h2 wrt Disturbing F1');

%%

responseThree = zeros(size(time));

responsThree(17: 24) = 1.75;

simResponseThree = lsim(G2, responseThree, time);

finalResponseThree = simResponseThree + h1ss;

figure;

plot(time, finalResponseThree);

xlabel('Time (s)');

ylabel('Height Tank One');

title('Rectangular Response of h1 wrt Disturbing F4');

responseFour = zeros(size(time));

responseFour(17: 24) = 1.75;

simResponseFour = lsim(G4, responseFour, time);

finalResponseFour = simResponseFour + h2ss;

figure;

plot(time, finalResponseFour);

xlabel('Time (s)');

ylabel('Height Tank Two');

title('Rectangular Response of h2 wrt Disturbing F4');

This is my attmpt to solve the problem, but I am sure I have mistakes, a lot.

[per isakson]

See How can I do a Simulink 3D Animation for a system which deals with the regulation of the water' level in a tank?

[Anass Naoushy]

But I need the solution as code, not as a simulation unfortunately.

[William Rose]

@Anass Naoushy,

I see you have converted "alpha1" to a2=.15, and "alpha2" to a5=.22, and "alpha3" to a3=.18. I like your subscripts better than the subscripts in the problem, because yours match the flows.

I think your equation for F2 is incorrect: it should be F2=(h1-h2)/a2.

I think a third material balance equation for the system is obtained by considering total inflow minus total outflow:

.

That is about all I will be able to contribute to this discussion, since I am not familiar with symbolic math in Matlab.

[Anass Naoushy]

@William RoseThanks a lot. I will try them now.

[William Rose]

@Anass Naoushy,

The third material balance equation which I proposed is not independent. It can be obtained as a linear combination of the two material balance equations which you identified.

[Anass Naoushy]

@William Rose

Yes, I got that.

One thing I am stuck at is that I keep obtaining a straight plot when perturbating F1 and checking the behavior of height 2. The same happens for F4 with height 2. Any ideas?

[William Rose]

@Anass Naoushy,

I ran your code and it runs without any errors. I see that figures 2 and 3 show no change in tank heights. I don;t know why. But I must admit I do not understand the code, because I am not familiar with ss2tf() or with Jacobian() or with symbolic math in Matlab.

I also see that the responses of heights to flow changes in figures 1 and 4 are immediate, and do not show signs of exponential build-up and decay as I woudl expect. Therefore I am not confident about those repsonses.

[William Rose]

I simulated the system using the differential equation solver. I notice that the time constant of decay after changes is a lot shorter than 1, which means a time step size of 1 will be too long to show the exponential decay after perturbations. That could explain why you do not see exponential decays in your results, but it does not explain why yyou do not see any change in response to certain perturbations. The so-called steady state values for h1 and h2 do not seem to be consistent with the default flow rates into the system. This is evident from the rapid decline of both heights from their initial "steady state" values. See attached plots. I run the simulation to t=20 sec, since nothing happens after t=20. This is a simulation of the F1 perturbation (problem 2a). Code atached.

Answer 1

[William Rose]

@Anass Naoushy,

The figures show simulaitons of responses to step changes in F1 and in F4. The attached code is for the seocnd simulaiton: step in F4 (problem 2C). To simulate the step in F1 (problem 2A), adjust the lines in the functions Flow1() and Flow4(), to get the step in F1 and to remove the step from F4. Adjust the plot title to say "F1" or "F4".

I realize this is not your assignment, because you are doing this symbolically. This shows a numerical solution. See my comments separately about the "steady state" heights, which are not really steady state.