TMA681 Tillämpad matematik K, 2001

STUDIO 5. THE TANK REACTOR: DESIGN FOR STABILITY

Important: write a readable report of your work in each studio session. You will need this at the end of the course when you prepare for the written exam; several exam questions will be based directly on the studio work which is only documented in your own notes and the instructions.

Some exercises are called homework in order to save time in the studio classroom.

Theory: S. Larsson, “Kompletterande föreläsningsanteckningar och övningar” (K).

1. Introduction

Begin by quickly repeating the exercises in Studio 4.

Recall that we determined the control variables Ū = [Ū_{1} ; Ū_{2}] so that the reactor would operate

at (c_{f} − c̄)/c_{f} = 0.5 (“50% omsättningsgrad”) and at reactor temperature T̄ = 99^{∘}C, i.e., at

X̄ = [X̄_{1} ; X̄_{2}] = [\begin{matrix}0.5 \\ (99 + 273.15)/T_{f} \end{matrix}].

We also found that (with the present choice of parameters) the desired operating point, X̄, is unstable with respect to perturbations of the initial value. This means that a small deviation from X = X̄ causes the tank reactor to depart from the desired operating point. Since these kinds of perturbations are inevitable in practice, the reactor does not remain in the desired state, which is therefore not stable. Rather, it (depending on the initial perturbation) approaches one of two other equilibrium points, which seem to be stable ones. These two are also stationary points, corresponding to the same control variables Ū.

We also performed a linear stability analysis based on the linearized equation

(1)  \begin{matrix}x′(s) = Ax(s) + Bu(s), s > 0, \\ x(0) = x_{0}, \end{matrix}

for the approximate perturbation x(s) ≈ ΔX(s) caused by the perturbations in input data x_{0} = ΔX_{0} and u(s) = ΔU(s). In (1),

(2)  A = [\begin{matrix}\frac{∂F_{1}}{∂X_{1}} (X̄, Ū) & \frac{∂F_{1}}{∂X_{2}} (X̄, Ū) \\ \frac{∂F_{2}}{∂X_{1}} (X̄, Ū) & \frac{∂F_{2}}{∂X_{2}} (X̄, Ū) \end{matrix}] = [\begin{matrix}−Ū_{1} − f(X̄_{2}) & −X̄_{1}f′(X̄_{2}) \\ αf(X̄_{2}) & −Ū_{1} + αX̄_{1}f′(X̄_{2}) − β \end{matrix}] ,

and

(3)  B = [\begin{matrix}\frac{∂F_{1}}{∂U_{1}} (X̄, Ū) & \frac{∂F_{1}}{∂U_{2}} (X̄, Ū) \\ \frac{∂F_{2}}{∂U_{1}} (X̄, Ū) & \frac{∂F_{2}}{∂U_{2}} (X̄, Ū) \end{matrix}] = [\begin{matrix}1 − X̄_{1} & 0 \\ 1 − X̄_{2} & β \end{matrix}] ,

are called Jacobi matrices of F (X, U) = [\begin{matrix}F_{1}(X, U) \\ F_{2}(X, U) \end{matrix}] at X̄, Ū.

We wrote a Matlab function for computing A and its eigenvalues. We found that it has two real eigenvalues, one positive and one negative. Because of the positive eigenvalue, one eigenmode in x(s) grows exponentially with time, and this explains the instability of X̄.

Our goal is now to re-design the tank reactor so that the desired operating point, X̄, becomes stable. We will do this by increasing the area A_{K} of the cooler or, equivalently, the parameter

β = \frac{κA_{K}}{ρc_{p}q_{ref}} ,

Date: April 27, 2001, Niklas Ericsson, Stig Larsson, Mathematics, Chalmers University of Technology.

until both eigenvalues of the Jacobi matrix A have negative real parts.

2. A parameter study

In the following exercise we make a parameter study of the eigenvalues of A as functions of the parameter β. Exercise 1. The following script file generates an “interval” of a thousand beta values and computes the eigenvalues of A for each β in the interval.

betamin = 0.1;  % left endpoint of "beta interval"
betamax = 10;   % right endpoint of "beta interval"
betavalues = betamin:(betamax - betamin)/1000:betamax;
N = length(betavalues);  % number of beta values in "beta interval"
lambda1 = zeros(1,N);    % initialize vectors containing eigenvalues,
lambda2 = zeros(1,N);    % each column corresponds to one beta value
for index = 1:N
  beta = betavalues(index); % current beta value
  A = jacobianA(Xbar);      % compute A; note that we do not need to
                            % recompute Ubar since A only depends on
                            % Ubar(1), which does not depend on beta
  eigenvalues = eig(A);     % compute eigenvalues of A
  lambda1(index) = eigenvalues(1);  % "first" eigenvalue of A
  lambda2(index) = eigenvalues(2);  % "second" eigenvalue of A
end

Add the following plot commands to the file and you will see two diagrams showing the real and imaginary parts of the eigenvalues.

clf   % clear current figure
subplot(2,1,1) % breaks the figure into a 2-by-1 matrix; selects top half
hold on
grid  % adds grid lines
plot(betavalues, real(lambda1), ':')  % dotted line
plot(betavalues, real(lambda2), '--') % dashed line
xlabel('beta')
ylabel('Real parts')
title('Dotted: \lambda_{1}   Dashed: \lambda_{2}')
hold off
subplot(2,1,2) % selects bottom half
hold on
grid
plot(betavalues, imag(lambda1), ':') % dotted line
plot(betavalues, imag(lambda2), '--') % dashed line
xlabel('beta')
ylabel('Imaginary parts')
hold off

Run the program first with a rather long β interval, e.g., betamin=0.1, betamax=10. Then zoom in on the interesting parts of the graphs by using shorter intervals. Make sure you can identify the following cases for increasing values of β:

• two real eigenvalues: one positive and one negative;

• two real negative eigenvalues;

• a conjugated pair of complex eigenvalues with negative real parts;

• two real negative eigenvalues. Select a value of β and repeat Exercise 4 of Studio 4. Do this with one β value from each case.

Since Ū_{2} depends on β, you must copy the line Ubar(2)= ... from the file data.m to the top of the file solve2.m. Then you can run the following:

>> data
>> S = 20;
>> X0 = Xbar + [0; 0.05];
>> beta=??
>> solve2
>> beta=??
>> solve2  % and so on

For which values of β is the operating point X̄ stable? Determine a stability threshold β_{0} such that X̄ is stable for β > β_{0}.

3. Stability with respect to perturbation of the control variables

So far we have only studied stability with respect to perturbations x_{0} of initial data X_{0}. Now we discuss stability with respect to perturbations u(s) of the control variables U(s). Exercise 2. Set X0=Xbar. Introduce small fluctuations in the control variables by making the following changes in the file tank2.m:

  U(1) = Ubar(1)+0.01*sin(3*s);
  U(2) = Ubar(2)+0.01*cos(s);

Select a value of β < β_{0} and repeat Exercise 4 of Studio 4. Is X̄ stable?

Select a value of β > β_{0} and repeat Exercise 4 of Studio 4. Is X̄ stable? Homework 1. (K 1.4, K 1.8) Solve the linearized equation (1) analytically. Assume that the matrix A is diagonalizable with eigenvalues λ_{1}, λ_{2} and a basis of normalized eigenvectors g_{1}, g_{2}. We assume that both eigenvalues have negative real part and that the eigenvectors are not nearly parallel, see (K 1.7).

Write x(s), x_{0} and Bu(s) in the eigenvector basis:

(4)  x(s) = y_{1}(s)g_{1} + y_{2}(s)g_{2}, x_{0} = c_{1}g_{1} + c_{2}g_{2}, Bu(s) = f_{1}(s)g_{1} + f_{2}(s)g_{2}.

Insert this into (1) to get

(5)  y′_{k}(s) = λ_{k}y_{k}(s) + f_{k}(s), s > 0; y_{k}(0) = c_{k}; k = 1, 2.

This equation is of the form (5) in Studio 1. The solution is

(6)  y_{k}(s) = c_{k}e^{λ_{k}s} + ∫_{0}^{s} e^{λ_{k}(s−t)}f_{k}(t) dt,

see Studio 1: (6) and Homework 2. To study stability with respect perturbations of the control variables we set x_{0} = 0, so that c_{k} = 0, and assume that u(s) is small, which implies that f_{k}(s) is small too:

(7)  max_{s≥0} |f_{k}(s)| ≤ C max_{s≥0} ∥u(s)∥,

where the constant C depends on the Jacobi matrix B and the angle between the eigenvectors, see Remark 1. Since both the eigenvalues have negative real part, we have |e^{λ_{k}s}| = |e^{(α_{k}+jω_{k})s}| = e^{α_{k}s} ≤ 1. For the absolute value of y_{k}(s) we then get

|y_{k}(s)| = | ∫_{0}^{s} e^{λ_{k}(s−t)}f_{k}(t) dt| ≤ ∫_{0}^{s} |e^{λ_{k}(s−t)}| |f_{k}(t)| dt

= ∫_{0}^{s} e^{α_{k}(s−t)}|f_{k}(t)| dt ≤ ∫_{0}^{s} e^{α_{k}(s−t)} dt max_{t≥0} |f_{k}(t)|

= \frac{−1}{α_{k}}(1 − e^{α_{k}s}) max_{t≥0} |f_{k}(t)| (because α_{k} < 0)

≤ \frac{−1}{α_{k}} C max_{t≥0} ∥u(t)∥.

This leads to

(8)  ∥x(s)∥ ≤ |y_{1}(s)| + |y_{2}(s)| ≤ 2\frac{1}{min |α_{k}|} C max_{t≥0} ∥u(t)∥.

Hence, x(s) stays small if u(s) is small. We conclude that X̄ is stable with respect to perturbation of U. Remark 1. The previous stability analysis was made under the assumption that both eigenvalues have negative real part, α_{k} < 0. However, the inequality (8) is rather useless if min |α_{k}| ≈ 0, i.e., if one of the real parts is close to 0. Therefore, in order to secure stable behavior we must not only make sure that both of the α_{k} are negative, but also that they are not too close to zero.

This is similar to what happens if the eigenvectors are nearly parallel, which we also need to avoid. In order to understand this we look at the inequality (7) in more detail. Let θ be the angle between the normalized eigenvectors g_{1} and g_{2} and assume that θ ≈ 0. Then (assuming for simplicity that the eigenvectors and the f_{k} are real; this is the case when the eigenvalues are real)

∥Bu(s)∥^{2} = f_{1}(s)^{2} + f_{2}(s)^{2} + 2f_{1}(s)f_{2}(s)(g_{1}, g_{2})

= f_{1}(s)^{2} + f_{2}(s)^{2} + 2f_{1}(s)f_{2}(s) cos(θ)

≥ (1 − cos(θ))(f_{1}(s)^{2} + f_{2}(s)^{2}) (because −2f_{1}f_{2} ≤ f_{1}^{2} + f_{2}^{2} and cos(θ) > 0)

≈ \frac{1}{2}θ^{2}(f_{1}(s)^{2} + f_{2}(s)^{2}) ≥ \frac{1}{2}θ^{2}f_{k}(s)^{2}

so that (approximately, when θ is small)

|f_{k}(s)| ≤ \frac{\sqrt{2}}{θ}∥Bu(s)∥ ≤ \frac{C_{1}}{θ} ∥u(s)∥.

We see that the constant in the inequality (7) becomes very large when the eigenvectors are nearly parallel.

4. Conclusion

We finally complete the design of a stable tank reactor. Exercise 3. Choose a value of β which is safely bigger than the stability threshold β_{0} (see Remark 1). Then compute the corresponding value of the cooler area A_{K}, the cooler temperature T_{K}, and the cooling power (how much heat that passes through the cooler per second).