Control Structure Selection Based on Linear Quadratic Regulator

1. Introduction

Chemical processing plants are composed of complex and integrated structures that meet economic criteria, subject to capacity, equipment, and environmental constraints. For these reasons, it is necessary to define and implement controls that meet performance requirements, while maintaining stability and reducing the effects of unmeasured disturbances [1]. The plant-wide control area aims to broadly assess the problem of control in an industrial plant [2]. The plant-wide control of a process unit is intended to define many issues related, for example, with the pairings between the manipulated and controlled variables, with the most suitable controller type and with controller’s tunings. The definition of control structure, or the pairing between manipulated and controlled variables, is the most basic definition of plant-wide control, whenever it is not always easy to solve. This follows from the fact that chemical process plants generally have more than hundreds of measured variables and dozens of manipulated variables, which can cause a combinatorial explosion of possible control structures. In these scenarios, many methodologies for the definition of control structures have been proposed since the pioneering work of Buckley[3]. These methodologies can be divided into two classes:

• Based purely in heuristics;
• Based in mathematical models and control theory;

Methodologies based only in process heuristics [4,5], without mathematical models, have great acceptance in industrial area, because are easy to understand and systematize the control structure creation process. Although, these methodologies are severely influenced by subjectivity of the professional whose apply it in your engineering problem. In the other hand, methodologies which are based in mathematical models, [6,7,8,9] and control theory [10,11] have a great academic acceptance, but poorly applied in industrial area. Generally, the main difficult to apply these techniques in industrial area is the unknown or unavailable mathematical models during design period. In order to reduce the control structure search space, some methodologies applying some kind of heuristics, as the use of static models [12,13] or linearized models [11]. In the mathematical view the complete problem of plant-wide control, without simplifications and heuristics, is a nonlinear, constraint, uncertain and dynamic optimization problem. Some works have already treated a plant-wide control with this less simplified formulation [14,15]. However, even with development of dynamic optimization [16], the best control structure searching is almost an infeasible problem in some cases. The plant-wide control problem does not yet have a definitive formulation, increasing the distance between industrial design practice and academic developments. One way to control chemical process plants is trough more complex controllers where the relation between structure and measurements is less clear, as in the case of LQR (Linear Quadratic Controller) controllers. The LQR controller was developed firstly in the paper [17] and actually are widely studied and applied [18], [19], [20], [21]. These controller can be used in the whole plant [18] or in the distributed form [20], [21], where subsystems are controlled with LQR controllers. The LQR controller problem can be formulated as a unconstrained optimization problem [17] or through a constrained optimization problem [22], [23]. The LQR controller is used in a plenty of areas, for example, robotics [24], aeronautics [25], power electric systems [26], construction [27], biomedical [28] and process control [29]. In the case of process control, one of most popular cases, LQR controller is used as a base of predictive controllers [30] , as, for example, the algorithms of DMC (Dynamic Model Controller) [31].

2. Control Structure in LQR Problems

In the case of control structure definition where performance metrics and the type of controller is already defined, a new type of optimization problem is formulated, whose the solution is suboptimal, when compared with unconstrained problem (without controller type definition). When the controller type is defined and the mathematical structure used to combine measured is defined too, occurs a reduction in the dimension of control structure problem and consequently in the manipulated variables possible trajectories. In the point of view of control process, these equality constraints make it possible to solve the problem with feasible controllers and control structures. In the application sense is interesting simple control structures, based on linear combination, and simple controllers, based on feedback actions. In the next section is formulated the control structure definition problem based on Linear Quadratic Regulator (LQR) [17].

2.1. Optimal Control Structure in Linear Systems

Proposition 1.

In the case when a linear system is defined by differential equation as,

$${\displaystyle \frac{dx}{dt}}=Ax+Bu,$$

(1)

$$\mathrm{\Psi }=\mathrm{\Xi }x,$$

(2)

whose performance criteria is based on the metric (LQR),

$$J={\int }_0^{tf}{x}^{T}Qx+u^{T}Rudt,$$

(3)

suppose a pure proportional controller,

$$\mathrm{{\rm Y}}=I.$$

(4)

The optimal control structure, which minimize the metric J is,

$$\mathrm{\Gamma }=-R^{-1}{B}^{T}P{\mathrm{\Xi }}^{-1}.$$

(5)

Proof. 

In linear systems, for quadratic regulatory problems (LQR), the response of the optimal controller is known [17]. In this case the optimal response has an intermediary step, the solution of Riccati equation,

$${\displaystyle \frac{dP}{dt}}=A^{T}P+PA-PB{R}^{-1}{B}^{T}P+H^{T}QH,$$

(6)

where the optimal transformation, $T_x:X\to U$, is defined by,

$$u=-R^{-1}{B}^{T}Px,$$

(7)

defining that controller type is proportional, as showed at 4, the induced mapping between measured variables and manipulated variables, $T_{opt}:\mathrm{\Psi }\to U$, is numerically equal to the mapping between measured variables and controlled variables, ${\mathrm{\Gamma }}_{opt}:\mathrm{\Psi }\to \Omega$. Using the relation 2, the possible solution to the control structure problem is

$${\mathrm{\Gamma }}_{opt}=-R^{-1}{B}^{T}P{\mathrm{\Xi }}^{-1}.$$

(8)

□

Thus, for each linear system or linearized system, it is possible to conclude that there is an optimal linear combination of states (optimal control structure).

2.2. Optimal Control Structure in Nonlinear Systems for LQR Problems

The LQR controller is a linear solution, it is not possible that only one gain matrix meets the defined optimality criteria. Thus, a possible solution for this problem is linearization for each segment of problem [24,32], creating a gain matrix, structure, for each segment of the operational point. Follow below the specific proposition for the case of linearized system parts, exemplifying the case of control structures using the LQR metric.

Proposition 2.

Suppose a nonlinear system,

$${\displaystyle \frac{dx}{dt}}=\mathfrak{A}x+\mathfrak{B}u,$$

(9)

$$\mathrm{\Psi }=\mathfrak{T}\left(x\right),$$

(10)

whose performance criteria is based on the metric (LQR),

$$J={\int }_0^{tf}{x}^{T}Qx+u^{T}Rudt,$$

(11)

defining that the controller is proportional,

$$\mathrm{{\rm Y}}=I.$$

(12)

The set of optimal control structures ${\mathrm{\Gamma }}_{opt}$ for each linearizable segment is defined by,

$$\mathrm{\Gamma }(\alpha ,\beta ,\mu )=-R^{-1}\mathfrak{B}{\left(\beta \right)}^{T}P(\alpha ,\beta )\mathfrak{T}{\left(\mu \right)}^{-1}.$$

(13)

Proof. 

Let the set of matrices $\mathfrak{A}$, be defined by the following polytope [33],

$$\mathfrak{A}=\{\mathfrak{A}\left(\alpha \right)=\sum _{i=1}^M{\alpha }_i{A}_i,\sum _{i=1}^M{\alpha }_i=1,{\alpha }_i\ge 0\}.$$

(14)

Suppose that the set of matrices $\mathfrak{B}$, composed by polytope,

$$\mathfrak{B}=\{\mathfrak{B}\left(\beta \right)=\sum _{i=1}^N{\beta }_i{B}_i,\sum _{i=1}^N{\beta }_i=1,{\beta }_i\ge 0\}.$$

(15)

Suppose the set of matrices $\mathfrak{T}$ which represents the transformations between the state variables and measured variables $\mathfrak{T}:X\to \mathrm{\Psi }$, composed by the polytope,

$$\mathfrak{T}=\{\mathfrak{T}\left(\mu \right)=\sum _{i=1}^K{\mu }_i{\mathrm{\Xi }}_i,\sum _{i=1}^K{\mu }_i=1,{\mu }_i\ge 0\}.$$

(16)

For each pair of polytopes, $\alpha$ and $\beta$, the solution of optimal control structure is equivalent to the equation 8,

$$\mathrm{\Gamma }(\alpha ,\beta ,\mu )=-R^{-1}\mathfrak{B}{\left(\beta \right)}^{T}P(\alpha ,\beta )\mathfrak{T}{\left(\mu \right)}^{-1}.$$

(17)

□

The existence of one control structure, considered an optimum solution, based on linear combination of measured variables, cannot be used as a general solution for all cases. A unique solution, based on linear combination of variables, is mathematically acceptable if and only if the model can be represented with $\alpha$ and $\beta$ constants in all domain of problem. In the next section will be presented a methodology based on LQR controllers to define the control structure in the process plant.

3. Control Structure and LQR Controllers

As showed in Section 2 to a linear system whose performance metric is based on average of quadratic error weighted by control actions and the weighted difference between controlled states and yours respective set points (LQR problem) the optimal control structure, in the case of proportional controller, is defined by the matrix gain (Equation (5)). Thus, for linear problems, where the LQR performance criteria is the metric interest, it is naturally defined a linear combination of measured variables which it has characteristic of control structure. Generally process plants are represented by nonlinear models, being necessary the analysis when used only one control structure. As defined in the Section 2.2 for nonlinear systems, systems represented by more than one linear model, according to performance metric LQR, there are more than one optimal control structure, whose structure is constructed based on linear combination of variables. Two possible solutions which arising from these conclusions are the control structure definition for a plenty of plant operational points or a unique suboptimal structure in some of plant operational points. Since, to apply a unique control structure is necessary stability and performance loss analysis due to non-optimality according LQR criterion. The LQR controllers need linearized models in the state space form. In the case of more complex dynamic systems, which has implicitness nonlinearities inside algebraic differential equations, arising from algebraic and iterative thermodynamic equilibrium relations, for example, make complicated the linearization through algebraic handles. In these cases application of numerical algorithms are more suitable. In the next section is discussed a new methodology based on LQR controllers and linearized models for the solution of control structure definition problems.

4. Minimum Control Structure Based on LQR Gain Matrix

Prior defined that the control structure is based on linear combination of variables, which the controllers are proportional type and the performance criterion is LQR, as showed, a natural optimal control structure is defined based on LQR gain matrix. The closed loop system which use a LQR controller reach the maximum performance, in terms of objective function, and consequently in dynamic behavior requisite. The dynamic behavior of a physics system is dependent of its eigenvalues, in the case of closed loop control system, control structure and controller, interfere in the resultants eigenvalues. For these systems, in the linear case, the follow matrix shall be analyzed,

$$A_{mf}=A-K.B.$$

(18)

Where in the associated transformation matrix K are aggregated in the transformation of control structure and the specified controllers. As showed in the Equation 18, in a closed loop system the control structure and the controllers integrate eigenvalue definition. Thus, some modification in control structure or parameters of controllers can modify the region of eigenvalues. For linear systems, whose performance criterion is based on LQR objective function type, often is established a control structure where full matrix is necessary, which means that all combined measured variables shall be used as controlled variables. However, in industrial applications this kind of variables pairing, generally, it is not the more suitable than a one by one paring, even with performance criterion loss. Thus, a good control structure choice shall be a control structure with low loss of performance and meets the systems constraints, for example, stability. In this sense, the performance of a control structure, based on linear combination of variables, has as a maximum limit of dynamic behavior which is associated to respective LQR controller. Thus, a control structure plus controller can be considered near of optimum performance, if the respective eigenvalues is near of closed loop system with its respective LQR controller. In this way, a good approximation can be defined by the proximity between eigenvalues of closed loop and of the eigenvalues of its respective system controlled by LQR. With the proposal of increase the number of null elements in the gain matrix, corresponding to control structure $\mathrm{\Gamma }$, is necessary to create an element exclusion process to evaluate each gain matrix element. One way to accomplish this process is through of a evaluation of closed loop eigenvalues sensibility in each element of gain matrix. Thus, can be made an analysis of gain matrix with the eigenvalues in closed loop in the following form,

$$d_{k_{ij}\lambda }=\prod _{i=0}^n{\displaystyle \frac{d{\lambda }_i}{d{k}_{ij}}}.$$

(19)

This operation is made for each matrix gain element, composing a sensitivity matrix,

$$D_{k_{\lambda }}=\left(\begin{array}{ccc}{d}_{k_{00}\lambda }& ...& d_{k_{0n}\lambda }\\ .& ...& .\\ .& ...& .\\ d_{k_{m0}\lambda }& ...& d_{k_{mn}\lambda }\end{array}\right).$$

(20)

For each matrix line $D_{k_{\lambda }}$ exist a manipulated variable correspondence, where the element of each that are most sensitize the eigenvalues must be maintained, because this element of gain matrix is the element that most degrades the performance, thus the new matrix is the following,

$$K_m=\left(\begin{array}{ccccc}0& ...& k_{0j}\left(max\left(d_{k_{0j}\lambda }\right)\right)& ...& 0\\ .& ...& .\\ .& ...& .\\ 0& ...& k_{mj}\left(max\left(d_{k_{mj}\lambda }\right)\right)& ...& 0\end{array}\right).$$

(21)

The product between the eigenvalues of closed loop system matrix can be interpreted as volume of an ellipsoid with dimension n [34], where each eigenvalue represents axes of the ellipsoid. In this sense, the gain matrix gain which must be selected is matrix gain with minimum difference between its ellipsoid volume and original ellipsoid volume.

4.1. Stability in Closed Loop

This new gain matrix does not guarantee that the constraint of stability is attended, thus it is necessary to verify the new region of eigenvalues in closed loop. Case the stability criterion is not attended, the control structure shall be re-evaluated. In this case, a new gain matrix ($K_m$), which represents the structure, shall be proposed. It is known that if all elements of matrix $K_m$ are refilled with the elements of K, the system stability is guaranteed and that stability condition for discrete systems is,

$$\underset{K_m\to K}{lim}\left|\right|\mathrm{\Lambda }\left(A_{mf}\right)\left|\right|<1.$$

(22)

In the limit, it is possible to conclude that,

$$K_m=K,$$

(23)

where stability is guaranteed. Thus, exist some less complex control structure between $K_m$ and K. This intermediary control structure can be found through the repopulation process of gain matrix $K_m$ by the elements of original gain matrix of LQR controller K, which approximates continuously the eigenvalues of the matrix of closed loop system, which use the matrix $K_m$, to eigenvalues of original closed loop system. Following in Proposition 3 the proof of continuously approximation, in each iteration of repopulation process, between eigenvalues of intermediary system and original system.

Proposition 3.

The distances between the eigenvalues of the original system and intermediary systems are continuously decreasing when the repopulation process of the gain matrix based on the LQR system is applied.

Proof. 

Knowing that the distance between the characteristic values of two matrices is defined by [35],

$$d_H\le k\left(\begin{array}{c}{k}^{-1}{\left(\begin{array}{c}{\displaystyle \frac{\left|\right|A_{lqr}-A_{lq{r}_i}\left|\right|}{1-\rho \left(A_{lqr}\right)\left|\right|A_{lq{r}_i}\left|\right|}}\end{array}\right)}^{{\displaystyle \frac{1}{2m}}}\end{array}\right),$$

(24)

where k is elliptic modulus defined in [35] as

$$k\left(q\right)={\left(\begin{array}{c}{\displaystyle \frac{{\sum }_{r=-\infty }^{\infty }{q}^{{(r+0.5)}^2}}{{\sum }_{r=-\infty }^{\infty }{q}^{r^2}}}\end{array}\right)}^2,$$

(25)

that in the worst case is equal to one. Where the coefficient m is the minimum degree of polynomial associated to matrix $A_{lrq}$ and $\rho$ is the spectral radius of closed loop matrix, which is always inside the unitary circle in the LQR case,

$$\rho \left(A_{lqr}\right)\le 1.$$

(26)

Knowing that the norm between the matrix $A_{lq{r}_i}$ is always increase in each interaction in the repopulation process,

$$\left|\right|A_{lq{r}_i}\left|\right|\le \left|\right|A_{lq{r}_{i+1}}\left|\right|,$$

(27)

and that the norm of difference between the original matrix in each iteration is decreasing,

$$\left|\right|A_{lqr}-A_{lq{r}_{i+1}}\left|\right|\le \left|\right|A_{lqr}-A_{lq{r}_i}\left|\right|,$$

(28)

is true said that,

$$d_{H_{i+1}}\le d_{H_i}$$

(29)

Thus is proven that the distance $d_{H_i}$ in the repopulation process is convergent and decrease continuously to the original LQR control structure. □

Thus, using Proposition 3, it is possible to suppose that there exists an intermediary matrix between $K_m$ and K that has the eigenvalues within the stability limits. The proposed procedure to define control structures based on LQR can be summarizded in the Figure 1. In this procedure, there exists the possibility of using a control structure with a full matrix ($\mathrm{\Gamma }$) or a minimum control structure.

The proposed procedure has some similarity with LQG controllers type (Linear Quadratic Gaussian) [36] and $H_2$ [37]. The first two steps of the proposed procedure are similar to the definition of LQG controllers or the procedure using metrics $H_2$, the identification step and the definition step of the controller. The difference is in the mode of solution of the problem because in these two methodologies the state estimation problem and the control structure definition problem have the same metric, unlike the proposed methodology where the method of state estimation is not defined. Another different point, when compared with these two controllers, is the simultaneity of estimation that occurs in controllers case; however, it is made only once in control structure definition. In a general manner, this procedure of gain matrix definition, or in other words control structure definition, try to guarantee through one by one pairing a closed-loop system with a performance close as possible of optimal closed-loop system in the LQR sense. A negative point of this methodology, similar to other linear methodologies, is the sensitivity to nonlinear systems. If the values of the matrix elements $A_{mf}$ are distant, in each operational segment, the differences between the eigenvalues of the design point to the operational point (distinct from the design point), and consequently the possibility of an unstable system or an inappropriate performance of the system are greater. In the next section, we will follow some application examples of the proposed methodology.

5. Application of Proposed Methodology

5.1. Application to Stabilizable Systems in Closed Loop

The LQR controller guarantees that the system is stabilizable in closed loop if the system is observable and controllable [36]. Gain matrix without some elements, which is the case of control structure with pairing one by one, the stability is not necessarily guaranteed. In this section are presented examples in which the process of repopulation is necessary to guarantee stability.

5.2. Stabilizable System in Closed-Loop through Repopulation

As commented previously in this paper, initial selection of controlled variables using the proposed methodology can lead to closed-loop unstable systems. In this example, initially the system is unstable in closed loop and the process of repopulation of gain matrix is applied to stabilize the closed-loop system.

5.2.1. Example

Suppose an arbitrary system, represented by a state space model,

$$\dot{x}=Ax+Bu,$$

(30)

where the matrix A is defined by,

$$A=\left(\begin{array}{cccc}-0.3286& -0.3384& -0.3692& -0.2536\\ -0.3384& -0.5081& 0.4105& 0.1170\\ -0.3692& 0.4105& 0.2445& -0.1155\\ -0.2536& 0.1170& -0.1155& -0.5342\end{array}\right).$$

(31)

and the matrix B is defined by,

$$B=\left(\begin{array}{cccc}-0.4419& 0& 0& -0.3313\\ 1.7727& 0.6247& 0.4195& 0\\ 0& 1.2614& 0.1706& -0.8419\end{array}\right).$$

(32)

The LQR matrix computing is made with the usual procedure, where initially is solved an algebraic Riccati equation to find P matrix [38] from Equation (5), necessary matrix in the definition of LQR controller gain matrix. Solving the matrix gain of LQR controller,

$$K_{lqr}=\left(\begin{array}{cccc}1.0051& 0.1259& -0.1010& 0.1347\\ -0.1206& 0.5193& -0.1328& -0.0955\\ -0.1302& -0.3257& 0.5781& 0.1900\end{array}\right).$$

(33)

The sensitivity analysis is made according to proposed Equations (19) and (20). The found sensitivity matrix $D_{k_{\lambda }}$ is,

$$D_{k_{\lambda }}=\left(\begin{array}{cccc}0.0408& 0.0001& 0.0001& 0.0004\\ 0.0794& 0.1096& 0.0251& 0.0602\\ 0.5776& 0.8468& 0.0968& 1.0707\end{array}\right).$$

(34)

In this case, the control structure matrix, without repopulation process is,

$$K=\left(\begin{array}{cccc}1.0051& 0& 0& 0\\ 0& 0.5193& 0& 0\\ 0& 0& 0& 0.1900\end{array}\right).$$

(35)

The system eigenvalues with one by one pairing are different of original eigenvalues of LQR closed loop controlled system. In this first approximation, the stability constraint is not observed, which leads the closed-loop system to be unstable. In Figure 2 follow the eigenvalues of the closed-loop system with the original LQR controller and the eigenvalues of the closed-loop system with a reduced one-by-one controller.

As shown in the diagram poles, the temporal response of the closed-loop system with one-by-one reduced controller presents unstable behavior. Follow in Figure 3 the temporal response of a closed-loop system with one-by-one reduced dimension.

The stabilization of the system is applied to the repopulation process, as shown in Section 4.1. In this procedure the next elements, ranked by sensitivity to eigenvalues, are inserted in the matrix until stability is not reached. Follow in Equation (36) the new matrix after the repopulation process,

$$K=\left(\begin{array}{cccc}1.0051& 0& 0& 0\\ 0& 0.5193& -0.1328& 0\\ 0& -0.3257& 0& 0.1900\end{array}\right).$$

(36)

In this case, the repopulation process computes two iterations to reach closed loop stability. After the repopulation process, as expected, the closed-loop poles were displaced into a unitary circle and the approximate closed-loop system with the original LQR gain matrix, as shown in the Figure 4.

As expected for the poles’ locations, the closed-loop temporal response has stable behavior. For the temporal test, the closed loop system is tested from initial state,

$$X_0=\left(\begin{array}{ccccccccc}0.1& 0.1& 0.1& 0.1& 0.1& 0.1& 0.1& 0.1& 0.1\end{array}\right),$$

(37)

to the state,

$$X_1=\left(\begin{array}{ccccccccc}0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right).$$

(38)

Follow in the Figure 5 the temporal response of the system after repopulation process.

The temporal response for the manipulated variables of closed loop system is illustrated in the Figure 6.

As expected, the closed-loop system has become stable, with an intermediary control structure, between the system controlled by one-by-one pairing and the system controlled by LQR. In this numerical example, the control structure has reduced dimension when compared with the original control structure, but as expected, the performance is suboptimal. In this regard, the methodology emphasizes more stabilization than performance.

5.3. Application to Van de Vusse Reactor

The Van de Vusse reactor [39] is reference control process example, where the cyclopentanol (B) is produced from the cyclopentadiene (A) and two by-products are generated, cyclopentanediol (C) and dicyclopentadiene (D), in follow form,

$$A\stackrel{\mathit{K}\mathit{1}}{\to }B\stackrel{\mathit{K}\mathit{2}}{\to }C,$$

(39)

$$2A\stackrel{\mathit{K}\mathit{3}}{\to }D.$$

(40)

This chemical process has 3 manipulated variables (cooling jacket heat transfer, $Qk$ in kJ/h; reactor inlet flow, $F_{in}$ em l/h; and reactor outlet flow, $F_{out}$ em l/h) and 6 process measurements (reactor temperature, T, in K; cooling liquid temperature , $T_K$ in K; concentration of component A in outlet flow, $C_A$ in mol/l; concentration of component B in outlet flow, $C_B$ in mol/l; volume of product inside of the reactor, $V_R$ in l). Figure 7 represents the described process.

The dynamic model has tree differential equations to formulate the mass balance and two equations for the energy balance. In the following, follow the dynamic model of the van de Vusse reactor based on paper [11]:

• Mass Balance:

  $${\displaystyle \frac{dVr}{dt}}=F_{in}-F_{out},$$

  (41)

  $${\displaystyle \frac{d\left(Vr{C}_a\right)}{dt}}=F_{in}{C}_{ain}-F_{out}{C}_a-Vr(k_1{C}_a+k_3{C}_a^2),$$

  (42)

  $${\displaystyle \frac{d\left(Vr{C}_b\right)}{dt}}=-F_{out}{C}_b-Vr(k_1{C}_a+k_2{C}_b).$$

  (43)
• Energy Balance:

  $$\begin{array}{c}\hfill {\displaystyle \frac{dT}{dt}}=F_{in}{T}_{in}-F_{out}T+{\displaystyle \frac{k_w{A}_r}{\rho C_p}}(T_k-T)-\\ \hfill {\displaystyle \frac{Vr}{\rho C_p}}(k_1{C}_a\Delta H_1+k_2{C}_b\Delta H_2+k_3{C}_a^2\Delta H_3),\end{array}$$

  (44)

  $${\displaystyle \frac{d{T}_k}{dt}}={\displaystyle \frac{1}{m_k{C}_{pk}}}(Q_k+k_w{A}_r(T-T_k)).$$

  (45)
• Perturbations:

  $$C_{ain}=f_4\left(t\right),$$

  (46)

  $$T_{in}=f_5\left(t\right).$$

  (47)
• Control action:

  $$F_{in}=f_1(V_r,C_a,C_b,T,T_k),$$

  (48)

  $$F_{out}=f_2(V_r,C_a,C_b,T,T_k),$$

  (49)

$$Q_k=f_3(V_r,C_a,C_b,T,T_k).$$

(50)

Unlike the non-phenomenological problem shown previously, where there is a state space model, a step of obtaining the state-space model is needed. This step can be realized through the linearization of the phenomenological model around an operational point or through an identification problem facing the system as a black box. In this case, the identification problem was the solution chosen. As in a usual black-box identification process, before parameter estimation, a system excitation is needed. System excitation was performed with the PRBS (pseudorandom binary sequence) signal type [40]. Follow in Figure 8 the applied excitation signals.

After the excitation process, an estimation of the parameters was made, based on PEM (Predictor Error Estimate) method [41]. This method has some convergence problems as a function of the initial estimative choice, which is typical of deterministic optimization algorithms based on the gradient descendent. Follow in the Equations (51) and (52) the initial solution for the matrix A,

$$A_{ini}=\left(\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0.5& 0.5& 0& 0& 0\\ 0.5& 0.5& 0.5& 0& 0\\ -0.5& 0& 0& -0.5& 0.5\\ 0& 0& 0& 0.5& -0.5\end{array}\right),$$

(51)

and for the matrix B,

$$B_{ini}=\left(\begin{array}{ccc}1& -1& 0\\ 5.2& 1& 0\\ 0& -1& 0\\ 1& -1& 0\\ 0& 0& 1\end{array}\right).$$

(52)

The linearized model found after identification process was:

$$A=\left(\begin{array}{ccccc}-0.1066& 0.1600& 0.2943& -0.1452& 0.2079\\ 0.7536& 0.2879& -0.2930& -0.1679& -0.0774\\ 0.1535& 0.1987& 0.6447& -0.2129& -0.0921\\ -0.2577& 0.1769& 0.3152& -0.3012& 0.5500\\ -0.1855& 0.0611& 0.1912& 0.5304& -0.5210\end{array}\right).$$

(53)

$$B=\left(\begin{array}{ccc}0.8188& -0.4044& -0.0175\\ 4.9948& 0.8004& 0.0412\\ -0.0088& -0.6672& 0.0289\\ 0.8749& -0.5620& -0.0951\\ 0.0410& -0.0063& 0.1643\end{array}\right).$$

(54)

With the linearized system available and weighing matrices (Q and R) defined as identity, the expression (8) can be applied. The following LQR gain matrix is found,

$$K_{lqr}=\left(\begin{array}{ccccc}0.1990& 0.0282& 0.0189& -0.0086& -0.0489\\ 0.1847& 0.4929& -0.2797& 0.2629& -0.4655\\ -0.0101& -0.0941& 0.9456& 0.1953& -0.2463\end{array}\right).$$

(55)

Given the gain matrix, the algorithm of variable selection can be started. To obtain the matrix $D_{k_{\lambda }}$, the Equations (20) and (19) are applied,

$$D_{k_{\lambda }}=\left(\begin{array}{ccccc}327.853& 33184.553& 15.101& 537669.368& 7329.291\\ 230.103& 1.341& 1.091& 3.4686e-07& 0.304\\ 1.468e-10& 3.545e-06& 8.866e-06& 3.001e-10& 2.395\end{array}\right).$$

(56)

Defined the matrix $D_{k_{\lambda }}$, the selection procedure of the most sensitivity elements can be made, obtaining the following control structure,

$$K=\left(\begin{array}{ccccc}0& 0& 0& -0.0086& 0\\ 0.1847& 0& 0& 0& 0\\ 0& 0& 0& 0& -0.2463\end{array}\right),$$

(57)

which corresponds to the following parameters:

$$F_{in}\to T;$$

(58)

$$F_{out}\to V;$$

(59)

$$Q_k\to T_k.$$

(60)

In this case, using a state-space model, the system presents a stable behavior, with closed loop poles are inside unitary circle. The position of the closed-loop poles, eigenvalues, is presented in Figure 9.

As predicted by the location of the poles, the state has a stable behavior, where the states migrate asymptotically to other points. The liquid volume in the reactor remains within the operational limits, and the limits of the reactor trip level have not been reached.

The concentrations are stable in a new position and are also maintained within operational limits.

Figure 11. Temporal response of concentrations inside of the reactor.

Figure 11. Temporal response of concentrations inside of the reactor.

Similarly to the case of concentration, the temperatures stabilize in new values and also within the operational limits.

Figure 12. Temporal response of the temperature control loops.

Figure 12. Temporal response of the temperature control loops.

The temporal response of the variables manipulated in the closed-loop system is shown in Figure 13.

In a general manner, the proposed methodology presented an adequate behavior, even with the use of the pure proportional controllers. The asymptotic behavior demonstrates the possibility of using the PID controller for offset correction, which provides a suitable regulatory control for industrial applications.

5.4. Application in the Tennessee Eastman Problem

The aim of the studied plant is to produce G and H from four reactants (A, C, D, E). More than two elements are present in this plant, a by-product F and an inert element B, completing a total of eight components. This chemical plant has five main equipment: reactor, product condenser, liquid vapor separator, recycle compressor, and product stripper column. In Figure 14 the process diagram is presented.

The gaseous reactants feeding the reactor, where its reacts and the products are formed in accordance with the following reactions:

$$A\left(g\right)+C\left(g\right)+D\left(g\right)\to G\left(liq\right)$$

(61)

$$A\left(g\right)+C\left(g\right)+E\left(g\right)\to H\left(liq\right)$$

(62)

$$A\left(g\right)+G\left(g\right)\to F\left(liq\right)$$

(63)

$$3D\left(g\right)\to 2F\left(liq\right)$$

(64)

All of these reactions are irreversible and exothermic, it is necessary a heat exchanger to chill this stream. The products and feeds that have not reacted leave the reactor as vapor, which are partially condensed in the condenser product and hence are sent to the liquid-vapor separator. In the separator, noncondensable components are recirculated to the reactor through a recycling compressor or purged. The liquid components are fed into the stripper column, where the remaining reactants are separated. The vaporized components are incorporated to recycle stream, and the liquid phase is sent as a process plant product. This process has 41 continuous and discrete measurements and 12 manipulated variables, in this problem is considered only the 22 continuous measurements. The variables are listed in Table 1 and Table 2.

The dynamic model of the Tennessee Eastman process has 20 perturbations, which are used in the evaluation of control strategies. In this paper, only the perturbation related with the ratio A/C step in stream 4 is used to evaluate the control strategy, keeping constant the composition in B.

5.5. Application of Methodology

The application process was similar to that performed in the other plants. Initially, open-loop identification was performed, except in inventories controls (reactor level, separator level, stripper level, and reactor pressure), which were tested in closed loop to avoid trips during identification. Similarly to the other cases, the plant excitation was made through PRBS signals. Follow in the Figure 15 the signals used during the excitation of the plant .

After the excitation, the state-space parameter estimation of the process was performed using the PEM method. To avoid numerical problems during the identification, the process data are normalized between 0 and 1, for this reason the state-space model is referent to the scaled model. To find the values of the real plants, the multiplication of the temporal response to the respective scales is necessary. In this process, the system matrices presented in Section 7 were found. Applying the LQR controller design process, the selection procedure, and the repopulation process of the gain matrix, the control structure matrix presented in the annex was found. The pairings found are presented in Table 3.

Some of the pairings are obviously as in the case of Feed A and C with the manipulated variable XMV4 (valve directly correlated to the feed flow A and C), the recycle flow with the manipulated variable XMV5 (valve directly correlated to recycle flow). The feeds A, D and E were paired as flow ratios. The feed A was paired with the manipulated variable XMV1 (valve direclty correlated to feed flow D), this pairing has a relation with the reaction 61. The objective of this pairing is keep the reaction 61 stoichiometry. The feed D was paired with XMV2 (valve directly correlated to the feed E), this pairing has the sense of keep the relation between H and G products. The stripper steam flow was paired with the manipulated variable XMV8 (valve directly correlated to stripper bottom flow), this pairing try to keep the ratio between stripper bottom flow and stripper steam flow, being this configuration usual in absorbers columns in the industry. Some of pairings are not of obvious interpretation, as in the case of separator level with the manipulated variable XMV12 (reactor agitator speed) and the compressor power rating was paired with the manipulated variable XMV9 (valve directly correlated with stripper steam), which has not process reason specifically, but was considered important to keep the eigenvalues in the vicinities of closed loop system controlled by LQR controller. In this procedure the optimization problem is faced as an unconstrained problem, which is one of the reasons why some pairings relates with the constraints are not appearing. Thus, it can be possible to conclude that the minimum pairing, with a proportional controller, which has a dynamic behavior closed as possible to the respective LQR controller, that not include necessarily the constraint variable. This method deficiency can be compensated with override strategies to the pairing with the constraints process variables.

6. Conclusions

The control structure based on LQR controller is a simple solution with easy application, since a state-space model is available. A positive aspect of this methodology is the prioritization of simple control structures (pairing one by one), using only linear combination of variables in the cases wherein the stability criterion is not met. In the cases wherein the systems are observable and controllable, the closed loop stability using the original gain matrix is guaranteed, being possible to find an intermediary matrix which became the closed loop system stable. The methodology was applied in three examples and feasible solutions were found, since in complex systems. The critical aspect of this methodology is the necessity of available state-space models, which is not available in every systems and in every operational point for nonlinear systems. The treatment of the constraints pairing is another point that can be investigated in future works. The robustness to the nonlinearity must be improved in the methodology, because in some cases performance loss can be greater than the tolerable limit and in worst case the instability can be reached. In these cases more than one control structure, based on LQR controller, can be used.

7. Matrices of Control Structures

Matrices of the Tennessee Eastman process

$$A=\left(\begin{array}{cccccccccccccccccccccc}0.227& 0.010& -0.013& -0.049& 0.060& 0.024& -0.039& -0.007& -0.012& -0.010& -0.019& 0.005& -0.075& 0.038& -0.005& -0.025& 0.015& -0.006& -0.004& 0.044& -0.024& -0.012\\ 0.033& 0.192& 0.029& -0.003& 0.005& 0.064& 0.052& 0.002& 0.014& 0.032& -0.010& 0.014& -0.019& -0.028& 0.034& 0.024& 0.024& 0.012& 0.027& -0.016& 0.006& -0.034\\ -0.053& -0.009& 0.171& 0.019& 0.035& 0.045& -0.063& -0.035& 0.106& 0.052& 0.077& 0.037& -0.048& -0.040& 0.100& -0.026& 0.027& 0.069& 0.013& 0.047& 0.016& 0.015\\ 0.007& 0.026& 0.027& 0.075& 0.044& 0.039& 0.083& 0.007& 0.002& -0.003& 0.063& 0.010& 0.087& 0.002& 0.036& 0.081& -0.008& 0.012& -0.032& 0.046& 0.069& 0.090\\ -0.013& 0.054& -0.010& 0.065& 0.025& 0.016& 0.063& 0.080& 0.010& 0.052& 0.067& 0.047& 0.083& -0.026& 0.068& 0.016& 0.009& 0.023& 0.012& 0.035& 0.000& -0.011\\ -0.049& 0.011& 0.026& 0.062& 0.075& 0.036& 0.114& 0.109& 0.018& -0.026& 0.034& 0.032& 0.105& -0.034& 0.001& 0.059& 0.004& -0.022& -0.026& 0.046& 0.018& 0.016\\ -0.009& 0.029& -0.073& 0.120& 0.075& 0.089& 0.099& 0.138& -0.003& 0.031& 0.004& 0.052& 0.138& 0.019& -0.068& 0.102& -0.089& 0.002& 0.070& 0.003& -0.005& -0.026\\ 0.003& -0.017& -0.059& 0.113& 0.047& 0.097& 0.093& 0.109& -0.068& 0.009& 0.062& -0.019& 0.079& -0.049& -0.001& 0.156& -0.062& 0.014& 0.018& 0.036& 0.050& -0.037\\ 0.020& 0.062& 0.053& 0.026& -0.021& 0.084& 0.037& -0.025& 0.019& 0.048& 0.014& 0.050& 0.002& 0.064& -0.010& -0.019& 0.087& 0.003& -0.026& -0.029& 0.089& 0.019\\ -0.009& 0.027& 0.018& 0.024& 0.060& 0.002& 0.075& 0.132& 0.035& 0.098& 0.011& -0.011& 0.077& 0.030& -0.032& 0.070& 0.010& -0.003& -0.012& 0.018& 0.044& 0.005\\ -0.026& 0.019& 0.063& 0.041& -0.007& 0.073& 0.028& 0.061& 0.066& 0.014& 0.013& -0.021& 0.028& 0.026& 0.038& -0.004& 0.068& -0.011& 0.023& 0.052& 0.004& 0.068\\ 0.072& 0.035& 0.068& -0.056& 0.016& 0.011& -0.071& -0.019& 0.090& 0.055& 0.025& 0.063& -0.097& 0.063& 0.104& -0.088& 0.052& 0.080& -0.037& 0.083& 0.025& 0.031\\ -0.034& 0.007& -0.022& 0.060& 0.111& 0.111& 0.053& 0.074& -0.021& 0.013& 0.071& 0.079& 0.055& -0.015& -0.023& 0.059& -0.050& -0.002& 0.050& 0.001& 0.056& -0.012\\ 0.038& 0.011& 0.007& -0.050& -0.078& -0.024& 0.034& -0.017& 0.082& 0.039& 0.055& -0.002& 0.045& 0.164& -0.016& 0.063& 0.094& -0.013& 0.043& -0.008& 0.013& 0.058\\ 0.062& 0.067& 0.083& -0.033& -0.023& -0.021& -0.046& -0.026& 0.069& 0.066& 0.021& 0.044& -0.070& 0.044& 0.063& -0.087& 0.118& 0.076& -0.044& 0.067& 0.051& 0.026\\ -0.047& 0.001& -0.072& 0.083& 0.061& 0.170& 0.174& 0.121& 0.006& 0.024& 0.046& 0.027& 0.165& -0.006& -0.019& 0.134& -0.095& 0.065& 0.052& 0.012& 0.023& 0.038\\ 0.034& -0.002& 0.054& 0.006& -0.066& 0.003& 0.024& -0.015& 0.023& 0.068& 0.031& -0.037& -0.007& 0.068& 0.011& 0.029& 0.222& 0.012& 0.008& -0.013& 0.049& 0.026\\ 0.041& 0.039& 0.129& 0.004& 0.037& -0.065& -0.065& -0.057& 0.065& 0.028& 0.057& 0.068& -0.031& 0.024& 0.125& -0.099& 0.059& 0.084& -0.021& 0.015& 0.013& 0.027\\ 0.014& 0.047& -0.039& -0.037& 0.035& -0.024& 0.077& 0.031& 0.000& 0.005& -0.033& 0.043& 0.021& 0.042& -0.046& 0.006& -0.013& 0.057& 0.272& 0.058& 0.051& 0.017\\ 0.082& 0.105& -0.006& 0.048& -0.008& -0.033& 0.048& 0.076& -0.008& 0.056& 0.065& 0.093& -0.041& -0.017& 0.053& -0.021& 0.045& 0.116& -0.013& 0.040& -0.002& 0.055\\ -0.010& 0.007& 0.024& 0.037& 0.025& -0.003& 0.028& -0.003& 0.006& 0.054& 0.071& 0.013& 0.032& 0.050& 0.015& -0.008& 0.043& 0.048& -0.039& -0.007& 0.006& -0.002\\ -0.021& -0.017& 0.039& 0.086& -0.032& 0.039& 0.031& 0.014& 0.037& 0.031& 0.027& 0.001& 0.051& 0.040& 0.031& -0.028& 0.064& 0.047& -0.011& 0.037& 0.012& -0.005\end{array}\right)$$

(65)

$$B=\left(\begin{array}{cccccccccccc}0.0233& -0.0611& 0.6625& 0.0361& -0.0138& 0.0032& 0.0045& -0.0059& 0.0218& 0.0352& 0.0139& 0.0029\\ 0.6321& -0.0024& 0.0289& 0.0289& 0.0222& 0.0156& -0.0315& -0.0168& 0.0063& -0.0413& -0.0175& 0.0135\\ 0.0303& 0.5262& -0.0374& -0.0382& -0.0052& -0.0378& 0.0408& -0.0146& 0.0436& 0.0270& -0.0013& -0.0589\\ -0.0342& -0.0385& 0.0090& 0.3627& -0.0418& -0.0381& -0.0216& -0.0135& -0.0370& -0.0160& -0.0472& -0.0360\\ -0.0069& -0.0522& -0.0260& 0.1311& 0.0953& 0.0567& -0.1014& -0.1086& 0.0229& 0.1227& 0.0315& 0.0603\\ 0.0253& -0.0120& -0.0772& 0.2055& 0.0082& -0.0436& -0.0141& -0.0184& 0.0207& 0.0507& -0.0121& 0.0432\\ 0.0043& -0.1061& -0.1378& 0.2629& 0.0294& -0.0291& 0.0179& -0.0014& 0.0110& 0.1512& 0.0125& 0.0893\\ 0.0188& -0.1677& -0.0959& 0.2035& 0.0968& 0.0775& -0.0499& -0.1078& 0.0182& 0.2104& 0.0715& 0.1261\\ -0.0061& 0.0838& -0.0261& 0.0860& -0.0198& -0.0176& 0.0589& 0.0572& -0.0126& -0.0507& -0.0145& -0.0241\\ -0.0297& -0.0535& -0.0297& 0.1150& -0.0249& 0.4012& -0.0374& 0.0043& 0.0060& -0.0393& -0.0643& 0.0256\\ 0.0152& 0.0582& -0.0305& 0.1478& 0.0075& 0.0229& 0.0301& 0.0580& -0.0067& 0.0018& -0.0062& 0.0263\\ 0.0417& 0.1823& 0.0564& 0.0040& 0.0326& 0.1550& -0.0041& -0.0360& 0.0444& 0.0586& 0.0568& 0.0418\\ -0.0361& -0.1026& -0.1094& 0.2234& 0.0546& -0.0268& -0.0275& -0.0481& -0.0036& 0.0865& -0.0085& 0.0443\\ -0.0740& 0.0124& 0.0192& -0.0265& -0.0604& -0.0106& 0.5654& 0.0870& -0.0083& -0.0084& 0.0075& 0.0195\\ 0.0242& 0.1702& 0.0776& -0.0441& -0.0168& 0.1317& 0.0796& 0.0773& -0.0595& 0.1004& 0.0338& -0.0057\\ 0.0218& -0.1361& -0.1309& 0.2617& 0.0446& -0.0609& -0.0208& -0.0563& -0.0078& 0.1440& 0.0039& 0.0702\\ 0.0314& -0.0118& 0.0117& -0.0270& 0.0445& 0.0013& 0.1354& 0.6072& -0.0338& -0.0256& 0.0173& -0.0358\\ -0.0124& 0.2608& 0.0335& -0.0193& 0.0118& 0.1406& 0.0319& 0.0149& 0.0525& 0.0665& 0.0641& 0.0351\\ -0.0269& -0.0704& -0.0675& -0.0626& -0.0267& -0.0191& -0.0051& -0.0108& 0.6968& 0.0330& 0.0066& -0.0638\\ 0.0558& 0.0869& 0.0669& 0.0133& 0.0430& 0.1914& 0.0146& -0.0193& 0.0379& 0.1792& 0.0996& 0.1181\\ -0.0072& 0.0720& -0.0175& 0.0861& -0.0247& -0.0147& 0.0365& 0.0827& 0.0202& -0.0632& -0.0245& -0.0221\\ -0.0001& 0.0233& 0.0018& 0.0746& -0.0241& -0.0293& 0.0406& 0.0533& -0.0016& -0.0708& -0.0231& -0.0226\end{array}\right)$$

(66)

Gain matrix of the LQR controller for Tennessee Eastman process

$$\mathrm{\Gamma }=\left(\begin{array}{cccccccccccccccccc}\mathbf{0}.\mathbf{932}& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000\\ 0.000& \mathbf{0}.\mathbf{856}& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000\\ 0.000& 0.000& \mathbf{0}.\mathbf{905}& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000\\ 0.000& 0.000& 0.000& \mathbf{0}.\mathbf{913}& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000\\ 0.000& 0.000& 0.000& 0.000& \mathbf{0}.\mathbf{994}& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000\\ 0.000& 0.000& -\mathbf{0}.\mathbf{020}& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000\\ 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& \mathbf{0}.\mathbf{012}& 0.000& 0.000\\ 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000\\ 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000\\ 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& \mathbf{0}.\mathbf{903}& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000\\ 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& -\mathbf{0}.\mathbf{003}& 0.000& 0.000& 0.000& 0.000\\ 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& \mathbf{0}.\mathbf{980}& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000\end{array}\right)$$

(67)