Ind. Eng. Chem. Res. 1995,34, 545-556
545
Optimization and Control of Semibatch Reactors Jyh-Shyong Chang* and Wen-Yen Hseih Department of Chemical Engineering, Tatung Institute of Technology, Taipei, Taiwan, ROC
We propose a n integrated method for optimization and control of semibatch reactors. Based on the desired control objective, dynamic programming is applied to obtain optimal operating trajectories. Yield optimization is assured for a real plant by tracking model-dependent optimal trajectories according to the proposed modified globally linearizing control (MGLC) structure. The behavior of the proposed MGLC structure is predictable and reliable, with tuning parameters based on the proposed tuning method if the manipulated variables are not constrained.
Introduction Batch or semibatch processes require control strategies different from those for continuous processes because they are commonly not operated at steady states. From an engineering point of view, the objectives of control of batch or semibatch chemical reactors are divided into two categories: optimization of yield and control of product quality (Kravaris et al., 1989). To achieve the desired control objective, the performance of a proposed control strategy is generally set to force the system output to track some optimal or reference trajectories efficiently. In this way, one can operate a batch or semibatch reactor t o meet the constraints of temperature and t o obtain the product at specification in least duration. Operation in a semibatch mode is typically used for strongly exothermic reactions, because one can balance the reaction heat with the available cooling of the reactor and the adjustment of federate (Bhat et al., 1991). Bhat et al. (1991) addressed temperature control of a semibatch reactor; synthesis of the hexyl monoester of maleic acid was involved. They proposed a control strategy to tackle the problem of coordination of reagent addition and coolant flow. The concentration of the reagent, hexanol, added to the reactor was one controlled variable and the temperature of the reacting mixture another. Reactor control requires typically specifications such as maximum allowable rate of temperature rise or overshoot (Bhat et al., 1991). According to the desired control specifications, a trajectory of temperature set point is specified a priori (Liptak, 1986). Then the resulting set-point policy of concentration of loading reagent was devised (Bhat et al., 1991). By proper tuning of control loops, they reported that such a set-point policy can limit the rate of temperature rise within the permissible value. Nevertheless, the calculated reference trajectory depends on the model and is unpredictable. A similar set-point policy is obtainable directly by the technique of dynamic programming (de Tremblay and Luus, 1989; Luus, 1990a; Luus, 1990b; Luss and Rosen, 1991). The appealing feature is that we need do no prior analyses of characteristics of reaction modes to find a suitable operating path. Furthermore, Bhat et al. (1991) considered no constraint of final total volume of reacting medium, the main drawback when one adopts the proposed approach of Bhat et al. (1991) in a practical setting. According to this discussion, we focus on the semibatch trajectory tracking problem addressed by Bhat et al. (1991) as a test problem according to an approach that we propose.
* Author t o whom correspondence should be addressed.
Optimization of dynamic programming is based on direct search and systematic decrease of the search region. M e r each iteration, the size of the region is diminished so that the optimum is found accurately. This optimization technique preserves the following attractive features in yielding a global optimal operating trajectory for batch or semibatch processes facing constraints of hardware: (a) the optimization routine is executed easily without excessive mathematical manipulation; (b) the computational load is moderate; (c) when the states of the system over the course of operation are chosen as optimal reference trajectories, one can devise a suitable control system to achieve properly the control objective. The way to formulate an optimal control problem and t o solve it by dynamic programming appears in Appendix A (Luus, 1989). Soroush and Kravaris (1993) proposed a framework for integrated design and operation of a single-stage batch or semibatch reactor. A systematic decoupling of optimization and design of reactor dynamics was proposed and applied successfully in an experimental demonstration (Soroush and Kravaris, 1992). One can obtain two subsystems with distinct characteristics. Alternatively, dynamic optimization and design are done in a single step (Soroush and Kravaris, 1993).This approach involves solutions of constrained dynamic optimization based on a complete dynamic model. However, bang-bang and singular control policies may result. One may find problems in implementing optimal control that are attributed to the open-loop form of the control law being sensitive to modeling uncertainties (Bhat et al., 1991). In this work we adopted the singlestep approach, but we used dynamic programming instead of dynamic optimization t o be immune from possible problems discussed above. Kravaris and Chung (1987) developed the GLC structure. They proposed a nonlinear transformation of the form u = 4 x 1 %x)u such that the closed-loop system has linear inputJoutput behavior. The state feedback control law, u = flu,x) can be then obtained by means of model inversion. In application of the inputloutput linearization method, there is invariably provision for state feedback. If the manipulated input appears linearly in the nonlinear model equation, analytic calculation of state feedback exists. However, closedloop stability by means of tuning controller parameters exists for only a minimum phase. Our results (Chang and Huang, 1994) showed that the modified globally linearizing control (MGLC) structure is applied efficiently to track the given trajectories of a batch reactor. Development of the MGLC structure was based on the GLC structure (Kravaris and Chung, 1987; Henson and Seborg, 1990). According to the MGLC structure, the
+
Q888-5885/95/2634-Q545$09.QQIQ 0 1995 American Chemical Society
546 Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995
Table 1. Model Parameters and Given Conditions for a Semibatch Reactor k = 1.37 x 10I2exp(-12628/13, m3/kmol-s -AH/&', = 16.92 m3*K/kmol CAO= 10.1kmol/m3, CBO= 0.0 kmol/m3, To = 328 K, VRO= 2.2 m3, VRf = 4.7 m3 U I , =~0 ~m3/s, ~ Ul,max= 0.01 m3is Uz,,,, = 0 s-', Uz,,,, = 0.000253 s-'
CBl = 9.7 kmol/m3
command variable u is modified to include the derivative terms of the desired trajectory to the relative order r of the process model. Based on frequency domain analysis, a transparent tuning method for the controller parameters was presented (Chang and Huang, 1994). Bhat et al. (1991) proposed a global inpuVoutput-linearized internal model control (GLIMC)structure. In this work, we examine the trajectory tracking performance between the proposed MGLC structure and GLIMC structure. In summary, achieving the desired successful control of semibatch processes hinges on integration of three important ingredients: an optimal operating trajectory, a suitable control law, and a suitable design of the control configuration. The optimal operating trajectory can be obtained based only on an available physical model, but there invariably exists some mismatch between plant and model; thus, the obtained optimal operating trajectory cannot be an exact optimal operating trajectory of a real process. If one designs a suitable control law t o track the obtained reference trajectory efficiently, whether the desired control objective, such as conversion, can be met is questionable. For semibatch reactor systems, the answer is affirmative if the cited three ingredients are obtained and designed properly. Details are reported in the following sections.
Reactor Control Problem. Consider the synthesis of the hexyl monoester of maleic acid according the following reaction:
+ hexanol-.
(A)
(B)
Objective: min[X(t,) - 99%12
+
hexyl monoester of maleic acid water Because of the large heat of reaction (AH = -33.5 MJI kmol), Bhat et al. (1991) addressed this reaction conducted in a semibatch reactor. In semibatch operation, reagent A (maleic anhydride) is first melted; then reagent B (hexanol)is added at a regulated rate so that heat generated is matched by cooling capacity. The mathematical model proposed by Westerterp et al. (1984) and adopted by Bhat et al. (1991) is
Model parameters and conditions of the semibatch reactor appear in Table 1. We simulated the actual process by eqs 1-4 but deviations of the rate coefficient k and parameter for rate of coolant flow Uz from values
(5)
Subject to system eqs 1-4 with constraints:
qmin IUJt) Iqmax, i = 1,2
Application to a Semibatch Reactor
maleic anhydride
in Table 1 were assumed. This viewpoint is practical, because t o design a model-based controller one must identify a suitable physical model. The identified physical model and its parameters are then intended to be adopted even if mismatch between plant and model occurs. Modeling errors may occur due to variation of the rate coefficient of the real process induced by varied content of impurity in the feed. Another possibility results from a decreased coefficient of thermal transfer induced by fouling of the equipment. Although the upper temperature of this reaction is less critical than for other reactions, because of reversibility of the reaction and formation of diester, temperatures exceeding 100 "C should be avoided (Hugo, 1981). In the control system studied, the control objective, that is, the minimum duration of operation to effect 99% conversion of reagent A, should be achieved. The rate of loading of reagent B and rate of flow of coolant are chosen as manipulated variables. The temperature of the reactor and concentration of reagent B in the reactor are chosen as controlled variables. Optimal Operating Trajectories of the Semibatch Reactor. Westerterp et al. (1984) showed that 99% reaction was completed in 8000 s on adopting a constant rate of addition of reagent and full cooling. We considered the problem of minimum duration of operation for 99% reaction by formulating the following constrained optimal control problem:
(6) (7)
v,,
v&"
(8) The formulated problem of constrained optimal control is solved by means of dynamic programming, whereas the minimum duration is obtained by iterative solution of the formulated problem. Optimal operating trajectories in two sets were determined by dividing the interval [O, tfl into 10 and 20 stages respectively. Figures 1 and 2 display the state trajectories, T(t) and CB(t) corresponding to these two cases. The control objectives and related constraints in these two cases were all met satisfactorily in 6570 s. In the work of Bhat et al. (19911,6570 s was also required to achieve the same target with the prescribed trajectories of temperature and composition set point (Westerterp et al., 1984); their other cases require much greater duration. As discussed in the Introduction, our optimal operating trajectories are similar to those based on on-line determination of a set-point policy. However, prior analyses of characteristics of reaction modes are unnecessary. Furthermore, we consider also the constraint of final total volume of reacting medium. In essence, the manipulated variables Ul(t) and Uz(t) depicted in Figures 1 and 2 are expected to be reproduced if these two control actions are not constrained during operation. This predictable performance is attributed t o the nonlinear model inversion method embedded in the proposed MGLC structure. The optimal operating trajectories (Figures 1 and 2) require hrther discussion. If there exists any modeling error or process noise, the obtained optimal trajectories 5 VR(tf)5
Ind. Eng. Chem. Res., Vol. 34,No. 2, 1995 547 0.30
1
I
0.25 0.20 0.15
ZY 0.10 0.05
0.16 0'0° 0.12
1
$?
3 O*OB 0.04 0.00
-u_ !?
$?
370
2.5
360
2.0
0"
350 340
1.0
g
330
0.5 o E
520
0
0.0 0
1000 2000 3000 4000 5000 6000 7000 Time Is)
Figure 1. Optimal control policy and system states CBand T for a semibatch reactor system obtained via dynamic programming; time stages, P = 10. 0.30
II
and 10, if VL and CB(tf)are met at t = tf according to the model-dependent operating trajectories, the desired conversion x(td of the real process is obtained exactly. This conclusion is correct for semibatch reactor systems generally considered even if there exists error in modeling of the real process. Further, if the model-dependent state trajectories T(t)and Cg(t) can be tracked tightly by any control law based on model inversion, then U1( t )and U d t ) can be generally reconstructed. Because the reconstructed feeding rate Ul(t)preserves stepwise variation like that obtained by means of dynamic programming, whether the desired total feeding volume VL is achieved is easily verified. In this way, the attainment of VLand CB(tf)is achieved properly and the desired conversion of plant is ensured. Because tight trajectory tracking of T(t)and C d t ) is desired, constraint control of states of the given problem is met concurrently. Based on this reasoning, the model-dependent optimal trajectories found according t o dynamic programming are really worthy of being tracked, because one can ensure obtaining the desired control target, i.e., the conversion of the real process at the end of a semibatch run. The states, T(t)and Cg(t),obtained from dynamic programming are then used as the optimal set points for the semibatch control system.
Modified Globally Linearizing Control Structure for the Control of MIMO Nonlinear System
I
0.20
0
0.15
The SISO MGLC structure (Chang and Huang, 1994) can be directly extended to a multiinputJmultioutput (MIMO) system. Consider MIMO nonlinear systems with inputs and outputs of equal number of the form
0.10
0.05
0.16 O.OO
rl
m
+ cgj(X)uj
;;;E I:& x = fix)
j=l
(11)
y i = hi(x), i = 1, ..., m
in which
; !,
,s
1.5 .-v)
340
1.0
g
330
0.5
5 V
320
0
0.0 1000 2000 3000 4000 5000 6000 7000 Time
(SI
Figure 2. Optimal control policy and system states CBand T for a semibatch reactor system obtained via dynamic programming; time stages, P = 20.
cannot be the real optimal operating path for a semibatch reactor. However, these two state trajectories are still adopted as optimal operating set points in closedloop control. In operation of a semibatch reactor, the desired final conversion x(tf)is conventionally given as
Given a multivariable nonlinear system of the form of eq 11,which is affne in the inputs, the ith output yi is defined t o have relative order ri if ri is the smallest integer for which
[LglL;-lhi(X)Lg2L;-lhi(X),...,LgmL;-lhi(X)l# [O 0 ... 01 (13) Assume that each output yi possesses its relative order ri. The matrix
1
in which
According to the above equation, one finds that whether the control objective can be achieved depends on two factors: (a) the total volume of component B fed into the reactor during operation VL; (b) the concentration of that component at the end, CB(tf). Based on eqs 9
is called the characteristic matrix of the system of eq 11. Ha and Gilbert (1986) considered the problem of finding a static method of state feedback design. On application of the proposed method to the system of eq 11, input and output of the closed-loop system become decoupled; i.e., the ith output depends on only the ith
548 Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995
external input:
=
(15)
lo b
for some P i k scalar and ri a positive integer. eq 11 can be decoupled under static state feedback if all its outputs possess relative orders and its characteristic matrix is nonsingular for all x. The required static state feedback for decoupled closed-loopinputloutput behavior is given by
(20)
The relative orders of the system considered are = 1 and r2 = 1. because
r1
U =
and
Moreover, the characteristic matrix
are scalar parameters. The inputs u, of MGLC are given by
KCL
u, = KcI(YISP(t) - y,(t))
+ -Jt(r,"p(t)
- y,(t)) dt
+
$ 0 (24
is nonsingular except at the initial conditions. Thus, the 2 x 2 system can be transformed by global linearization to a decoupled linear system. The control law for state feedback is then obtained:
TI1
[Z]
whereas those for GLIMC are
= Y(x,v)=
in whichyLsP is the set point. The computational method of trajectory-based command signals, yisP and dyisPldt, for the MGLC structure is presented in Appendix B.
I.
Simulation Results and Discussion
-kC,C, -kCACB f(x)= -AH
(-)kc
,o
(19) ACB d
in which the command signal ui are chosen according to eq 17 or 18 to correspond to the MGLC or GLIMC structure. In general, P I 0 and PZOare chosen to be unity, eq 25 becomes expressed explicitly as
Ind. Eng. Chem. Res., Vol. 34,No. 2, 1995 549 I
Table 2. Simulation Conditions for Control Tests of a Semibatch Reactor
I
I
U 'I
f
dk
Ea,7
YI
State
I
y2
Y,
Figure 3. Modified globally linearizing control (MGLC)structure for a MIMO system.
u2 = - -
V2 - T
[
P21
+
"(
T - 327 vR
1-
($)kcAcB]T
modeling error
yq
figure no. 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
UZ,,,
rate of reaction
x 1.5
no no no no no no no Yes no Yes no yes yes no yes
no no no no no no Yes Yes no no Yes Yes Yes
heat transfer coeff (+20%)* no no no no no no no no Yes Yes no no no
C
C
C
C
measurement noise T 0.5 "C no no no no no Yes no no no no Yes Yes Yes Yes Yes
The rate coefficient of the process is assumed to exceed that in the model. The coefficient of thermal transfer of the process is assumed to be smaller than that in the model. Process mean noise w = 0.028 is added to the process, eq 3. a
- 327 (27)
To execute the state feedback of the MGLC or GLIMC structure, all process state variables should be mea-
Table 3. Tuning Parameters for Control Tests of a Semibatch Reactor
sured or estimated continually. Continual measurement of all states is commonly impracticable. In such cases it is necessary to estimate unavailable states from a dynamic model and output measurement. As shown in the state feedback control law, eqs 26 and 27, u depends on four states CA, CB,T, and VR. Among these states, CAand VR cannot be measured continually. A state observer is required to estimate CAW and V&). Based on eqs 1 and 4,this process involves continual integration of the following differential equations:
figure no. PIO PII 4 1 1 5 1 3600 6 1 1 7 1 100 8 1 3 6 1 3600 9 10 1 36 11 1 36 12 1 36 13 1 36 14 100 3600 15 1 36 16 100 3600 17 1 36 18 1 36
in which CA, and VR, denote respectively estimates of CAand VRusing available measurements of T and
CB. A block diagram of the process and controller designed based on the MGLC structure is depicted in Figure 3. This figure shows the various blocks of the controller and process and their interconnections. If the sampling period is 1s, the actual duration (CPU time) needed to execute computation tasks (PI controller, state observer, state feedback) is 0.05 s. A microcomputer (IBM PC/ AT type) was used. This small CPU time on a modest microcomputer shows the computational efficiency of the nonlinear control method. Therefore, applicability of the proposed method t o active control in an actual setting is acceptable. For the following tests, simulation conditions and tuning parameters appear in Tables 2 and 3, respectively. For all tests shown in Table 3, P control was sufficient to tackle the problem of modeling errors occurring in the tests. Therefore, the integral control was almost closed by setting the value of TI to be lo6. The performance of the designed control structure was set to be tight trajectory tracking and measurement noise attenuation. Theoretically, one can obtain optimal tuning parameters for either the MGLC or GLIMC structure by analyzing their linear input/output open-
tuning parameters Pzo
PZI
1 1 1 3600 1 1 1 100 1 3 6 1 3600 1 36 1 36 1 36 1 36 100 3600 1 36 100 3600 1 36 1 36
Ke1
TII
KCz
512
z
1 x 10-6 106 1 x 10-6 106 0.990 106 1 1 106 0.991 0.990 0.988 1 106 1 lo6 0.987 106 1 1 106 0.991 106 1 1 lo6 0.995 1 106 1 lo6 0.994 106 1 1 lo6 0.993 1 106 1 106 0.990 1 106 1 lo6 0.995 106 1 1 lo6 0.994 1 106 1 lo6 0.994 106 1 1 106 0.990 1 106 1 106 0.990
loop transfer functions (Figures 19 an1 20). Furthermore, if the control action is not saturated during operation, tuning parameters based on the method in Appendix C can be applied to the real process to effect the control structure satisfactorily. Based on analysis in Appendix C, the property of perfect tracking is inherently preserved in the MGLC but not the GLIMC structure. Therefore, if tight trajectory tracking is desirable in batch or semibatch control, we recommend the MGLC structure. As noise attenuation is improved in step 2 of the proposed tuning method (Appendix C), we conclude that sensitivities of U I and UZ are also diminished with respect to noise-corrupted measurements. Furthermore, if variation of control action is t o be minimized with respect t o the effect of noise-corrupted measurement, the sensitivity of each term in the control law given by eqs 26 and 27 must be analyzed individually. For the MGLC structure, the sensitivity of the first term in the right side of eq 27 with respect to noise decreased with increased values of K,z andI%!, maintaining the ratio Kc2/t'321. Hence, variation of the control action UZwith respect to noise was diminished but the function of the P or PI control embedded in the control law was preserved. However, sensitivity analysis of this kind may be done case by case. Therefore the tuning method (Appendix C) may be adequate for application.
550
Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995 0.30
0.30
0.25 FI
0.20
0.20
E 0.15
E
5. 0.10
5-0.10
0.15
0.05
0.05
I
0.00 0.16
0.16
--
"b 5
1
0.12 0.08 0.04 0.00
2.5
370
-Y
360
? $'
350 340
2.0
I
C
1.5:.
.-
1.0
330
0.5
320
0.0 0
0
1000 2000 3000 4000 5000 6000 7 Time
g
E
0.20 0.15
-
? $
350 340
u" C
1.5 ; .
.-
1.0
2
0.5 o E 0
0.0
(SI
(SI Figure 6. Tracking of trajectories CB and T in Figure 1 via a GLIMC structure: for states, (-1 desired trajectories, (. simulated responses; for manipulated variables, (--). Simulation conditions and tuning parameters are given in Tables 2 and 3. Time
0.30 I
11
2.0
0
a*)
-
2.5
360
330
Figure 4. Tracking of trajectories CB and T in Figure 1 via a MGLC structure: for states, (-1 desired trajectories, (. simulated responses; for manipulated variables, (-1. Simulation conditions and tuning parameters are given in Tables 2 and 3. 0.25
370
-Y-
e)
1
0.30
I
I
0.16
-/
I
v ,
c
0.00
-
0.16
1
r"
I
0.12
0.04
0.04
0.00
0.00
370
2.5
-Y-
360
2.0
E" E
350 340
1.5: . 0.5 o E
330
320
370
u"
-
360
E" $
350 340
1.0
330
0.5 o E
Y-
C
1.5
0.0 1000 2000 3000 4000 5000 6000 7000
320
g
0
0
0
.-2 v1
0
0.0
1000 2000 3000 4000 5000 6000 7000 Time is)
Time (s)
Figure 5. Tracking of trajectories CB and T in Figure 1 via a
Figure 7. Tracking of trajectories CB and T in Figure 1 via a GLIMC structure: for states, (-1 desired trajectories, (. * simu-
MGLC structure: for states, (-) desired trajectories, (. * simulated responses; for manipulated variables, (-1. Simulation conditions and tuning parameters are given in Tables 2 and 3.
lated responses; for manipulated variables, (-). Simulation conditions and tuning parameters are given in Tables 2 and 3.
Based on the preceding discussion, the following simulation results are all predictable and reliable. The results shown in Figures 4-8 show this important feature. In Figures 4-8, no modeling errors, measurement noise, and possible constraints on the manipulated variables throughout operation were assumed. The performances of tracking trajectories of C B ~ P ( ~and ) TBP( t ) shown in Figures 4 and 5 were observed despite varied tuning parameters ( / 3 1 1 , Pzl) being chosen in the MGLC structure. In essence, Kcl can be chosen as 1 instead of 1 x in Figure 4 to give a similar performance of trajectory tracking. However, the control action is wilder than that shown in Figure 4, because a ratio Kc21p21 = 1 weights too much the P control action in the control law (eq 27). As shown in Figure 6, the GLIMC structure gives a performance
equivalent to that of the MGLC structure on choosing (/311, p 2 1 ) = (1,1). However, if one chooses (/311, Pzl) as (100, loo), the trajectory tracking performance of the GLIMC structure depicted in Figure 7 becomes degraded, because the magnitude ratio (MR)of eq C3 cannot invariably preserve the value unity like that of the MGLC structure when the tuning parameters are located in a wider domain. In Figures 4-6, the control action Ul(t)was reconstructed almost by the value given in Figure 1; the desired was also attained. Therefore, the desired conversion of the reaction was achieved. Another feature is that control actions shown in this work are much more predictable than those based on continual set-point determination proposed by Bhat et al. (1991). Therefore, it is much more reliable to apply the proposed integrated optimization and
e)
e)
Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995 551
*I
E
0.30
0.30
0.25
0.25
0.20
0.20
p
0.15
>" 0.10
z?
0.15 0.10
0.05
0.05
0.00
0.00 0.16
i
i
n
1
0.04
0.00 370
1
w4
2.5 2.0
8
1.5
.$
?
$
f 340 E
330 320
-Y
0.5 o 0
0.0 1000 2000 3000 4000 5000 6000 7000 Time
,.....
370
2.5 2.0
8
350 340
1.0
g
330
0.5
320
0.0 1000 2000 3000 4000 5000 6000 7000
360
0
5
u
Time (s)
(s)
Figure 8. Tracking of trajectories CB and T in Figure 2 via a MGLC structure: for states, (-) desired trajectories, (. .) simulated responses; for manipulated variables, (--). Simulation conditions and tuning parameters are given in Tables 2 and 3.
Figure 10. Tracking of trajectories CB and T in Figure 1 via a MGLC structure: for states, (--) desired trajectories, (. simulated responses; for manipulated variables: (-). Simulation conditions and tuning parameters are given in Tables 2 and 3. e.)
The effects of corrupting noise in the measured temperature on performance of the MGLC structure were examined as follows. We simulated the measured temperature Tr according to the equation (Bhat et al., 1991)
+
Tr = T ak - 0 . 8 6 6 ~ ~ ~ ~ (30) in which {ah} is a random noise with zero mean and 0.5 variance. The noisy temperature was filtered to give a filtered temperature TF:
0.00 0.16
-
"b 5
1
I
0.12
TF(k
0.08 0.04 0.00
-u_
370 360
c, 350 E 340 330 320
0
1000 2000 3000 4000 5000 6000 7000'Time
(sl
Figure 9. Tracking of trajectories CB and T in Figure 1 via a MGLC structure: for states, (-1 desired trajectories, (. * simulated e)
responses& for manipulated variables, (-1. Simulation conditions and tuning parameters are given in Tables 2 and 3.
control strategy to operate the reactor system. The trajectory tracking performance depicted in Figure 7 may not achieve the desired conversion, because Ul(t) is not well reconstructed according to the value given in Figure 1. In essence, whether the desired conversion is achieved depends on sensitivity of characteristics of the process considered. Figure 8 depicts the performance of the MGLC structure in tracking the trajectories given in Figure 2 for which the total duration is divided into 20 stages instead of 10 stages; the result is also acceptable. This test proves the feasibility of the optimization technique based on dynamic programming. We use the tuning parameters (/311, p21> = (36, 36) to show that operators have a wide domain to choose the tuning parameters t o operate a semibatch reactor system (2.2 m3).
+ 1)= 0.8TF(k)+ 0.2Tr(k)
(31)
In Figure 9, the temperature response in the presence of noise appears acceptable for the MGLC structure. The tuning parameters ($11, /321) = (3600,3600)were chosen based on step 2 of the proposed tuning method t o give a least erratic performance of the control action, Uz(t). In the following tests, Figures 10-13, no measurement noise was assumed. Sensitivity assessment with respect to uncertainty in the kinetics of the reaction is considered for the model-based MGLC structure. The rate of reaction of the real process is assumed to be onetenth greater than the value used in the model. Figures 10 and 11present results of closed-loop simulation for the MGLC structure. The temperature constraint is violated (Figure 10) because the rate of reaction of the real process is one-tenth greater than that assumed in the model over the course of operation. Also, UZ,, = 0.000 253 s-l is assumed. Even if Udt) maintains its maximum value when t E [ZOOO, 35001 s, the resource of thermal transfer is still limited to provide sufficient cooling. In practice, one can improve the efficiency of thermal transfer t o increase the heat transfer resource by half (Baker and Walter, 1979). The trajectory tracking performance of the MGLC structure was simulated for such a case test (Figure 11). The desired control objective x ( t f ) was achieved suitably (Table 3) even if there existed model mismatch as described above. The result is predictable according to previous reasoning. Increasing the total duration of operation is an alternate means to overcome the limited cooling resource, as the intensity of released heat by reaction becomes diminished.
552 Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995 0.40 I
0.05 0.00
-
0.16
r
I
I
P
0.00
I
0.16
t
I
1
I
0.12
!;!E 1:
0.04
0.04
0.00
0.00
;
1.5 .-,g
-Y
370
2.5
380
2.0
u"
340
1.0
g
330
0.5 o E
n 350
z
M
340
1.0
330
0.5 o
E
V
320 0
320
0.0 1000 2000 3000 4000 5000 6000 7000
0.0 0
Time (si
1000 2 0 0 0 3000 4000 5 0 0 0 6000 7000
Time
Figure 11. Tracking of trajectories CB and T in Figure 1 via a MGLC structure: for states, (-1 desired trajectories, (. * simulated responses; for manipulated variables, (-1. Simulation conditions and tuning parameters are given in Tables 2 and 3. a)
(s)
Figure 13. Tracking of trajectories CB and T in Figure 1 via a MGLC structure: for states, (-) desired trajectories, (. *.) simulated responses; for manipulated variables, (-). Simulation conditions and tuning parameters are given in Tables 2 and 3.
0.30 I
0.30 0.25 cI
0.20
0
0.15
5 0.10 0.05
0.00 0.16
1
1
KT E 0.16
c,
-
"b 7
5
0.12 0.08 0.04
0.00 370
Y
E"
380
0.04
0.00 370
......
2.5
350 340
C
Y
M
E
2.0
J
1.5
.g .-
1.0
g E
2.5
360
0.5 a
330
0.0 IO
320
1000 2000 3000 4 0 0 0 5000 6000 7
J C
320
0
2.0
1.5 ; .
._
350 340
330
0
,..-...
1.0
g Y)
0.5 o E 0
0
1000 2000 3000 4000 5000 6 0 0 0 7
0.0 IO
Figure 12. Tracking of trajectories CB and T in Figure 1 via a MGLC structure: for states, (-1 desired trajectories, (. * simulated responses; for manipulated variables, (-1. Simulation conditions and tuning parameters are given in Tables 2 and 3.
Time (SI Figure 14. Tracking of trajectories CB and T in Figure 1 via a MGLC structure: for states, (-1 desired trajectories, (. simulated responses; for manipulated variables, (-1. Simulation conditions and tuning parameters are given in Tables 2 and 3.
The effect of uncertainty in the coefficient for thermal transfer was also examined. This coefficient of the process was assumed t o be one-fifth less than that in the model. Therefore, a value of UZoverestimated by one-fifth was assumed in the control law. The control action U2(t)is saturated €or most operation (Figure 12) to induce violation of the temperature constraint. Similarly, if the value of uzmaxis increased by half, a perfect trajectory tracking performance is achieved (Figure 13); the reasoning is similar t o that for Figure 10 or 11. To prove the robustness of the proposed MGLC structure when encountering both modeling error and measurement noise, we present Figures 14-16. The rate of reaction was assumed decreased one-tenth in the control law; the temperature measurement was concurrently corrupted. As depicted in Figure 14, a response
was obtained similar to that shown in Figure 10. If the cooling resource is improved, the control performances (Figures 15 and 16) depend on the ratio Kci/Dilli=l,a. Based on step 2 of the proposed tuning and the discussion of eqs 26 and 27, the results are predictable because the control action Uz(t)(Figure 16) is less wild than that shown in Figure 15 with increased values of Kci and ,& although maintaining the ratio Kci/pi1. To test further the robustness of the proposed MGLC structure, we took other process noise into account by adding a random disturbance with standard deviation w = 0.028 t o the process dynamics (eq 3). Noise of temperature measurement was also added. Figures 17 and 18 depict the simulated results. In Figure 17, the cooling resource is limited such that perfect tracking cannot be achieved completely during operation. How-
Time
Is) a)
e.)
Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995 563 n
L"
0.30 I
1.7"
*)
p
0.35 0.30 0.25 0.20
0.25 *)
p
>" 0.15
0.20
0.15
>" 0.10
0.10 0.05 0.00
0.05 0.00
0.16
1
0.16
k
-
'b 0.12 r 5
-
0.04
0.08 0.04
0.00 !.5
!-
370
!.O
8
.5
.g .-
I\
,,w 8 2.5 2.0
C
C
1.5
H
.o
1.0
1.5 o E
0.5 o E
0
0
1000 2000 3000 4000 5000 6000 ; Time
p
0
1.0 3
0.0 0
Is1
1000 2000 3000 4000 5000 6000 i
0
(SI Figure 17. Tracking of trajectories CB and T in Figure 1 via a MGLC structure: for states, (-1 desired trajectories, (. simulated responses; for manipulated variables, (-). Simulati conditions and tuning parameters are given in Tables 2 and B Time
Figure 15. Tracking of trajectories CB and T in Figure 1 via a MGLC structure: for states, (-) desired trajectories, (. .) simulated responses; for manipulated variables, (-). Simulation conditions and tuning parameters are given in Tables 2 and 3.
*)
.g .-U
.e)
n An -.-"
0.40 0.35 0.30 0.25 0.20
0.35 0.30 0
p
>" 0.15
0.25
0.20
>" 0.15
0.10 0.00
0.10 0.05 0.00
0.16
0.16
0.05
1
:::E
0.04
0.00
-
370
2.5
360
2.0
Y,
E f
I
1.5, 350 340
1 .o
!.5
;
u
!.O 8 .5
.o
340
330
0.5
330
390
0.0
320
Time
(5)
a)
ever, perfect trajectory tracking is obtained if the cooling resource is sufficient (Figure 18). The control behavior was again demonstrated to be predictable and reliable in such a noisy environment. Conclusions We propose an integrated framework of optimization and control to operate semibatch reactors. Based on the identified physical model, model-dependent optimal operating trajectories were determined effectively by formulating the optimization problem as one of openloop constrained optimal control. The use of penalty functions t o manage state constraints is incorporated into the iterative dynamic programming algorithm by employing accessible grids and region contraction. The
g
0
1000 2000 3000 4000 5000 6000 7 Time
responses; for manipulated variables, (-). Simulation conditions and tuning parameters are given in Tables 2 and 3.
U
1.5 o E 0
Figure 16. Tracking of trajectories C g and T in Figure 1 via a MGLC structure: for states, (-1 desired trajectories, (. simulated
c
.$ .-
1.0 1
(s)
Figure 18. Tracking of trajectories CB and T in Figure 1 via a MGLC structure: for states, (-1 desired trajectories, (. *) simulated responses; for manipulated variables, (-1. Simulation conditions and tuning parameters are given in Tables 2 and 3.
obtained optimal reference trajectories are more reliable and predictable than those found from the set-point policy proposed by Bhat et al. (1991). Development of the MGLC structure based on a state feedback linearization technique proved to preserve the capability of perfect tracking if the manipulated variables were not constrained throughout operation. Besides a perfect trajectory tracking property, the MGLC structure provides a domain of tuning parameters to perform trajectory tracking in a noisy environment wider than that for the GLIMC structure. If resources of manipulated variables are not constrained for an unacceptable duration, the behavior of the MGLC structure in trajectory tracking is predictable and reliable according to this tuning method. The robust
554 Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995
~
State Observe
Y
I - v
Figure 19. Modified globally linearizing control (MGLC)structure for a SISO system.
I';1
.+
Model
Figure 20. Global inpuffoutput-linearized internal model control (GLIMC) structure for a SISO system.
I
Figure 21. Open-loop transfer function block diagram for MGLC and GLIMC structure.
performance of the MGLC structure was demonstrated through simulations. In operating semibatch reactors, one generally considers a yield optimization problem. In one semibatch reactor system, minimum duration with 99% conversion of reagent A is desired. If the system is provided with sufficient resources for control actions, then modeldependent optimal trajectories can be tightly tracked. According to the analysis and simulation results, we demonstrate that achievement of the desired control objective is ensured despite mismatch between the identified model and the real process.
Acknowledgment We thank the National Science Council (NSC 820402-E-036-067)for support and Dr. T. S. Lin, President of Tatung Institute of Technology, Taipei, Taiwan, ROC. Nomenclature = random noise signal C(x)= characteristic matrix of the system, eq 11 CA,Cg = concentration of reagent A and B in the reactor,
ak
kmoVm3 CBL= concentration of reagent B in the feed stream, kmoV
m3 C, = heat capacity of the reaction mixture, kJ/(kgK) f, g = vector fields that characterize the state model fi' = function of state estimator, eqs 28 and 29 H(s) = transfer function h = state/output map I = performance index, eq A3
I, = Lagrange polynomial J
=
augmented performance index, eq A5
K, = proportional gain k = rate coefficient, m3/(kmol.s) L&h) = kth order Lie derivative of h with respect to f N = number of internal interpolation points n = number of state variable; zero mean random disturbance n' = output variable of zero mean random disturbance, n filtered by the process P = number of stages p i = penalty function for a particular constraint p = penalty function introduced to augment the performance index p ~ - 1= node polynomial R = set of real numbers R, = state constraint over the course of operation r = relative order; number of constrained state variables S , = final time constraint t = time, s T = reactor temperature, "C or K TF(k) = filtered value of T at kth sampling instant T,(k) = temperature of noisy reactor at kth sampling instant Ui = reagent addition rate, m3/s UZ= parameter for rate of flow of coolant, l/s V, = overall heat transfer coefficient of reactor jacket, kJ/ (m2-K-s) u , u = manipulated variables VL = total volume of component loaded over the course of operation, m3 VR = reactor volume, m3 u , v = vector of transformed control variables W = number of constrained state variables w = mean process noise x = conversion of reagent A, eq 9 x = n-dimensional state vector y = output variable ylSp= set point of ith output variable y' = internal state
Greek Letters a = lower bound on the control p = upper bound on the control pL= tunable parameters of inpuuoutput linearized system ynsct,= autocovariance function of n' (-AH) = heat of reaction, MJkmol 6 = penalty function factor K, ; I = transformed vector fields Q = density of reaction mixture, kg/m3 T I = reset time, s q5 = objective function \v = transformation, eq 25 o = frequency, Hz Set-Theoretic Symbol E = belongs to Subscript sp = set point Subscripts 0 = initial f = final L = loading m = model; number of inputs or outputs max = maximum min = minimum
Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995 555
Appendix A. Formulation of the Optimal Control Problem to be Solved by Dynamic Programming A continuous dynamic system is described by means of the differential equation d d d t = f(x,u) (A1) in which the initial states x(0)is given; x is an ( n x 1) state vector and u is an ( m x 1) control vector bounded by
aIuIp The problem of optimal control is to find the control policy u(t)in the interval 0 It Itf that minimizes the performance index
I[x(O),t,I = L:qXx,u) dt
043)
At final time tf, constraints are placed on r of the state variables, for which r In, as either equality constraints = Si,i = 1, 2, ...,?or inequality constraints Xi(tf)
+
J = I[~(O),tfl p(x(tf)) (A61 in which p(x(tf))is a penalty function introduced to ensure that the restriction in eq A4 or A5 is satisfied. When equality constraints specified by eq A4 are imposed, the penalty function is chosen as r
i=l
When inequality constraints specified by eq A5 are imposed, the penalty function
is chosen in which r
(-49)
i=l
If there are state variables to be constrained as shown in eq A10, the penalty function (eq A l l ) should be imposed.
x,(t) I R,, i = 1, 2, ..., r, 0 I t I t,
Appendix B. Computational Method of Trajectory-BasedCommand Signals for the MGLC Structure The required signals for the MGLC structure, ,&(dkysP/ dtk),k = 0, ..., r, are obtained according to two methods: (a) If the given trajectory is expressed in analytic form, then the command signals /3k(dky5P/dtk) are obtained directly. (b) If the trajectory is t o be obtained numerically by applying Pontryagin's principle or other optimization technique, the guessed trajectory, ysP(t)of each subsection of the duration of the batch or semibatch reactor may be expressed as a Lagrange interpolation polynomial as follows: NC 1
(A41
xi@,)ISi,i = 1, 2, ..., ?(A5) To satisfy the constraints specified at final time tf, the augmented performance index was constructed
P M t f ) )= > i ( X ( t / ) )
For a detailed computational algorithm to execute dynamic programming, the reader is referred to the literature (Luus, 1989).
(A10)
in which W P
The total duration of the process is divided into P stages. Combination of the penalty functions is expected t o depend on characteristics of the individual test. A positive constant Bi is clearly t o be chosen such that during an iterative procedure a correct answer t o the original problem is obtained when the augmented performance index J is to be minimized.
i=l
The building blocks are N degree N
+ 1 polynomials Zi(t) all of
p ~ + l ( t=) ( t - tl)(t - t 2 ) ...(t - t ~ + 1 is ) a polynomial of degree N 1with leading coefficient 1, called the node polynomial, whereas Z&) are Lagrange polynomials. ti are N 1 interpolation points and ysP(ti) of eq B1 are 1 ordinates at these points (Villadsen and the N Michelsen, 1978). p ( l ) ~ +is1the first derivative OfPN+l. In applying the orthogonal collocation method, the interpolation points ti are chosen as the zeros of any Jacobi polynomial. The ordinates of optimal trajectory ysPP(ti)of each subsection of the duration of the batch or semibatch reactor are determined by the chosen optimization technique. These data are then stored for computation of command signals /3k(dkysp(t)/dtk) by the interpolation technique of the orthogonal collocation method.
+
+ +
Appendix C. Development of a Tuning Method Based on Frequency Domain Analysis for SISO MGLC and GLIMC Structure If the manipulated variables are not constrained during operation, one can obtain the frequency domain analyses for SISO MGLC and GLIMC control structures. An equivalent open-loop transfer function of their individual closed-loop block diagram (Figures 19 and 20) is depicted in Figure 21. If the condition for internal stability is met by testing the zero dynamics of the openloop system, eq A l , and if the manipulated variable u is unconstrained, then the input/output response between ys%) and y'(t) and that between n(t) and n'(t) becomes well behaved by means of the transfer functions, Hl(s) and Hz(s), respectively. For the MGLC structure
556 Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995
For the GLIMC structure
H l ( s )=
1
pp'
H ~ ( s=)
+ p r - p + ... + p1s + Po 1
+
+
+
((23)
((24)
pp' p 7 - p ... + pls po if the transfer function of the control law of GLIMC is chosen (Bhat et al., 1991) pp'
+ p&-l + pJ-l
+ ... + pls + po + ... + pls + po
((25) pp' For any single-inputhingle-output trajectory tracking problem satisfying these assumptions, the input'output responses from MGLC and GLIMC structures are governed exactly by eqs Cl-C4. Therefore, the magnitude ratio of IHl(iw)l = 1 for the MGLC structures and the optimal tuning parameters K,, TI, and P k are obtained based on suitable response analyses of eqs C1 and C2. However, such a characteristic is not embedded in the design of the GLIMC structure. Although y Y t ) generally varies with time, and the spectral density of its energetic power is primarily located about the low-frequencyregion (o:O-1 Hz), the measurement noise possesses a full spectrum of spectral density energetic power. The following tuning rules are proposed for any control structure based on global linearizing method. Although it is here proposed for a system with linearizability index r = 1,extension of the proposed tuning method to r > 1 is simple. Step 1. Choose PO = 1 and close I control by setting a large value of TI. Step 2. Depending on the measurement noise t o be attenuated and to what extent, a suitable value of the ratio of K,/p1 is chosen. For example, if K, = 1 and p1 = 10, then KJpl = 1/10. One can attenuate a random noise (standard deviation AT = 5 "C) to corrupt the state of temperature about f0.8 "C. These data can be estimated using the equation of the autocovariance function, ynytlllz = ATKC/(2P1(P0 Kc))1/2,which is derived based on eq C2 (Wu and Pandit, 1979). Step 3. Using chosen values of K, and PI,plot the Bode diagram of Hl(s) to test whether values of MR approach unity about the frequency domain ( ~ 0 - 1Hz). When values of MR are not near unity, increase the numerical values of K, and in the same proportion to maintain the ratio KJPl t o be the same as that determined in step 2. Then, replot the Bode diagram of Hl(s) t o discover any improvement of values of MR. If this procedure succeeds, then a set of optimal tuning parameters is obtained; otherwise, return to step 1,but open the action of I control by decreasing the value of TI. For the control scheme of the MGLC structure, the values of MR equal unity in the frequency region considered. Hence, the tracking property is almost perfect (eq Cl), and only step 2 of the above tuning rule is required. For the control scheme of the GLIMC structure, further trial and error may be needed to obtain optimal tuning parameters, because eq C3 fails t o preserve the same tracking property as that inherent in the MGLC structure. Accordingly, we conclude that for the same trajectory tracking problem, the optimal tuning parameters of the GLIMC structure are invariably more limited than those for the MGLC structure.
HGLIMC(S)=
+
Our main contribution is the proposal of the MGLC structure to tackle nonlinear batch or semibatch reactor systems by introducing command signals, CPk(dkysp/dtk). For the trajectory tracking problem, if input and output of a control system can be made linear, then one can determine optimal tuning parameters via analyses of transfer functions Hl(s) and H2(s). Thus we obtain desirable tuning parameters without resort to active trial-and-error tuning if manipulated variables are not constrained during operation.
Literature Cited Baker, C. K.; Walter, G. H. Application of Dynamic Modeling Technique to the Analysis and Prediction of Heat Transfer with Particular Reference to Agitated Vessel. Heat Transfer Eng. 1979,1,28-40. Chang, J.S.;Huang, K. L. Performance Study of Control Strategies for Trajectory Tracking Problem of Batch Reactors. Can. J . Chem. Eng. 1994,72,906-919. Bhat, J.; Madhavan, K. P.; Chidambaram, M. Multivariable Global InputiOutput Linearized Internal Model Control of a Semibatch Reactor. Ind. Eng. Chem. Res. 1991,30,1541-1547. De Tremblay, M.; Luus, R. Optimization of Non-Steady-State Operation of Reactor. Can. J . Chem. Eng. 1989,67,494-502. Ha, I. J.; Gilbert, E. G. A Complete Characterization of Decoupling Control Law for a General Class of Nonlinear System. ZEEE Trans. Autom. Control 1986,31, 823-829. Henson, M.;Seborg, D. Input-Output Linearization of General Nonlinear Process. AIChE J . 1990,36, 1753-1757. Hugo, P. Start-up and Operation of Exothermic Batch Process. Ger. Chem. Eng. 1981,4,161-173. Kravaris, C.; Chung, C. B. Nonlinear State Feedback Synthesis by Global Input/Output Linearization. AIChE J . 1987,4,592603. Kravaris, C.; Wright, R. A.; Carrier, J. F. Nonlinear Controllers for Trajectory Tracking in Batch Processes. Comput. Chem. Eng. 1989,13,73-82. Liptak, B. G. Controlling and Optimizing Chemical Reactors. Chem. Eng. 1986,26,69-81. Luus, R. Optimal Control by Dynamic Programming using Systematic Reduction in Grid Size. Int. J . Control 1990a,51,9951013. Luus, R. Application of Dynamic Programming to High-dimension Non-linear Optimal Control Problem. Int. J . Control 199Ob,52, 239-250. Luus, R.; Rosen, 0. Application of Dynamic programming to Final State Constrained Optimal Control Problems. Ind. Eng. Chem. Res. 1991,30,1525-1530. Soroush, M.; Kravaris, C. Nonlinear Control of a Batch Polymerization Reactor: an Experimental Study. AIChE J . 1992,38, 1419-1428. Soroush, M.; Kravaris, C. Optimal Design and Operation of Batch Reactors. 1.Theoretical Framework. Znd. Eng. Chem. Res. 1993, 32, 866-881. Villadsen, J.; Michelsen, M. L. Solution of Differential Equation Models by Polynomial Approximate; Prentice-Hall: Englewood Cliffs, NJ, 1978; Chapter 3. Westerterp, K. R.; Van Swaaij; W. P. M.; Beenackers A. A. C. M. Chemical Reactor Design and Operation; Wiley: New York, 1984; pp 278-281. Wu, S. M.; Pandit, S. M. Time Series and System Analysis Modeling and Application; Michigan Technology University: Houghton, MI, 1979; Chapter 6.
Received for review J a n u a r y 28, 1994 Revised manuscript received October 7, 1994 Accepted October 20, 1994@ IE9400614 Abstract published i n Advance ACS Abstracts, December 15, 1994. @
