PID Control of DAE Systems - American Chemical Society

Ind. Eng. Chem. Res. 2003, 42, 4599-4610

4599

PID Control of DAE Systems Nammi Madhusudana Rao,† Prashant Vora,‡ and Kannan M. Moudgalya*,†,‡ Department of Chemical Engineering and Systems and Control Group, Indian Institute of Technology, Bombay, Powai, Mumbai 400 076, India

It is shown through two nontrivial examples that it is possible to control DAE systems with PID controllers. A two-phase reactor and a CSTR with a heating jacket are linearized about their operating points, and standard linear DAE control techniques are applied to arrive at transfer function matrices. PI controllers derived using Ziegler-Nichols and BLT methods are found to control these systems at least as well as geometric-theory based control techniques. Considering its simplicity, the proposed method should be added to the toolkit of control engineers. 1. Introduction The most popular tool used in industry even today is the PID controller.1 PID controllers are easy to implement. The infrastructure required for their implementation is minimal. Their workings can easily be explained to plant operators. These controllers also allow operators to tune one loop at a time. It is fair to say that the PID controllers have held their ground despite the several new control design techniques that have become available. Differential algebraic equation (DAE) systems naturally arise in dynamic systems. They are more difficult to handle than ordinary differential equations (ODEs) because of issues such as index and consistent initial conditions. Although one way to handle DAEs is to first convert them into ODEs and then apply the standard ODE techniques, there are compelling reasons for handling DAEs directly.2 The physical meaning of the variables could be lost during the conversion process. Also, one might obtain different models with solution manifolds of different dimensions. Finally, solving the original DAEs directly is usually much more efficient. Some of the popular methods of controlling DAE systems involve geometric control techniques. Extensive studies using this approach have been carried out by Kumar and Daoutidis3 and Krishnan and McClamroch.4-6 Although the latter did not consider the case of manipulated inputs appearing in algebraic constraints, the former did allow for this possibility: when the inputs are present in algebraic constraints, dynamic state feedback compensation is used. Except for this difference, the two approaches are similar. We address the problem of controlling DAE systems in a slightly different way. The nonlinear DAEs describing the dynamic system of interest are first linearized about their operating point. The linear DAE control technique of Byers et al.7 is applied to arrive at an index * To whom correspondence should be addressed. E-mail: [email protected]. Tel.: +91 22 2576 7039. Fax: +91 22 2572 3702. Internet: www.che.iitb.ac.in/faculty/km.html. † Department of Chemical Engineering. ‡ Systems and Control Group.

1 system. PID controllers, tuned with the industryfriendly method of BLT developed by Luyben1 are then used. This approach is illustrated with a two-phase reactor and a CSTR with a heating jacket, as studied by Kumar and Daoutidis.3,8 The efficacy of the proposed controller is demonstrated through a detailed comparative study. This article is organized as follows: In section 2, we summarize linear DAE control theory. In section 3, the algorithm used in this work is presented. Sections 4 and 5, respectively, address the control design problems of the two-phase reactor and the CSTR with a heating jacket. The last section provides conclusions. 2. Linear DAE Control Theory Consider the following linear time-invariant system of differential algebraic equations (DAEs)

E ˜ x˘ (t) ) A ˜ x(t) + B ˜ u(t) y(t) ) C ˜ x(t), t g 0

(1)

˜ ∈ Rn×n, and B ˜ ∈ Rn×m; with system matrices E ˜ ∈ Rn×n, A state x ) x(t) ∈ Rn; input u ) u(t) ∈ Rm; and output y ) y(t) ∈ Rp. The solution of DAE system 1 in the sense of distribution is given by9

x(t) ) P

P

[]

I A1t e [I 0 ]P-1x(0) + 0 t A1(t-τ) I e B1u(τ) dτ P 0 0

[ ]∑ 0 I

µ-1

[ ]∫

[ ]∑

δ(i-1)(t)Ni[0 I ]P-1x(0) - P

i)1

0 I

µ-1

NiB2u(i)(t)

i)0

(2) where

[]

[ ]

[ ]

x B A 0 P-1x ) x1 , PA ˜Q ) 1 , PB ˜ ) B1 0 I 2 2

x1 ∈ Rn1, x2 ∈ Rn2, n ) n1 + n2, and δ(t) is the Dirac delta function. Note that the solution includes impulse terms due to either initial conditions or possible jump behavior in the control input and input derivatives.

10.1021/ie0208671 CCC: $25.00 © 2003 American Chemical Society Published on Web 06/20/2003

4600 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003

Therefore, it is necessary to analyze the effect of control on impulse terms in the state response, which can be done through impulse controllability. Definition 1 (Impulse Controllability9). The system of equations 1 is said to be impulse-controllable if, for any initial conditions x(0) ∈ Rn, τ ∈ R, and w ∈ Rn2, such there exists an admissible control input u ∈ Cµ-1 p that

[ ]

0 ) Iτ(w,t) ) x2τ(t)

[

0 µ-1

δ(i-1)(t - τ)Niw ∑ i)1

]

3. Algorithm for Control Law Design

(3)

[

]

E ˜ 0 0 ) n + rank E ˜ A ˜ E ˜ B ˜

Lemma 1.9 If deg(|sE ˜ -A ˜ |) ) rank E ˜ , then system of equations 1 is impulse-controllable and has index 1. 0 If the system of equations 1 is impulse-controllable, then it can be reduced to an index 1 system by using the infinite pole placement technique.9 In many applications, however, the assumption of impulse controllability does not hold. The condensed staircase form is proposed to identify the modes that are not controllable at infinity and that cannot be reduced to index 1 by feedback.7 Theorem 1.7 For the DAE given in eqs 1, the following results hold: There exists a unitary state transformation by which eqs 1 can be reduced to the form

Eζ˙ ) Aζ + Bu y ) Cζ where

[

]

E11 0 E13 0 E23 , E ) UE ˜V ) 0 0 0 E33

(4)

[

A11 A12 A13 A ) UA ˜ V ) A21 A22 A23 0 0 A33

[] []

We consider a nonlinear DAE system of the form

x˘ ) f(x) + b(x)z + g(x)u

where N is nilpotent with the nilpotent index denoted by µ, the differential index of the DAE system. Impulse controllability characterizes the ability to generate impulse terms by admissible control. A necessary and sufficient condition for impulse controllability is9

rank

Central to this computation procedure is singular value decomposition. Note that the only way to handle uncontrollable modes at infinity is to ensure consistent initial conditions. We next present the algorithm used in this work to control nonlinear DAE systems.

]

ζ1 B1 B ) UB ˜ ) B2 , ζ ) ζ2 , C ) C ˜ V ) [C1 C2 C3] ζ3 0 with rank(E11) ) t1; rank(B2) ) t2; A33 being block upper triangular; and E33 being block upper triangular, having zero diagonal blocks, and being partitioned conformally with A33. There exists a state feedback gain matrix F ∈ Rm×n such that λE - (A + BF) is invertible, a requirement for unique solution, if and only if A33 is nonsingular. Then, for t > 0, ζ3(t) ) 0 is independent of the control input u. If E33ζ3(0-) ) 0, then ζ3 ) 0 and is impulse-free. 0

0 ) k(x) + l (x)z + c(x)u yi ) hi(x), i ) 1, 2, ..., m where x ∈ Rn; u ∈ Rm; z ∈ Rp; and k(x), f(x), h(x) and b(x), g(x), l(x) are vector fields and matrices of appropriate dimension, respectively. We also assume that all states are observable. The algorithm is presented next. 1. Linearize the nonlinear system around some operating point. We assume that the linear DAE system is either regular [i.e., det(sE ˜ -A ˜ ) * 0] or regularizable by state feedback [i.e., det(sE ˜ - (A ˜ +B ˜ K)) * 0]. 2. Apply the linear DAE control law as follows: (a) Check the impulse controllability of the resulting linearized system. If it is not impulse-controllable, then some modes of this system are not controllable at infinity. These modes do not change the system dynamics if consistency conditions are satisfied. (b) Separate modes that are not controllable at infinity by using theorem 2. These modes can be removed from the system. The reduced-order system will be controllable at infinity. (c) Move the poles that are at infinity to finite positions by using the infinite pole placement technique.9 The resulting system will be of index 1. 3. The index 1 DAE system can be converted to a standard state-space system, which can be controlled using any conventional technique. The approach used in this work to arrive at the PI controller can be summarized as follows:1 (a) Evaluate the transfer function matrix. (b) Pair inputs and outputs. (c) Find the Ziegler-Nichols settings for each loop. (d) Detune the loops through the BLT method. 4. Control of a Two-Phase Reactor Consider a reactor with liquid and vapor phases, as discussed in ref 8. Reactants A and B are fed to the CSTR as pure vapor and liquid streams, respectively, at molar flow rates FAo and FBo, while the two outlet streams from the liquid and vapor phases have molar flow rates FL and FV, respectively. It is assumed that the individual phases are well-mixed, that they are in physical equilibrium at pressure p and temperature T, and that the chemical reaction is slow compared to the mass transfer across the interface. A detailed description of this process is given in ref 8.

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4601

4.1. Process Model. The process description takes the form of the following DAE model

dNV ) FAo - NA + NC - FV dt

wi ) xi - x/i , wj ) zm - z/m, vk ) uk - u/k

( ) ( )

dyA FAo(1 - yA) 1 - yA yA ) NA N dt NV NV NV C

where i ) 1, ..., 6; j ) 7, ..., 9; m ) 1, ..., 3; k ) 1, 2; x/i , z/j , u/k are equilibrium points, wi, wj are new state variables; and v1, v2 are new input variables. The linearized model is as follows

dNL ) FBo - FL - RC + NA - NC dt

w˘ 1 ) -w7 + w8 - v1

xA dxA FBoxA + RC(1 - xA) 1 - xA + NA + N )dt NL NL NL C

] ( )

[

w˘ 2 ) 8.1565 × 10-7w1 - 0.0039w2 - 0.0221w7 0.0558w8

dxB FBo(1 - xB) - RC(1 - xB) xB N + ) dt NL NL A

(

w˘ 3 ) -0.0036w3 - 0.1944w4 - 0.0684w5 - 0.052w6

( )

xB N NL C

dT FAo(TAo - T) + FBo(TBo - T) + ) dt (NL + NV) RC

5 is linearized through a Taylor series expansion around the operating point using Matlab. Let

(5)

)

∆Hr ∆Hv T+ (N - NC) + cp (NL + NV) (NL + NV)cp A 1 Q (NL + NV)cp

w˘ 4 ) -2.1525 × 10-4w3 - 0.0408w4 - 0.0041w5 3.1236 × 10-4w6 + 0.0595w7 + 0.0186w8 w˘ 5 ) -9.0511 × 10-5w3 - 0.0049w4 - 0.031w5 1.3241 × 10-4w6 -0.0529w7 + 0.0524w8 w˘ 6 ) 2.695 × 10-5w1 - 0.0399w3 - 2.1494w4 -

(6)

0.7556w5 - 0.0746w6 + 9.7481w7 - 9.748w8 + 4.8741 × 10-4v2

0 ) -xAPsA + pyA

0 ) 16 985w2 - 51 099w4 - 504.4365w6 + 0.716w9

0 ) -(1 - xA - xB)PsC + p(1 - yA)

0 ) -16 985w2 + 56 751w4 + 56 751w5 -

(VTF - NL) 0 ) -NVRT + p F

142.4291w6 + 0.284w9 0 ) -2839.3w1 - 1132.4w3 - 106.7434w6 +

y1 ) yA

2.1462w9

y2 ) T In the above equations, NA is the molar rate of transfer of reactant A from the vapor to the liquid phase, NC is the molar rate of transfer of product C from the liquid to the vapor phase, Q is the heat input to the reactor, and VT is the reactor volume. The rate of formation of the product C is given by (-Ea/RT)

RC ) koe

CACBVL ) koe

y1 ) w2 y2 ) w6 This can be written as

E ˜ w˘ ) A ˜w + B ˜v y)C ˜w

(-Ea/RT)

NLFxAxB

(7)

Refer to Kumar and Daoutidis8 for the nominal values of the parameters. It is desired to control the composition of the vapor phase yA and the temperature T using the vapor-stream outlet flow rate FV and the heat input Q as the manipulated inputs. Define the differential variables (x), the algebraic variables (z), and the controlled (y) and manipulated (u) inputs as given below

4.3. Controller Design Using the Proposed Approach From

x1 ) NV, x2 ) yA, x3 ) NL, x4 ) xA

we see that the system of equations 7 is not impulsecontrollable. Using theorem 1, the staircase form of eqs 4 is attained. The resulting matrices are given in ref 10. It can be verified that, for the two-phase reactor system, det(sE - A) * 0. As a result, A33 is nonsingular from theorem 1. It follows that, in the two-phase reactor, if the consistency conditions are satisfied, the uncontrollable infinite poles of the system play no role in the dynamics and, hence, can be removed.

x5 ) xB, x6 ) T, z1 ) NA, z2 ) NC, z3 ) p y1 ) yA, y2 ) T, u1 ) FV, u2 ) Q 4.2. Linearized Model of the Two-Phase Reactor. We assume that the variables deviate only slightly from the nominal operating point. The nonlinear model in eqs

rank

[

]

E ˜ 0 0 ) 13 < n + rank E ˜ ) 15 A ˜ E ˜ B ˜


Figure 1. Closed-loop profiles of the controlled outputs and manipulated inputs for a 15% decrease in the set point for the mole fraction of species A in vapor phase: proposed approach (solid line) and that of ref 8 (dashed line).

The part of the system that is controllable at infinity is given by

[

E11 0 0 0

][ ] [

][ ] [ ] []

ζ˙ 1 A11 A12 ) ζ˙ 2 A21 A22

ζ1 B1 + ζ2 B2 v

ζ y ) [C1 C2] ζ1 2

u1 ) 18.4434(y10 - y1) + (8)

|[

][

∫

18.4434 (y10 - y1) dt + 231 0.05

u2 ) 132.7298(y20 - y2) +

∫

132.7298 (y20 - y2) dt + 75.96 100 (9)

where y10 and y20 are the outputs of low-pass filters

It can be verified that

A A sE11 0 deg - A11 A12 0 0 21 22

to this procedure, a PI controller of the following form was designed

]|

) 4 ) rank E

and hence, from lemma 2, it can be concluded that the system of equations 8 is impulse-controllable and that it has index 1. From the index 1 DAE system of equations 8, we arrive at a 2 × 2 transfer function matrix between v and y, listed in Appendix A. Any standard control design technique can now be applied to this system. In this work, we use the Ziegler-Nichols method to design controllers for individual loops and the BLT method for detuning.1 The detuning factor F is chosen to be 2. In addition, the reference signal is passed through a lowpass filter, as the controller gains are high. According

y10 )

ysp1

ysp2 , y20 ) 1500s + 1 1500s + 1

Here, ysp1 is the set point for y1, and ysp2 is the set point for y2. The efficacy of this controller in handling the nonlinear system of equations 5 is verified next. 4.4. Comparison of Controller Performance. The performance of the PI controller is now compared with that of the controller used in ref 8, which was derived through a two-step procedure: (1) The DAE system was first converted into an ODE system through differentiation and simplification. (2) A global linearization technique was then employed to determine the controller. We will refer to this controller as the geometric controller. The state-of-the-art symbolic mathematical package Mathematica11 was used for all of the simulations


Figure 2. Closed-loop profiles of the controlled outputs and manipulated inputs for a 2.5% increase in the set point for the temperature: proposed approach (solid line) and that of ref 8 (dashed line).

reported in this work. The simulations were carried out on a Compaq Tru64 UNIX system, consisting of 16 833MHz Alpha Processors with 32 GB of RAM. As suggested by Kumar and Daoutidis,3 the performance of the controllers for set-point changes and disturbance rejection were evaluated. The first three simulation runs address the set-point tracking capabilities of the PI controller and the geometric controller. In the first run, at time t ) 5 min, the process at its nominal steady state is subjected to a 15% decrease in ysp1 or, equivalently, an increase in the desired mole fraction of the product C in the vapor phase. The corresponding input and output profiles under the two controllers are shown in Figure 1. Clearly, the PI and geometric controllers provide good responses to setpoint changes in y1 and also decoupled responses, as required. In the second run, the process is subjected to a 2.5% increase in ysp2, that is, an increase in the desired reactor temperature T. Again, both controllers provide a good response for y2 and also a decoupled response, as required. The output and manipulated input profiles are shown in Figure 2. In these two runs, the performances of the two controllers are more or less similar. Moreover, these results are also in agreement with those of ref 8. The controller performances were also compared for larger changes in set point, as shown in Figure 3. It was found that, as the set-point changes get larger, the performance of the geometric controller gets worse. In fact, even for a 4% change in set point, the system with

the geometric controller becomes unstable. In contrast, the PI controller works fine for a change of this magnitude. It is to be noted that the maximum set-point change reported in ref 8 is 2.5%. We next study the issue of disturbance rejection. At t ) 5 min, the process at nominal steady state is subjected to a 5% increase in vapor-stream inlet flow rate FAo. It can be seen from the closed-loop profiles in Figure 4 that the performances of the two controllers are more or less similar. It was found that the new controller works well for all other perturbations reported in ref 8 as well. In the interest of space, these results are not reproduced here. The controller proposed by Kumar and Daoutidis8 is computationally involved, as it involves several symbolic differentiations. As a result, there are concerns in making a comparative study based on simulation results, especially through a general purpose language such as Fortran. To overcome this shortcoming, the work of Vora et al.10 is completely reproduced here through Mathematica, a more appropriate environment, especially for implementation of the control technique of Kumar and Daoutidis. Indeed, a minor discrepancy was observed: whereas Vora et al.10 reported that the closed-loop system became unstable for a 6% change in temperature set point, Mathematica reports that the instability occurs at 4%, as shown in Figure 3. Whereas the powerful symbolic computing facilities of Mathematica allowed the geometric controller to be


Figure 3. Closed-loop profiles of the controlled outputs and manipulated inputs for a 4% increase in the set point for the temperature: proposed approach (solid line) and that of ref 8 (dashed line).

developed with relative ease as compared to the Fortran implementation, the former is much slower in performing computations. For example, Mathematica takes 30 min to simulate the system with the geometric controller in place and to plot the dashed lines in Figure 1. The corresponding time required by the Mathematica, geometric controller combination for the simulation results in Figure 3 is 90 min. Whether Mathematica or Fortran, the simulations with the PI controller proceed much more quickly. For example, Mathematica requires less than 1 min to carry out any simulation of the closed-loop system with the PI controller reported in this work. 5. Control of a CSTR with a Heating Jacket The second example we present is a CSTR with a heating jacket. This is an interesting problem because the system includes both fast and slow modes, which makes its control somewhat complicated.3 5.1. Process Description. Reactant A is fed to the CSTR at a flow rate FA, molar concentration CAo, and temperature TA. An irreversible endothermic reaction A f B yields the product B, and the product stream is withdrawn at flow rate F ) FA, i.e., the reactor holdup volume V is constant. The reaction rate RA is given by the expression

activation energy, respectively; T is the reactor temperature; and CA is the molar concentration of A in the reactor. Heat is provided to the reactor from the jacket, where a heating fluid is fed at flow rate Fh and temperature Th. The modeling equations for the process include the mole balances for the two components in the reactor and the energy balances in the reactor and the jacket. 5.2. Process Model. The rate-based model of the process, including an explicit expression for the heattransfer rate, takes the form of the following DAE model3

C˙ A )

FA (C - CA) - koe(-Ea/RT)CA V Ao

C˙ B ) -

FA C + koe(-Ea/RT)CA V B

T˙ )

FA ∆Hr Q (T - T) - koe(-Ea/RT)CA + V A Fcp FVcp

T˙ j )

Fh Q (T - Tj) Vh h FhVhcph

0 ) Q - UA(Tj - T)

RA ) koe(-Ea/RT)CA

y 1 ) CB

where ko and Ea are the reaction rate coefficient and

y2 ) T

(10)


Figure 4. Closed-loop profiles of the controlled outputs and manipulated inputs for a 5% increase in the inlet flow rate FAo: proposed approach (solid line) and that of ref 8 (dashed line).

It is assumed that the density and specific heat capacities of the two liquids are the same, i.e., Fh ) F and cph ) cp, and that the liquid holdup in the jacket at temperature Tj has a constant volume Vh. Furthermore, it is assumed that the heat transfer between the jacket and the reactor, Q ) UA(Tj - T), is fast compared to the reaction, i.e., (UA/Fcp) ) (1/) has a large value. Note that the dynamic model in eq 10 includes the rate expression for the rapid heat transfer Q. In this process, it is desired to control the product concentration CB and reactor temperature T using the reactant flow rate FA and the heating fluid flow rate Fh as the manipulated inputs. Define the differential (x) and algebraic (z) variables and the controlled (y) and manipulated (u) inputs as

x1 ) CA, x2 ) CB, x3 ) T, x4 ) Tj, z ) Q, y1 ) CB, y2 ) T, u1 ) FA, u2 ) Fh The system given by eqs 10 takes the form

1 x˘ ) f(x) + g(x)u + b(x) k(x) yj ) hi(x), i ) 1, 2, ..., m

(11)

where x ) [CA CB T Tj]T, u ) [FA Fh]T, k(x) ) x4 - x3,

and

[

-koe(-Ea/Rx3)x1

f (x) )

koe(-Ea/Rx3)x1

( )

0

]

[ ]

-koe(-Ea/Rx3)x1

∆Hr , Fcp

(CAo - x1) 0 V x2 0 V , b(x) ) g(x) ) (T - x ) A 3 0 V (Th - x4) 0 Vh

[] 0 0 1 V

-

1 Vh

5.3. Linearized Model of a CSTR with a Heating Jacket. We assume that the variables deviate only slightly from the nominal operating point. The nonlinear model in eqs 10 is linearized using a Taylor series expansion around the operating point. Let

wi ) xi - x/i , wj ) zm - z/m, vk ) uk - u/k where i ) 1, ..., 4; j ) 5; m ) 1; k ) 1, 2; x/i , z/j , u/k are equilibrium points; wi, wj represents the new state


Figure 5. Comparison of the closed-loop profiles of the controlled outputs and manipulated inputs for a 15% increase in the product concentration CB in the nominal system: proposed approach (solid line) and that of ref 3 (dashed line).

variables; and v1, v2 are new input variables. The linear model is as follows

w˘ 1 ) -0.939 76w1 - 0.0761w3 + 0.3404v1 w˘ 2 ) 0.639 76w1 - 0.3w2 + 0.0761w3 - 0.3404v1 w˘ 3 ) -3.5544w1 - 0.7227w3 + 2.7778 × 10-5w5 + w˘ 4 ) -0.1w3 - 2.7778 × 10-4w5 - 89.63v2

1.592v1

0 ) w5 - 25 000(w4 - w3) y1 ) w2 y2 ) w3

(12)

system of equations 13 is impulse-controllable and that it has index 1. One can also directly verify that the linear system given by eqs 13 is of index 1. It can be converted to a regular state-space model by substituting the values of algebraic variable w5 into differential equations. From this regular state-space model, we arrive at a 2 × 2 transfer function matrix between v and y, as given in Appendix B. Any standard control design technique can now be used for this system. In this work, we use the ZieglerNichols method to design controllers for individual loops. There is no need to detune these loops. In addition, the reference signal is passed through a lowpass filter, as the controller gains are high. According to this procedure, the PI controller of the following form was designed

This linear singular system can be written as

u1 ) -1.487 91(y10 - y1) -

Ew˘ ) Aw + Bv y ) Cw

(13)

5.4. Controller Design Using the Proposed Approach. From the impulse controllability test

rank

[

]

E 0 0 ) 9 ) n + rank E A E B

it can be seen that the system of equations 13 is impulse-controllable. It can be verified that deg(|sE - A|) ) rank E ) 4, and hence, from lemma 1, it can be concluded that the

u2 ) 72.875(y20 - y2) +

∫

1.487 91 (y10 - y1) dt + 8.2325 3.0

∫

72.875 (y20 - y2) dt + 8.727 × 10-4 0.1 (14)

where y10 and y20 are the outputs of low-pass filters

y10 )

ysp1 ysp2 , y20 ) 10s + 1 15s + 1

where ysp1 is the set point for y1 and ysp2 is the set point


Figure 6. Comparison of the closed-loop profiles of the controlled outputs and manipulated inputs for a 10% increase in the reactor temperature T in the nominal system: proposed approach (solid line) and that of ref 3 (dashed line).

for y2. We next briefly discuss the controller presented by Kumar and Daoutidis.3 5.5. Geometric Controllers. As in the previous example, we compare the performance of the PI controller with those of geometric controllers reported in the literature. It is not straightforward to design a controller for the CSTR through the geometric control approach. Indeed, Kumar and Daoutidis3 tried three different approaches to finally arrive at a control configuration that works. Their three approaches are as follows: Direct Approach. An attempt is made to directly control the system presented in section 5.2 through the global linearization technique. For the rate-based model in eqs 10, the relative orders of the outputs y1 and y2 with respect to the manipulated input vector u1 are r1 ) 1 and r2 ) 1, respectively. However, the other manipulated variable u2 is absent in both. This results in the characteristic matrix being singular.12 Thus, the direct approach is not possible in this case. Dynamic Extension. The above difficulty suggests a dynamic feedback controller, which is designed through dynamic extension by defining v1 ) u˘ 1 as a new manipulated input to obtain a system with an extended state vector xj ) [xT u1]T, new input vector v ) [u˘ 1 u2]T, and relative orders r1 ) 2 and r2 ) 2. Although the characteristic matrix is now nonsingular, more differentiations are required to compute this matrix as the relative orders have increased. An input/output linearization dynamic state feedback controller is designed on the basis of the extended system to induce a wellcharacterized closed-loop response. Nevertheless, the relative orders of the outputs in the extended system

are both 2. Consequently, the resulting controller involves terms multiplied by large factors, (1/) and (1/2), the consequence of which is that the controller is severely ill-conditioned. Extraction of Slow Modes. As both of the above methods fail,3 Kumar and Daoutidis consider the twotime-scale behavior of the model and convert the DAE system presented in section 5.2 into a standard singular perturbation form through nonlinear coordinate transformation. From this formulation, a three-dimensional slow subsystem is extracted. This slow subsystem consists of only ODEs. The controller is designed for this system only, as the fast system is stable and hence is left in open-loop form. We refer to the controllers designed through the above approaches as geom-1, geom-2, and geom-3, respectively. 5.6. Comparison of Controller Performance. In this section, we compare the performance of the PI controller with that of geom-3 (see the previous section for nomenclature). As suggested by Kumar and Daoutidis,3 the performances of the controllers for set-point changes and parameter uncertainties are evaluated. The first two simulation runs address the set-point tracking capabilities and the next two are for parametric uncertainties. In the first run, at time t ) 30 min, the process at its nominal steady state is subjected to a 15% increase in ysp1 or, equivalently, to an increase in the concentration of product B. The corresponding input and output profiles under the two controllers are shown in Figure 5. Clearly, the controller based on the new approach and


Figure 7. Comparison of the closed-loop profiles of the controlled outputs and manipulated inputs for a 15% increase in the product concentration CB in the presence of 0.5% error in ∆Hr, F, and cp: proposed approach (solid line) and that of ref 3 (dashed line).

geom-3 provide good responses for set-point changes in y1 and also a decoupled responses, as required. In the second run, the process is subjected to a 10% increase in ysp2, i.e., an increase in the desired reactor temperature T. Again, both controllers provide a good response for y2 and a decoupled response, as required. The output and manipulated input profiles are shown in Figure 6. In these two runs, the performances of two controllers are more or less similar. Moreover, these results are also in agreement with those of Kumar and Daoutidis.3 The third run is performed for a 15% increase in ysp1 in the presence of small (+0.5%) errors in the process parameters ∆Hr, F, and cp. The corresponding input and output profiles are shown in Figure 7. Clearly, the PI controller yields a performance virtually indistinguishable from that in the nominal case and gives results as good as those of geom-3. In contrast, geom-2 is sensitive to these errors.3 The fourth run is identical to the third but for the following difference: a larger error of +5% is applied to the process parameters. The results are shown in Figure 8. It is seen from this figure that the performance of the PI controller is comparable to that of geom-3. It is interesting to note that geom-2 leads to instability in this case.3 DAEs represent systems of extreme variation in the speeds of different modes. For example, algebraic equations represent dynamics of infinite speed. Thus, controllers designed on the basis of DAE control theory are likely to be capable of handling systems consis-

ting of modes of different speeds well. Possibly as a result, the proposed controller works well for this example. 6. Conclusion It is shown that PID controllers can provide good results for complicated dynamic systems described by DAEs. Linearization and linear DAE control law were used as vehicles to arrive at a formulation suitable for PID control design. The PI controllers used in this work are simple to design and implement. From a computational point of view, they are also quite inexpensive. For the examples of a two-phase reactor and a CSTR with a heating jacket, it is shown that these controllers produce at least as good a performance as given by more advanced geometric controllers. Mathematica11 was used to carry out all of the comparative studies reported in this work. As no single controller can be the best one for all processes, there is room for a variety of controllers. As shown in this work, in fact, considering the simplicity of the design and implementation of PID controllers, the controller designer might often benefit from giving the PID controller a chance before examining more advanced controllers. It should also be noted that the first part of the control design technique suggested in this work is not restricted to PID control design. For example, it is possible to apply some other controller to the transfer function matrix determined the end of the first step of linear DAE control design.


Figure 8. Comparison of the closed-loop profiles of the controlled outputs and manipulated inputs for a 15% increase in the product concentration CB in the presence of 5% error in ∆Hr, F, and cp: proposed approach (solid line) and that of ref 3 (dashed line).

Last, but not least, is the issue of using linearized models for controller design. This approach is usually valid for only small perturbations about the operating point. Quantification of the maximum permissible amount of deviation, through perhaps robust control methodologies applied to nonlinear DAE systems, is an open problem.

The entries for the CSTR with heating jacket are given by

G11 )

-0.3404s3 - 2.8613s2 - 1.7309s - 0.271 s4 + 9.7014s3 + 15.6624s2 + 6.7969s + 0.8833

G12 )

1.592s3 + 11.98s2 + 5.467s + 0.6049 s + 9.7014s3 + 15.6624s2 + 6.7969s + 0.8833

G21 )

1.776 × 10-15s3 + 4.737s + 1.421 s4 + 9.7014s3 + 15.6624s2 + 6.7969s + 0.8833

G22 )

62.24s2 + 77.17s + 17.55 s4 + 9.7014s3 + 15.6624s2 + 6.7969s + 0.8833

A. Description of Matrices Let the transfer function matrix be given by

[

G G G ) G11 G12 21 22

]

The entries for the two-phase reactor are given by G11 ) -6.07 × 10-4s3 + 4.564 × 10-4s2 + 1.398 × 10-5s + 1.528 × 10-8 s4 + 0.070 29s3 + 9.789 × 10-4s2 + 3.445 × 10-6s + 2.426 × 10-9 G12 ) -1.368s3 - 0.049 14s2 - 4.292 × 10-4s + 1.898 × 10-7 s + 0.070 29s3 + 9.789 × 10-4s2 + 3.445 × 10-6s + 2.426 × 10-9 4

G21 ) 1.442 × 10-7s3 - 3.598 × 10-8s2 - 1.183 × 10-10s + 1.239 × 10-13 s4 + 0.070 29s3 + 9.789 × 10-4s2 + 3.445 × 10-6s + 2.426 × 10-9 G22 )

1.849 + 4.765 × 10-6 s2 + 3.419 × 10-8 s + 7.131 × 10-11 s4 + 0.070 29s3 + 9.789s2 + 3.445 × 10-6s + 2.426 × 10-9

4

Acknowledgment The encouragement received from Prof. Bill Luyben at the early stages of the career of K.M.M. is acknowledged. We also acknowledge important advice from Bill. Literature Cited (1) Luyben, W. L. Process Modeling, Simulation and Control for Chemical Engineers; McGraw-Hill: New York, 1990. (2) Brenan, K. E.; Campbell, S. C.; Petzold, L. R. Numerical Solution of Initial-Value Problems in Differential Algebraic Equations; SIAM: Philadelphia, PA, 1996. (3) Kumar, A.; Daoutidis, P. Control of Nonlinear Differential Algebraic Equation Systems; Chapman & Hall/CRC: Boca Raton, FL, 1999. (4) Krishnan, H.; McClamroch, N. H. A new approach to position and contact force regulation in constrained robot systems.

4610 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 In Proceedings of the IEEE International Conference in Robotics and Automation; IEEE Press: Piscataway, NJ, 1990; pp 13441349. (5) Krishnan, H.; McClamroch, N. H. Computation of state realizations for control systems described by a class of linear differential-algebraic equations. Int. J. Control 1992, 51 (4), 1425. (6) Krishnan, H.; McClamroch, N. H. Tracking in nonlinear differential algebraic control systems with applications to constrained robot systems. Automatica 1994, 30, 1885. (7) Byers, R.; Geerts, T.; Mehrmann, V. Descriptor system without controllability at infinity. SIAM J. Control Optimization 1997, 35, 462. (8) Kumar, A.; Daoutidis, P. Feedback control of nonlinear differential-algebraic equation systems. AIChE J. 1995, 41 (1), 619. (9) Dai, L. Lecture Notes in Control and Information Sciences, Singular Control Systems; Springer-Verlag: Berlin, 1988.

(10) Vora, P.; Moudgalya, K. M.; Pani, A. K. Control of a high index DAE system through a linear control law. In American Control Conference; IEEE Press: Piscataway, NJ, 2002, pp 465470. (11) Wolfram, S. The Mathematica Book; Cambridge University Press: Cambridge, U.K., 2002. (12) Daoutidis, P.; Kumar, A. Structural analysis and output feedback control of nonlinear multivariable processes. AIChE J. 1994, 40, 647.

Received for review November 4, 2002 Revised manuscript received April 17, 2003 Accepted April 17, 2003 IE0208671

PID Control of DAE Systems - American Chemical Society

Recommend Documents