The Kronecker Summation Method for Robust Stabilization Applied to a Chemical Reactor

Abstract

The paper focuses on robust stabilization where the suitable parameters of a simple continuous-time PI controller are determined through a combination of the Kronecker summation method, sixteen plant theorem, and an algebraic approach to control design in the ring of proper and stable rational functions. The initial theoretical background is followed by an illustrative experiment which includes computation of the controller and verification of control results for a continuous stirred tank reactor with exothermic reaction modelled as a fourth-order interval system.

1. Introduction

A popular and practically preferred approach to appropriate control synthesis for plants with complex properties consists of construction of an uncertain model, which should cover all possible operating points, parameter variations and nonlinear behavior, and consequent robust control design [1, 2]. Although the resulting controllers can have a simple structure and fixed parameters, which is even acknowledged from the practical application viewpoint, both process modelling and control design are generally nontrivial tasks.

Chemical reactors, which belong among the most interesting and critical processes in all chemical engineering, represent the class of systems suitable for robust control applications. Their control is usually affected by very complex behavior and, moreover, bounded with potential safety problems. A common type of reactor is known as a continuous stirred tank reactor (CSTR) [3]. The mathematical model of CSTR, robustly stabilized in this paper, has been constructed in [4]. Moreover, the same work has presented stabilization of the CSTR using technique from [5, 6] embellished with a polynomial control. Besides, robust static output feedback control has been utilized to this CSTR in [7]. The idea of robust stabilization applied in this paper is similar to [4], but the contribution is mainly in the use of alternative methods, that is, combination of the Kronecker summation method [8] and an algebraic approach to control design under the ring of proper and stable rational functions (R_PS) allowing the elegant tuning [9–14] which has been already investigated in [15, 16].

This paper deals with design of robustly stabilizing continuous-time PI controllers for a continuous stirred tank reactor (CSTR) in which exothermic reaction occurs. The controlled plant is assumed as a fourth-order interval plant, and easily tunable PI controller is designed in order to robustly stabilize the closed control loop. As a synthesis method, the combination of Kronecker sum method, sixteen plant theorem and an algebraic approach, is utilized. This compound and its application to the chemical reactor model represent the key contribution of the work. The robust stabilization is verified on a simulative example.

The paper is organized as follows. In Section 2, the Kronecker summation method is described. Section 3 then extends the idea for an interval plant. Subsequently, Section 4 briefly outlines the algebraic approach to controller design itself. Next, a chemical reactor description and specific control experiment are provided in the extensive Section 5. And finally, Section 6 offers some conclusion remarks.

2. Kronecker Summation Method

An interesting technique for computation of stabilizing PI controllers based on the Kronecker summation has been presented in [8]. The main purpose of the method is to find possible variations of PI controller parameters which ensure stability of the classical one-degree-of-freedom closed control loop according to Figure 1, where𝐺(𝑠)=𝐵(𝑠)𝐴(𝑠) (1) is transfer function of a controlled system with fixed coefficients, and𝐶(𝑠)=𝑘_𝑃+𝑘_𝐼𝑠=𝑘_𝑃𝑠+𝑘_𝐼𝑠 (2) is a PI controller.

Figure 1 

The classical closed control loop.

Remember that Kronecker summation of two general square matrices 𝑄 (with size 𝑞-by-𝑞) and 𝑅 (𝑟-by-𝑟) is defined as [17]𝑄⊕𝑅=𝑄⊗𝐼_𝑟+𝐼_𝑞⊗𝑅, (3) where 𝐼_𝑞,𝐼_𝑟 are identity matrices of size 𝑞-by-𝑞 and 𝑟-by-𝑟, respectively, and where ⊗ stands for the Kronecker product. For example,𝑄⊗𝐼_𝑟=𝑞₁₁𝐼_𝑟⋯𝑞_1𝑞𝐼_𝑟𝑞⋮⋱⋮_𝑞1𝐼_𝑟⋯𝑞_𝑞𝑞𝐼_𝑟. (4) The significant feature of the final square matrix (3) (𝑞𝑟-by-𝑞𝑟) is that it has 𝑞𝑟 eigenvalues which are pair-wise combinatoric summations of the 𝑞 eigenvalues of 𝑄 and 𝑟 eigenvalues of 𝑅. In other words, the operation of Kronecker summation induces the “eigenvalue addition” feature to the matrices. One can exploit this property to obtain the equation for which all pairs (𝑘_𝑃,𝑘_𝐼) leading to purely imaginary roots comply.

The characteristic equation of the closed control loop (Figure 1) has the form𝑃_𝐶𝐿=𝐴(𝑠)𝑠+𝐵(𝑠)𝑘_𝑃𝑠+𝑘_𝐼=𝑓_𝑛𝑘_𝑃,𝑘_𝐼𝑠^𝑛+⋯+𝑓₁𝑘_𝑃,𝑘_𝐼𝑠+𝑓₀𝑘_𝑃,𝑘_𝐼=0. (5) Now define𝑥^′₁=𝑥₂𝑥^′₂=𝑥₃⋮𝑥^′_𝑛=𝑓₀𝑘_𝑃,𝑘_𝐼𝑓_𝑛𝑘_𝑃,𝑘_𝐼𝑥₁−𝑓₁𝑘_𝑃,𝑘_𝐼𝑓_𝑛𝑘_𝑃,𝑘_𝐼𝑥₂𝑓−⋯−^𝑛−1𝑘_𝑃,𝑘_𝐼𝑓_𝑛𝑘_𝑃,𝑘_𝐼𝑥_𝑛, (6) and transform (5) into matrix differential equation𝑋^′=𝑀𝑋, (7) where 𝑀 is matrix of size 𝑛-by-𝑛−𝑓𝑀=0100⋯00010⋯0⋮001⋯⋮⋮⋮⋮⋮⋱0000⋯01₀𝑘_𝑃,𝑘_𝐼𝑓_𝑛𝑘_𝑃,𝑘_𝐼−𝑓₁𝑘_𝑃,𝑘_𝐼𝑓_𝑛𝑘_𝑃,𝑘_𝐼−𝑓₂𝑘_𝑃,𝑘_𝐼𝑓_𝑛𝑘_𝑃,𝑘_𝐼𝑓⋯⋯−^𝑛−1𝑘_𝑃,𝑘_𝐼𝑓_𝑛𝑘_𝑃,𝑘_𝐼, (8)and 𝑋^′=[𝑥^′₁,𝑥^′₂,…,𝑥^′_𝑛]^T, 𝑋=[𝑥₁,𝑥₂,…,𝑥_𝑛]^T. Equations (5) and (7) are connected by relation𝑃_𝐶𝐿=𝑓_𝑛𝑘_𝑃,𝑘_𝐼det(𝑠𝐼−𝑀)=0. (9)

Evidently, the same complex variable 𝑠 is both the root of (5) and the eigenvalue of 𝑀. Thanks to the fact that 𝑀 is a constant matrix, the complex conjugates of 𝑠 must satisfy also (9), that is,det𝑠^∗𝐼−𝑀=0. (10) Consequently, as it was presented in [8], if 𝑠=𝑗𝜔 is the root of (5) it must be the eigenvalue of 𝑀. Furthermore, 𝑠^∗=−𝑗𝜔 is also the root of (5) and the eigenvalue of 𝑀. As the summation of two eigenvalues 𝑠=𝑗𝜔 and 𝑠^∗=−𝑗𝜔 is equal to zero, the Kronecker summation of two matrices must be singular when such correspondence of 𝑘_𝑃, 𝑘_𝐼, and 𝜔 occurs. Thus,det(𝑀⊕𝑀)=0 (11) determines the stability boundary in (𝑘_𝑃,𝑘_𝐼) plane, because every couple (𝑘_𝑃,𝑘_𝐼) satisfying (11) means that the same couple inserted into (5) will lead to the pair of conjugate purely imaginary roots or zero roots. These positions together with the line 𝑘_𝐼=0 are the only ones where the system stability can shift. Generally, the stability boundary splits the (𝑘_𝑃,𝑘_𝐼) plane into the stable and unstable areas. The selection of the stabilizing region(s) can be performed through a test point and corresponding representative polynomial within each area.

3. Robust Stabilization of Interval Systems Using PI Controller

So far, we could apply the Kronecker sum method to calculate a region of stabilizing PI controller parameters for plant with fixed coefficients. However, the papers [5, 6, 8] extended this (or an alternative) stabilization technique also for interval systems. The simple idea consists of the combination with the so-called sixteen plant theorem [1, 18, 19]. This proposition says that a first-order controller (such as PI controller) robustly stabilizes an interval plant𝐺()=𝐵()𝑠,𝑏,𝑎𝑠,𝑏𝐴()=∑𝑠,𝑎^𝑚𝑖=0𝑏⁻_𝑖,𝑏⁺_𝑖𝑠^𝑖𝑠^𝑛+∑^𝑛−1_𝑖=0𝑎⁻_𝑖,𝑎⁺_𝑖𝑠^𝑖,𝑚<𝑛, (12) if and only if it stabilizes its sixteen Kharitonov plants. The values 𝑏⁻_𝑖,𝑏⁺_𝑖,𝑎⁻_𝑖, and 𝑎⁺_𝑖 represents, respectively, lower and upper bounds for parameters in numerator and denominator.

Remember that the Kharitonov plants are defined as𝐺_𝑖,𝑗(𝑠)=𝐵_𝑖(𝑠)𝐴_𝑗(𝑠), (13) where 𝑖,𝑗∈{1,2,3,4}; 𝐵₁(𝑠) to 𝐵₄(𝑠) and 𝐴₁(𝑠) to 𝐴₄(𝑠) are the Kharitonov polynomials for the numerator and denominator of the interval system (12), that is [20],𝐵₁(𝑠)=𝑏⁻₀+𝑏⁻₁𝑠+𝑏⁺₂𝑠²+𝑏⁺₃𝑠³𝐵+⋯₂(𝑠)=𝑏⁺₀+𝑏⁺₁𝑠+𝑏⁻₂𝑠²+𝑏⁻₃𝑠³𝐵+⋯₃(𝑠)=𝑏⁺₀+𝑏⁻₁𝑠+𝑏⁻₂𝑠²+𝑏⁺₃𝑠³𝐵+⋯₄(𝑠)=𝑏⁻₀+𝑏⁺₁𝑠+𝑏⁺₂𝑠²+𝑏⁻₃𝑠³+⋯, (14) and analogically𝐴₁(𝑠)=𝑎⁻₀+𝑎⁻₁𝑠+𝑎⁺₂𝑠²+𝑎⁺₃𝑠³𝐴+⋯₂(𝑠)=𝑎⁺₀+𝑎⁺₁𝑠+𝑎⁻₂𝑠²+𝑎⁻₃𝑠³𝐴+⋯₃(𝑠)=𝑎⁺₀+𝑎⁻₁𝑠+𝑎⁻₂𝑠²+𝑎⁺₃𝑠³𝐴+⋯₄(𝑠)=𝑎⁻₀+𝑎⁺₁𝑠+𝑎⁺₂𝑠²+𝑎⁻₃𝑠³+⋯. (15)

Thus, robust stabilization of interval system directly follows from the simultaneous stabilization of all sixteen fixed Kharitonov plants. The final stability region for the original interval system is given by the intersection of sixteen partial regions obtained individually via the Kronecker summation method.

4. Algebraic Tuning of PI Controller

Recapitulate that for now we are able to compute all possible robustly stabilizing variations of proportional and integral parts in PI compensator. However, the final choice of the controller from the obtained stability region represents another task. A very good survey on PI(D) control issues is provided, for example, in [21]. For the purpose of this paper, a simple but effective solution offers an algebraic approach to control synthesis [9–11]. This fractional design is grounded in general solutions of Diophantine equations in R_PS, Youla-Kucera parameterization of controllers, and conditions of divisibility in R_PS. One of main advantages of this method can be seen in the existence of the only tuning parameter 𝑚>0 which serves for influencing the control behavior. This paper is not going to explain the details of this approach. Interested readers can find them, for example, in [11–14]. This work only takes advantage of one simple tuning rule. The parameters of PI controller (2) can be computed in compliance with𝑘_𝑃=2𝑚−𝑎₀𝑏₀,𝑘_𝐼=𝑚²𝑏₀, (16) where 𝑎₀ and 𝑏₀ come from the first-order nominal controlled plant𝐺_𝑁(𝑠)=𝑏₀𝑠+𝑎₀, (17) and where the tuning parameter 𝑚>0 can be selected, for example, using the recommendation from [12]𝑚=𝑘𝑎₀. (18) Suitable coefficient 𝑘 depends on the size of first overshoot. For example, the choice 𝑘=2.14, which is applied in the following simulation experiments, leads to 3% overshoot.

5. Application to a Chemical Reactor

5.1. CSTR Description

The controlled process adopted from [4, 7, 22] represents hydrolysis of propylene oxide to propylene glycol in a CSTR. More specifically, the chemical reaction of the process isC₃H₆O+H₂O⟶C₃H₈O₂. (19) Besides propylene oxide and water, methanol is also added to the CSTR in order to improve the solubility of propylene oxide in water. The excess of water ensures higher selectivity to propylene glycol and eliminates consecutive reactions of propylene oxide with nascent propylene glycol. Dependence of the rate constant of chemical reaction on the temperature can be described by the well-known Arrhenius equation𝑘=𝑘_∞𝑒^{−𝐸/𝑅𝑇𝑟}, (20) where 𝑘 means reaction rate constant, 𝑘_∞ is the pre-exponential factor, 𝐸 represents the activation energy, 𝑅 signifies the universal gas constant, and 𝑇_𝑟 is the temperature of the reaction mixture.

Under assumption of ideal mixing in the CSTR, constant reacting volume, and the identical volumetric flow rates of the inlet and outlet streams, the mass balance of the system can be given by𝑉_𝑟𝑑𝑐_𝑖𝑑𝑡=𝑞_𝑟𝑐_𝑖0−𝑐_𝑖+𝑉_𝑟𝜈_𝑖𝑟𝑖=1,2,3, (21) where 𝑉_𝑟 stands for the reacting volume, 𝑐_𝑖 means the molar concentration of the 𝑖th component, 𝑐_𝑖0 is the feed molar concentration of the 𝑖th component, 𝑞_𝑟 represents the volumetric flow rate of the reaction mixture, 𝜈_𝑖 determines the stoichiometric coefficient of the 𝑖th component, and 𝑟=𝑘𝑐_C3H6O is the molar rate of the chemical reaction.

Further, independency of the specific heat capacities, densities, and volumetric flow rates on temperature or mixture composition has been supposed. Moreover, the mixing volume and the heat of mixing have been neglected. So, the simplified enthalpy balance of the reaction mixture and the simplified enthalpy balance of the cooling medium introduced in the monograph [23] and subsequently in the papers [4, 7] can be formulated as, respectively,𝑉_𝑟𝜌_𝑟𝑐_𝑝𝑟𝑑𝑇_𝑟𝑑𝑡=𝑞_𝑟𝜌_𝑟𝑐_𝑝𝑟𝑇_𝑟0−𝑇_𝑟−𝑈𝐴𝑇_𝑟−𝑇_𝑐+𝑉_𝑟−Δ_𝑟𝐻^𝑜𝑟𝑉_𝑐𝜌_𝑐𝑐_𝑝𝑐𝑑𝑇_𝑐𝑑𝑡=𝑞_𝑐𝜌_𝑐𝑐_𝑝𝑐𝑇_𝑐0−𝑇_𝑐+𝑈𝐴𝑇_𝑟−𝑇_𝑐, (22) where 𝑇 means the temperature, 𝜌 represents the density, 𝑐_𝑝 is the specific heat capacity, Δ_𝑟𝐻⁰ stands for the reaction enthalpy, 𝑈 is the overall heat transfer coefficient, and 𝐴 is the heat exchange area. Furthermore, meaning of subscripts is as follows: 0 is for the feed, 𝑐 for the cooling medium, and 𝑟 for the reaction mixture. Equation (22) represents a standard in CSTR design. Interested readers can find the specific values of all constant parameters and steady-state inputs of the CSTR in tables in [4, 7].

In fact, the physical parameters (such as reaction enthalpy, pre-exponential factor, and overall heat transfer coefficient) of this CSTR are not known exactly but they are supposed to vary within some intervals (see again [4, 7]). However, this paper particularly takes advantage of the final mathematical model of the CSTR introduced in [4], where it is obtained in the linearized form of an interval transfer function:𝐺()=𝑏𝑠,𝑏,𝑎²𝑠²+𝑏₁𝑠+𝑏₀𝑠⁴+𝑎₃𝑠³+𝑎₂𝑠²+𝑎₁𝑠+𝑎₀, (23) with parameters which can vary within the following bounds:𝑏₂∈⟨⟩𝑏−0.0291,−0.0245₁∈⟨⟩𝑏−0.0199,−0.0127₀∈⟨⟩𝑎−0.0005740,−0.0003549₃∈⟨⟩𝑎0.5801,0.9030₂∈⟨⟩𝑎0.1002,0.2299₁∈⟨⟩𝑎0.0062,0.0142₀∈⟨0.0001094,0.0002412⟩. (24)

5.2. Control Experiments

The interval transfer function (23) with parameters (24) describing the CSTR is considered to be the controlled plant. The first of its sixteen Kharitonov plants (13) is𝐺_1,1(𝑠)=−0.0245𝑠²−0.0199𝑠−0.000574𝑠⁴+0.903𝑠³+0.2299𝑠²+0.0062𝑠+0.0001094. (25) Corresponding closed-loop characteristic (5) can be computed as𝑠⁵+0.903𝑠⁴+0.2299−0.0245𝑘_𝑃𝑠³+0.0062−0.0199𝑘_𝑃−0.0245𝑘_𝐼𝑠²+0.0001094−0.000574𝑘_𝑃−0.0199𝑘_𝐼𝑠−0.000574𝑘_𝐼=0. (26) which means that the matrix (8) takes the following form:−𝑓𝑀=01000001000001000001₀𝑘_𝑃,𝑘_𝐼𝑓₅𝑘_𝑃,𝑘_𝐼−𝑓₁𝑘_𝑃,𝑘_𝐼𝑓₅𝑘_𝑃,𝑘_𝐼−𝑓₂𝑘_𝑃,𝑘_𝐼𝑓₅𝑘_𝑃,𝑘_𝐼−𝑓₃𝑘_𝑃,𝑘_𝐼𝑓₅𝑘_𝑃,𝑘_𝐼−𝑓₄𝑘_𝑃,𝑘_𝐼𝑓₅𝑘_𝑃,𝑘_𝐼, (27)where−𝑓₀𝑘_𝑃,𝑘_𝐼𝑓₅𝑘_𝑃,𝑘_𝐼=0.000574𝑘_𝐼−𝑓₁𝑘_𝑃,𝑘_𝐼𝑓₅𝑘_𝑃,𝑘_𝐼=−0.0001094+0.000574𝑘_𝑃+0.0199𝑘_𝐼−𝑓₂𝑘_𝑃,𝑘_𝐼𝑓₅𝑘_𝑃,𝑘_𝐼=−0.0062+0.0199𝑘_𝑃+0.0245𝑘_𝐼−𝑓₃𝑘_𝑃,𝑘_𝐼𝑓₅𝑘_𝑃,𝑘_𝐼=−0.2299+0.0245𝑘_𝑃−𝑓₄𝑘_𝑃,𝑘_𝐼𝑓₅𝑘_𝑃,𝑘_𝐼=−0.903. (28)

The stability boundary is determined by (11). The positions of such pairs (𝑘_𝑃,𝑘_𝐼) which fulfill (11) are plotted in Figure 2. The decision on area of stability and instability can be simply done with the assistance of an arbitrary testing point from the appropriate set. Moreover, the figure is supplemented with the half line 𝑘_𝐼=0 bordering the “upper” part of the stabilizing area.

Figure 2 

Areas of stability/instability for the Kharitonov plant (25).

Further, we must repeat the analogical procedure for all sixteen Kharitonov plants. The stability areas for the Kharitonov plants are depicted in Figure 3, and its zoomed version with highlighted intersection of all particular areas of stability is shown in Figure 4. The highlighted area represents the final region of robustly stabilizing PI controller parameters for the original interval transfer function of CSTR (23), (24).

Figure 3 

Areas of stability for sixteen Kharitonov plants.

Figure 4 

Zoomed robust stability region for the interval system (23), (24).

The following question is how to find the practically convenient PI controller from the obtained robust stability region. This paper utilizes the algebraic-based approach outlined in the Section 5. Nevertheless, this method requires a first-order fixed nominal model (17) of controlled plant in order to compute the final controller of appropriate (first) order, that is, with PI structure. So, the paper employs very simple but effective method of obtaining such model. In the first step, the fourth-order model with fixed parameters has been computed using the average values of interval parameters in (23), (24):𝐺_𝐴(𝑠)=−0.0268𝑠²−0.0163𝑠−0.000464450.74155𝑠³+0.16505𝑠²+0.0102𝑠+0.0001753. (29) Next, the simplest possible approximation takes advantage of neglecting the higher than zero-order powers of 𝑠 in numerator and higher than first-order powers of 𝑠 in denominator, that is, the corresponding coefficients are supposed to be zero. This approximation is easy but efficient enough as can be seen in the following steps. Obviously, it results in a first-order nominal model suitable for the applied synthesis:𝐺_𝑁(𝑠)=−0.00046445=0.0102𝑠+0.0001753−0.045534𝑠+0.017186. (30)

Further, (16) has been used for the calculation of PI controller parameters while the tuning parameter 𝑚 has been chosen according to??(18)𝑚=2.14⋅𝑎₀=2.14⋅0.017186≐0.036778, (31) in order to obtain 3% first overshoot for the nominal case. Consequently, the resulting controller is given by𝐶(𝑠)=−1.238𝑠−0.029706𝑠. (32)

This controller is located inside the stability region as can be seen in Figure 5. It means that the regulator robustly stabilizes the CSTR (23), (24). Besides, the controller lies on the curve hypothetically connecting the other potential controllers tuned by various parameters 𝑚>0.

Figure 5 

Position of controller (32) in stability region.

Finally, the robust stability and control performance is confirmed and demonstrated in Figure 6. It shows the output signals of the control loop with designed PI controller (32) and 128 “representative” systems from the interval family describing the CSTR (23), (24). The minimum and maximum values of each interval parameter have been used. It results in 2⁷=128 systems for simulation. On top of that, the red curve represents the output for the nominal system (30). Furthermore, Figure 7 depicts corresponding 128 + 1 control (actuating) signals. As far as control quality is considered, the results from Figure 6 might not be impressive at first sight, but there is a tradeoff between the simplicity of applied control algorithm and the performance here. As can be seen, the only one off-line tuned feedback PI controller with fixed coefficients has been utilized for controlling the chemical reactor with all the possible variations of uncertain parameters. We can observe that the CSTR is robustly stabilized successfully.

Figure 6 

Output signals of 128 “representative” plants and nominal system.

Figure 7 

Control signals for 128 “representative” plants and nominal system.

6. Conclusion

The main aim of the paper has been to present a possible approach to robust stabilization of a CSTR with exothermic reaction modelled as a fourth-order interval system. The developed easy but effective method of PI controller design has combined the Kronecker summation method with sixteen plant theorem and the algebraic tools. The proposed technique is applicable to a wide range of real processes, provided that they can be expressed by means of interval system and subsequently temporarily approximated by a first order model. The demerit of the method can be seen in a missing guarantee of coincident nominal performance and robust stability before the design process itself. They have to be verified during or after the design. Nevertheless, the applicability has been clearly demonstrated on the example where the CSTR has been successfully robustly stabilized.

Acknowledgments

The authors would like to acknowledge the support from the Ministry of Education, Youth, and Sports of the Czech Republic under Research Plan no. MSM 7088352102, the European Regional Development Fund under Project CEBIA-Tech no. CZ.1.05/2.1.00/03.0089, the Scientific Grant Agency of the Slovak Republic under Grant 1/0537/10, and the Slovak Research and Development Agency under Project VV-0029-07.