Temperature-controlled 3 tank system

sharing on inlet flow using on-off logical control modelled with the hybrid dmss class.
Open Live Script in Matlab
clear, close all
Example for CoDIT 2024 Conference [1].
Table of Contents

Model creation

The system is illustrated in figure 1. It consist of 3 water tanks sharing one inlet flow u2. The water tanks transfer heat between each other, as well as losing heat to the environment with the ambient temperature u1. If multiple valves (z1, z3 and z5) are opened at the same time the inlet flow will be shared uniformly between those. The control is designed to maintain the temperatures of the tanks (x1, x2 and x3) between lower (u3, u5 and u7) and upper bounds (u4, u6 and u8). There is one 2-Point controller per valve/tank.
Figure 1: Illustration of the temperature-controlled 3 tank system.
The system can be described by the following difference equations and inequality constraints:
eq = {"m*cp*(dx1-x1)/ts = k0A0*(u1-x1) + kinAin*(x2-x1) + kinAin*(x3-x1) + z1/y1*(u2-x1)*cp*fv";...
"m*cp*(dx2-x2)/ts = k0A0*(u1-x2) + kinAin*(x1-x2) + kinAin*(x3-x2) + z3/y1*(u2-x2)*cp*fv";...
"m*cp*(dx3-x3)/ts = k0A0*(u1-x3) + kinAin*(x1-x3) + kinAin*(x2-x3) + z5/y1*(u2-x3)*cp*fv";...
"y1 = z1+z3+z5+(1-z1)*(1-z3)*(1-z5)";...
"0 >= (u3-x1)*(1-z1+z2)";...
"0 >= (u5-x2)*(1-z3+z4)";...
"0 >= (u7-x3)*(1-z5+z6)";...
"0 >= (x1-u4)*(1-z2+z1)";...
"0 >= (x2-u6)*(1-z4+z3)";...
"0 >= (x3-u8)*(1-z6+z5)"};
where the first three equations represent the dynamic energy balance of the three water tanks, the fourth equation calculates the share of water flow each tank with an active valve is getting and the inequlities represent the activation of the upper and lower bounds.
The set of equations and inequalities can be transformed into a dmss model by
ts = 0.1;
hybrid3TankSys = sym2dmss(eq, ts, ["dx", "x", "u", "y", "z"]);
Reduced by 26 column(s) and 0 equation(s) with trivial Reduction due to duplications.
with a time step size of 0.1 seconds. The set of equations and inequalities can be viewed:
disp(hybrid3TankSys.symbolicEquations)
To simulate the system the symbolic parameters must be replaced by numeric values
m = 5;
cp = 4200;
k0A0 = 250* (0.5*1);
kinAin = 250* (0.5*1);
fv = 0.09;
 
symbolicParameters = [sym("m") sym("cp") sym("fv") sym("ts") sym("k0A0") sym("kinAin")];
numericParameters = [m cp fv ts k0A0 kinAin];
hybrid3TankSys = hybrid3TankSys.replaceSymbolicParameters(symbolicParameters, numericParameters);

Simulation

For simulation the initial states and input signals have to be defined
x0 = [35, 30, 40];
t = 0:ts:1000;
u = [0, 80, 30, 40, 25, 35, 35, 45].*ones([length(t), hybrid3TankSys.m]);
which are held constant for this example.

Sparsity pattern analysis

Simulation can be performed by using the dmsim class, with by default using a equation reordering algorithm and identification of sub-problems and explicit solvable equations by sparsity pattern analysis.
The sparsity pattern can be calculated and viewed by the following procedure (also automaticly computed inside the simulation routine).
The sparsity pattern is calculated by
I = hybrid3TankSys.incidenceMatrix;
I = [I.equality.stateDerivative, I.equality.algebraic, I.equality.boolean;...
I.inequality.stateDerivative, I.inequality.algebraic, I.inequality.boolean];
f1 = figure();
f1.Position(3:4) = [140 105];
spy(I)
and resorted using the Dulmage-Mendelsohn decomposition into a block lower triangular representation:
I = hybrid3TankSys.orderIncidenceMatrix(I);
f2 = figure();
f2.Position(3:4) = [140 105];
spy(I)

Simulation with and without sparsity pattern analysis algorithm

Simulation time can be compared with and without using the sparsity pattern analysis algorithm:
tic
result = dmsim(hybrid3TankSys, x0, t, u);
time = toc()
time = 6.6791
 
opt.equationReordering = false;
tic
resultOff = dmsim(hybrid3TankSys, x0, t, u, opt);
timeOff = toc()
timeOff = 58.5345
The simulation times per time step are:
disp(time/length(t))
6.6784e-04
disp(timeOff/length(t))
0.0059

Simulation results

The simulated temperatures can be accessed via the state property of the simulation object and can be plotted by
figure()
plot(result.tsim, result.x)
hold on
ylabel('temperatures [°C]')
xlabel('time [sec]')
hold off
as well as the binary valve signals (with an according time signal)
figure()
plot(result.zt, result.z(:,[1,3,5]))
hold on
ylabel('valve position [-]')
xlabel('time [sec]')
hold off
With the simulated state signals the phase diagram can be plotted, revealing a period-2 orbit
figure()
plot3(result.x(:,1), result.x(:,2), result.x(:,3))
grid on
xlabel("x1")
ylabel("x2")
zlabel("x3")
view([-67.5,29.5])

References

[1] T. Warnecke and G. Lichtenberg, "Hybrid implicit multilinear simulation using difference algebraic equations reordering by sparsity patterns", 2024 10th International Conference on Control, Decision and Information Technologies (CoDIT), Valletta, Malta, 2024. In press.
Author(s): Torben Warnecke