Optimal Reference Trajectory Shaping and Robust Gain-Scheduling

Aug 31, 2009 - Informatica e Sistemistica, UniVersita` “Federico II” di Napoli, Via Claudio 21, 80125 Napoli, Italy, and. Dipartimento d'Ingegneri...
0 downloads 0 Views 541KB Size
9128

Ind. Eng. Chem. Res. 2009, 48, 9128–9140

Optimal Reference Trajectory Shaping and Robust Gain-Scheduling for Transition Control of Nonlinear Processes Pietro Altimari,*,† Lucia Russo,‡ Erasmo Mancusi,§ Mario di Bernardo,| and Silvestro Crescitelli⊥ Dipartimento di Ingegneria Chimica e Alimentare, UniVersita` di Salerno, Via Ponte Don Melillo, 84084, Fisciano (SA), Italy, Istituto sulle ricerche sulla combustione, CNR, Piazzale Tecchio 80, Napoli, 80125, Italy, Facolta` di Ingegneria, UniVersita` del Sannio, Piazza Roma, 82100, BeneVento, Italy, Dipartimento di Informatica e Sistemistica, UniVersita` “Federico II” di Napoli, Via Claudio 21, 80125 Napoli, Italy, and Dipartimento d’Ingegneria Chimica, UniVersita` “Federico II” di Napoli, Piazzale Tecchio 80, 80125 Napoli, Italy

A novel method for transition control of nonlinear processes is presented. Gain scheduling is implemented to ensure stability and desired output behavior over the operating region of interest while transition between the initial and the final steady state is achieved by predictive reference control (or reference governor). In this framework, the search of a feasible reference input sequence is performed so that the closed-loop system moves along an optimal curve of steady states. This curve is constructed off-line so that its points lie far from the border of the process constraint set and correspond to satisfactory controllability properties of the uncontrolled plant. Moreover, bifurcation analysis of the closed-loop system is performed by varying the reference signal according to the selected curve. In this way, regions of state multiplicity and/or instability are identified enabling the choice of the controller parameters values preventing transitions to undesired solution regimes. The method is validated on a simulated problem of start-up control of a jacketed continuous stirred tank reactor. 1. Introduction Large changes of the operating conditions typically occur during the operation of chemical processes. In polymerization reactors, for example, operating conditions are modified to attain product specifications responding to different market demands.1,2 Products obtained during these transients do not meet generally the required specifications and are recycled or sold off at considerably lower market price. Hence, the period needed for the plant to reach the desired steady state must be minimized to reduce the amount of off-specification products. On the other hand, moving nonlinear processes over wide regions of the state and parameter space imposes severe control difficulties. These primarily result from the need to account for system nonlinearities during the design of the closed-loop system ruling out the possibility to directly apply established linear control tools. Furthermore, obstacles may occur due to the presence of process constraints. These impose bounds on state and/or manipulated variables and, if violated, may lead to undesired phenomena as, for example, plant shut-down and thermal runaways.3,4 The design of controllers performing constrained transitions of nonlinear processes has been the object of much research effort. Model predictive control is recognized, in this context, to provide an effective approach to handle input and/or state constraints within an optimal control setting.5 Currently, model predictive control schemes are available stabilizing nonlinear processes around one a priori known steady state (see, for example, (Chen and Allgower)6). However, no stability guarantees are provided when this strategy is implemented to control transitions between different steady states. Furthermore, the implementation of model predictive control laws relies on the * To whom correspondence should be addressed. E-mail: paltimar@ unina.it. † Universita` di Salerno. ‡ CNR. § Universita` del Sannio. | Universita` degli Studi di Napoli “Federico II”. ⊥ Universita` “Federico II”.

online solution of a constrained optimization problem which can be cumbersome or even impossible to solve. Hence, model predictive controllers are typically used, in industrial practice, to control nonlinear processes around a fixed steady state and are introduced only when plant stabilization has been preliminarily achieved via standard feedback regulators as, for example, PID controllers. Alternatively, the problem of controlling process transitions can be handled by adaptive control techniques.7 These are characterized by the use of control schemes involving a feedback controller and an adaptation strategy and have been proved to be effective for controlling nonlinear processes over wide operating regions. Nevertheless, they can be difficult to implement because of the complexity in the controller design. A simpler approach relies on the implementation of gain scheduling8 (GS). In its standard formulation, this technique is based on the use of a family of linear feedback controllers, each of them guaranteeing stability and desired output behavior around a different steady state. Plant steady states and local feedback controllers are, in this context, parametrized by a suitable set of reference variables, typically function of output variables. Hence, the transition between two given steady states is achieved by a step change of the reference signal issued to the closed-loop system. This can be thought of as a switch between the local feedback controllers associated to the initial and the final steady state and may cause constraints violation due to large changes in process dynamics. This can be prevented by adding to the closed-loop system a predictive reference (PR) controller (or reference governor (RG)) which modifies online the reference signal so as to fulfill the constraints over a finite time horizon.9-12 PR-control entails computing a reference input sequence enforcing process constraints based on the prediction of the future evolution of the closed-loop system. This can be done, for example, by constrained minimization of a process cost function over a finite time horizon. However, this may require significant computational resources and does not ensure, in

10.1021/ie9001553 CCC: $40.75  2009 American Chemical Society Published on Web 08/31/2009

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

general, transition to the desired steady state. A more effective approach relies on the parametrization of the reference inputs. In this context, the plant is constrained to move along a prescribed curve of steady states. This restricts the search of a feasible reference input sequence to a one-dimensional subset of the reference variables domain. As a result, a significant reduction of the computational burden is achieved. Moreover, convergence to the desired steady state can be, in this context, proved provided that the steady states of the selected curve fulfill process constraints.9,10 However, it must be emphasized that, with this approach, the evolution of the controlled plant is primarily determined by the choice of a specific curve of steady states. In fact, if steady states close to the border of the process constraint set are selected, undesired reductions of the allowed variations of the manipulated variables are observed. This may in turn lengthen the transient needed for the plant to reach the desired steady state. Despite the practical relevance of these arguments, no studies have been performed of methodologies to adequately compute curves of steady states for PR-control applications. Typically, reference inputs are restricted to the segment connecting the points of the reference variables domain associated to the actual and the desired steady state.9,10 In general, no guarantee about the optimality of the steady states identified by this segment can be provided. These steady states may result, indeed, near or even crossing the boundaries of the process constraint set and/or correspond to unsatisfactory controllability properties of the uncontrolled plant. A further obstacle to the combined application of GS and PR-control is the lack of methodologies for systematically assessing the global stability features of the GS-closed-loop system. In practice, the global performance evaluation of GScontrollers is based on simulation studies.8 However, this approach may fail to detect the presence of multiple stable solution regimes or identify robustness margins of the desired steady states from instability boundaries. In these conditions, sudden variations of the reference signal may cause transition to undesired regimes. In this paper, a novel approach to control transitions of nonlinear processes is presented. The proposed method makes use of GS and PR-control ideas and extends approaches from nonlinear dynamics and flight control which appeared in the literature.13-15 In this context, feedforward control is implemented to force the trajectories of nonlinear systems to track a prescribed curve of steady states and parametric continuation is performed to compute a parametrized family of feedback controller matrices ensuring closed-loop global stability. Inspired from this idea, we here propose to implement PR-control by performing the selection of a feasible reference input sequence so that the plant moves along an optimal curve of steady states. In particular, as in Wang et al. and Richardson et al.,14,15 we use parametric continuation to select the gains of the GScontroller but, at the same time, use optimization techniques to select (or tailor) the curve of steady states to be tracked in order to fulfill a set of control requirements. Thus, an optimal curve is computed so that its steady states lie far from the border of the process constraint set and correspond to satisfactory controllability properties of the uncontrolled plant. Moreover, bifurcation analysis is performed by varying the reference signal according to the computed curve. In this way, regions of state multiplicity and/or instability are identified providing guidelines to select the controller parameters values preventing transitions to undesired regimes.

9129

The paper is structured as follows. In section 2, guidelines to design PR and GS-controllers are summarized and limits of the current approach to the combined application of these two techniques are analyzed. In section 3, a qualitative description of the proposed approach is presented. In section 4, the problem of finding an optimal curve of steady states is formulated. In section 5, the application of bifurcation analysis to the design of the GS-closed-loop system is discussed. In section 6, results of the application of the proposed approach to a simulated problem of start-up control of a jacketed continuous tank reactor (CSTR) model are described. Final remarks end the paper. 2. Control Structure Design In this section, guidelines to design PR and GS-controllers are described. In particular, the main elements of classical GS are summarized in subsection 2.1. Then, the PR-control scheme proposed by Bemporad9 and its application to the class of problems considered in this paper are discussed in subsection 2.2. 2.1. Gain Scheduling. Let the plant to be controlled be described by the following nonlinear system:

{

x˙(t) ) f(x(t), u(t)) y(t) ) g(x(t), u(t))

(1)

where x(t) ∈ Rn is the state vector, y(t) ∈ Rp is the output vector, and u(t) ∈ Rm is the input vector. State and input variables are assumed to be constrained within the following compact subset of the state-parameter space: C ) {(x, u) ∈ Rn × Rm: xmin < x < xmax, umin < u < umax} (2) We define a plant steady state or operating point as a point (xs,us) such that f(xs,us) ) 0 and refer to the set E ) {(xs, us) ∈ Rn+m: f(xs, us) ) 0}

(3)

as the set of steady states of system 1. If f is continuously differentiable and rank[Dxf, Duf] ) n, the set E is described by an m-dimensional manifold of the state-parameter space. In the following, we will refer to this manifold as the bifurcation manifold. To simplify notation, we introduce the single variable p ) (x,u) to denote the n + m state and input variables. We assume that the following parametrization for the set E is available: Ω(σ):D ⊂ Rm f E ⊆ Rn+m

(4)

with σ referred to as the vector of scheduling variables. Typically, the set E is parametrized in terms of plant outputs, that is σ ) σ(y). However, situations might occur where this is not possible. In these cases, differential geometry can be exploited to construct a parametrization of E.16 Once the parametrization (4) is obtained, the typical method of GS is to linearize system 1 about ps ) Ω(σ) leading to the following linear parameter varying model:8

{

·

[ ∂x∂f (x (σ), u (σ))]δx + [∂u∂f (x (σ), u (σ))]δu ∂g ∂g δy ) [ (x (σ), u (σ))]δx + [ (x (σ), u (σ))]δu ∂x ∂u

δx )

s

s

s

s

s

s

s

(5)

s

where δx ) x - xs(σ), δu ) u - us(σ) are the deviations of state and input variables from Ω(σ) ≡ (xs(σ),us(σ)). Then, a linear feedback control strategy can be implemented upon system 5 to guarantee stability and desired output behavior as

9130

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

domain. In particular, the reference input σk provided to the GS-closed-loop-system is chosen as a function of σk-1 and σf(kTp) as follows:9,10 σk ) σk-1 + R[σf(kTs) - σk-1],

R ∈ [0, 1]

(8)

On the basis of eq 8, only reference inputs lying on the segment S(σk-1,σf(kTp))⊂D connecting σk-1 and σf(kTp) are considered as candidates when selecting a feasible reference input σk. We here assume that σf(kTp) ) const:) σf. Under this assumption, eq 8 restricts the search of a feasible reference input sequence {σk}k∈N to the segment S(σ0,σf) connecting σ0 and σf. This segment can be parametrized as follows: Γ(R) ) σ0 + R[σf - σ0], Figure 1. Block scheme of the controlled plant: (a) GS-closed-loop system; (b) GS-closed-loop system and PR-controller.

σ varies. In this paper, we consider the following state feedback GS-control law: (6) where K(σ) denotes a parametrized family of feedback gain matrices. Substituting (6) into (1) gives the governing equations of the GS-closed-loop system:

{

x˙ ) f(x, u(x, σ)) y ) g(x, u(x, σ))

(7)

Since closed-loop dynamics are only affected by the scheduling variables σ, we will use p(t,x0,w) to denote the evolution of (7) when σ ) w and x(0) ) x0. It is worth noting that the control law (6) is composed of two distinct terms: a feedforward uff(σ) and a feedback ufb(x,σ). From this perspective, the closed-loop system (7) can be represented by the block scheme displayed in Figure 1a. Here, σ is the reference fed to the closed-loop system and identifies, based on the filtering action of the feedforward control unit, the desired steady state ps ) Ω(σ). Furthermore, feedback controller gains are scheduled according to the evolution of σ to enforce stability of the desired steady states as plant operating conditions vary. 2.2. Predictive Reference Control. According to the block scheme of Figure 1a, transition between an initial steady state ps0 ) Ω(σ0) and a final steady state psf ) Ω(σf) can be achieved by performing a step change of the reference signal σ between σ0 and σf. However, if large changes in process dynamics occur between ps0 and psf, such a step change may cause process constraints violation. To avoid this, the evolution of the reference signal can be adequately controlled. This can be done by adding to the GS-closed-loop system a discrete time device, called reference governor (RG), which modifies online the reference signal so as to enforce the fulfilment of the constraints.9 Reference control is, in this context, performed in a predictive manner: a feasible sequence of reference inputs is selected based on the prediction of the future evolution of the GS-closed-loop system. Since this requires a finite computational time, the reference governor operates in discrete time modifying the reference signal every period Tp. In particular, during the time interval [kTp, (k+1)Tp] a constant reference input σk is applied guaranteeing that process constraints are not violated over the horizon interval t ∈ [kTp, kTp+Th]. To reduce the computational burden, the search of a feasible reference input sequence is restricted to a one-dimensional subset of the reference variables

(9)

In this setting, a feasible reference input sequence {σk}k∈N ≡{S(Rk)}k∈N can be computed by solving at each instant t ) kTp the following scalar nonlinear program:9 Rk )

u(x, σ) ) us(σ) - K(σ)(x - xs(σ)): ) uff(σ) + ufb(x, σ)

R ∈ [0, 1]

{

max R

R∈[Rk-1,1]

s.t. p(t, x(kTp), Γ(R)) ∈ C,

∀t ∈ [kTp, kTp + Th] (10)

It is important to note that, when applying eq 10, only the restriction of the control law (6) to S(σ0,σf) is implemented. Therefore, closed loop dynamics can be described as follows:

{

x˙(t) ) f(x(t), u(x(t), Γ(R))) y(t) ) g(x(t), u(x(t), Γ(R)))

(11)

Accordingly, the overall controlled plant can be represented by the block scheme reported in Figure 1b. Here, Rf ) Γ-1(σf) is fed to the reference governor which computes online a feasible reference input R. Remark 1. In accordance with Bemporad,9 indispensible prerequisite to achieve constrained transition to psf when applying the PR-control scheme (9-10) upon (11) is fulfilling the following conditions: (a) Ω(Γ(R)) ⊂ C ∀R ∈ [0, 1] (b) if the reference signal R(t) is monotonically convergent to R*, ∀x(0) lim p(t, x(0), Γ(R(t))) ) Ω(Γ(R*)) tf+∞

Condition a is needed to prevent that the border of C is approached in steady state while condition b ensures that no transition to undesired solution regimes takes place when the reference signal R is varied. Unfortunately, situations can likely occur where conditions a and b are not fulfilled. Concerning condition a, it must be noted that steady states ps ∈ Ω(S(σ0,σf)) violating process constraints may be found due to the nonlinearity of Ω(σ). An example showing this possibility is reported in Figure 2. On the other hand, the presence of undesired solution regimes coexisting with ps ∈ Ω(S(σ0,σf)) may prevent the fulfilment of condition b. To prove it, observe that if undesired stable solution regimes coexisting with ps ∈ Ω(S(σ0,σf)) and fulfilling process constraints are found, sudden changes in R might cause the trajectories of the GS-closed-loop system to leave the stability regions of ps ∈ Ω(S(σ0,σf)) without violating process constraints. In this case, application of the PR-control scheme (9-10) might move the GS-closed-loop system to an undesired solution regime even if an infinite horizon Th is used.

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

9131

Figure 2. Projections of the segment S(σ0,σf) onto the bifurcation manifold. Because of the shape of the bifurcation manifold, steady states ps ∈ Ω(S(σ0,σf)) violating process constraints are found.

Remark 2. To prevent transition to undesired solution regimes, Lyapunov functions are employed in Gilbert and Kolmanovsky17 to get an explicit characterization of the stability regions of ps ∈ Ω(S(σ0,σf)). In this framework, only steady states ps ∈ Ω(S(σ0,σf)) whose stability regions include the current state of the GS-closed-loop system are considered as candidates when selecting a feasible reference input σk ∈ S(σ0,σf). However, Lyapunov functions can be derived only for few classes of nonlinear systems restricting the range of applicability of this approach. Therefore, the global stability features of GScontrollers are frequently investigated by means of extensive numerical simulations.8 Obviously, this approach cannot, in general, reveal the presence of undesired solution branches.15 Remark 3. It must be stressed that even if a systematic analysis of the global stability features of the GS-closed-loop system is possible based, for example, on the use of Lyapunov functions, the presence of multiple stable solution regimes may significantly restrict the stability regions of the desired steady states ps ∈ Ω(S(σ0,σf)). In this case, the reference signal must be slowly varied to prevent transitions to undesired solution regimes eventually lengthening the transient needed for the GSclosed-loop system to reach psf. Remark 4. Steady states ps ∈ Ω(S(σ0,σf)) can present characteristics preventing fast transition to psf even if conditions a and b hold. In particular, steady states ps ∈ Ω(S(σ0,σf)) close to the border of C may be found. In these conditions, undesired reductions of the allowed variations of the manipulated variables are observed forcing the reference governor to slowly vary the reference signal to ensure the fulfilment of the constraints. Moreover, the closed-loop system may exhibit undesired output behavior around ps ∈ Ω(S(σ0,σf)). This can happen even when the closed-loop system is carefully designed due to unsatisfactory controllability of the open-loop system about ps ∈ Ω(S(σ0,σf)).18-22 For example, unstable steady states ps ∈ Ω(S(σ0,σf)) may be found requiring large control effort to be stabilized. 3. Problem Formulation and Outline of the Solution Procedure In accordance with remarks 1-4, the following limits of application of the PR-control scheme (9-10) can be identified: (1)The steady states ps ∈ Ω(S(σ0,σf)) may be near or crossing the border of C and the uncontrolled plant may exhibit

Figure 3. Projections of curves of the reference variables domain onto the bifurcation manifold.

unsatisfactory controllability properties about ps ∈ Ω(S(σ0,σf)). (2)Undesired stable solution regimes coexisting with the steady states ps ∈ Ω(S(σ0,σf)) may be found. The first limit arises as a result of restricting a priori the reference inputs to S(σ0,σf) without taking into account the characteristics of the corresponding steady states ps ∈ Ω(S(σ0,σf)). To overcome this limit, we extend the application of the PR-control scheme (9-10) to curves of the reference variables domain D connecting σ0 and σf. Infinite curves Γ ⊂ D connecting σ0 and σf are, indeed, available when dim(D) ) m > 1 allowing to optimize the characteristics of the steady states ps ∈ Ω(Γ) around which the plant is operated during the transition (see, for example, Figure 3). In particular, the curve Γ ⊂ D should be selected so that no steady states ps ∈ Ω(Γ) close to the border of C are found and the uncontrolled plant exhibits satisfactory controllability properties around ps ∈ Ω(Γ). Here, we consider the set of continuously differentiable curves Z(σ0,σf) ⊂ D connecting σ0 and σf. This means that it is always possible to construct a continuously differentiable and invertible parametrization mapping the interval [0,1] onto a curve Γ ∈ Z(σ0,σf): Γ(R):[0, 1] f D ⊆ Rm Γ(0) ) σ0, Γ(1) ) σf

(12)

Once a parametrization of the type (12) is selected, a feasible reference input sequence {wk}k∈N ) {Γ(Rk)}k∈N can be computed by solving the nonlinear program (10) at each period t ) kTp.

9132

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

The second limit imposes to account for the global stability features of the GS-closed-loop system. In particular, both the selection of the curve Γ and the design of the GS-closed-loop system should be performed so that no undesired solution regimes coexisting with the steady states ps ∈ Ω(Γ) are detected as R varies. It is, indeed, clear that depending on how the GSclosed-loop system is designed, it might be impossible to find a curve of globally stable steady states Ω(Γ). This happens, for example, when the GS-closed-loop system is designed so that undesired solution regimes coexisting with ps0 and/or psf are detected. Therefore, Γ(R) and K(Γ(R)) should be computed so that the following requirements are fulfilled as R varies: • The steady states Ω(Γ(R)) lie as far as possible from the border of C. • No undesired solution regime coexisting with Ω(Γ(R)) is detected. • The closed-loop system exhibits desired output behavior about Ω(Γ(R)). It must be noted that the requirement of satisfactory controllability properties of the uncontrolled plant is here implicitly imposed by demanding that the closed-loop system does exhibit desired output behavior about Ω(Γ(R)). It is also worth observing that computing Γ(R) and Κ(Γ(R)) so as to achieve the addressed control objectives results in tailoring the bifurcation diagram (the nature, the number, and the stability characteristics) of the GS-closed-loop system as R varies. To solve the previous problem, information about closedloop performance and distances of plant steady states from the border of the process constraint set and from regions of the parameter space characterized by state multiplicity must be derived. In this framework, we observe that several methods are available in literature to estimate closed-loop performance and can be used to structure the formulated problem.22 Also, for the class of process constraints considered in this paper, functions providing information about the distance of plant steady states from the border of the process constraint set can be easily derived. On the contrary, difficulties might arise when attempting to obtain information about the distance of plant steady states from regions of the parameter space characterized by state multiplicity. In principle, nonlinear constraints can be formulated imposing that the distances of ps ∈ Ω(Γ) from the boundaries of such regions do not exceed a prescribed value.23,24 This approach has been applied for the optimization based design of nonlinear processes.25,26 Nevertheless, its application to compute a parametrized family of steady states, such as Ω(Γ), may become cumbersome. In this case, a set of nonlinear constraints should be imposed for each ps ∈ Ω(Γ), significantly increasing the complexity of the problem. To overcome this difficulty, we separately compute Γ and K(Γ). In particular, the following two-step solution procedure is proposed: (1) Optimal reference trajectory shaping. A curve Γ is computed so that Ω(Γ(R)) lies as far as possible from the border of C and the open-loop system linearization about Ω(Γ(R)) exhibits satisfactory controllability properties as R varies. (2) Robust GS-closed loop design. K(Γ(R)) is computed so that the GS-closed-loop system does not exhibit undesired solution regimes and desired output behavior is observed about Ω(Γ(R)) as R varies. In accordance with step 1, the curve Γ is computed before designing the closed-loop system by only exploiting information about the uncontrolled plant. The selection of Γ is performed so to enforce satisfactory controllability of the open-loop system linearization about ps ∈ Ω(Γ). Fulfilling this requirement is,

indeed, necessary to ensure the existence of a family of feedback gain matrices K(Γ(R)) ensuring desired output behavior around ps ∈ Ω(Γ) as R varies. This approach to select Γ sidesteps the need to compute the distances of the steady states Ω(Γ) from bifurcation boundaries. Conversely, step 2 aims at computing a family of feedback gain matrices K(Γ(R)) ensuring satisfactory output behavior and ruling out the presence of undesired solution regimes. In sections 5-6, a simple method to compute such a K(Γ(R)) is presented ensuring feasibility of the proposed solution procedure. 4. Optimal Reference Trajectory Shaping In this section, the problem of computing an optimal reference trajectory is described. A rigorous formulation of this problem is presented in subsection 4.1 while computational issues concerning its numerical solution are discussed in subsection 4.2. 4.1. Optimization Problem. To solve the optimal reference trajectory shaping problem, functions providing information about the distance of plant steady states from the border of C and open-loop controllability must be defined. Integrating these functions over the curve Ω(Γ) as Γ varies gives a functional defined in the set Z(σ0,σf). Thus, an optimal Γ can be obtained by minimization of such a functional. Therefore, given the steady states ps0 and psf ∈ E ∩ C, a vector function L(p):E ⊂ Rn+m f Rl whose elements provide measures for open-loop controllability about p, a function M(p):E ⊂ Rn+m f R quantifying the distance of p from the border of C, the optimal reference trajectory shaping problem can be formulated as the problem of computing a parametrization Γ(R):[0,1] f D ⊂ Rm verifying that Γ(R) )

min Φ(Γ) )

Γ∈Z(σ0,σf)

min

Γ∈Z(σ0,σf)

subject to



1

0

l

(Ψ′((

∑w

1,iLi(Ψ))

+ w2M(Ψ))) dR

(13)

i

Ψ(R) ∈ C L(Ψ(R)) e η

∀R ∈ [0, 1] ∀R ∈ [0, 1]

(14) (15)

where Ψ(R)) Ω(Γ(R)), Ψ′(R))DRΨ(R), w1 and w2 are weighting constants and η ∈ Rl is a vector of bounds. It is worth noting that, with constant M and Li, the functional Φ defined in eq 18 increases as the length of the curve spanned by Ψ(R) increases. This is in line with the objective of reducing the transient needed for the plant to reach the final steady state. To prove it, consider the problem of performing a transition between two steady states lying close to each other. In this case, the use of a short curve restricted to the bifurcation manifold and connecting the two considered steady states would be probably sufficient to guarantee fast transition while fulfilling process constraints. Nevertheless, if the steady states ps ∈ Ω(Γ) are selected only based on controllability measures, a longer curve may be found moving far away from the region surrounding the initial and the final steady state. As a result, the transient needed for the plant to reach the final steady state may be lengthened even though the open-loop system linearization about the selected steady states ps ∈ Ω(Γ) exhibits satisfactory controllability properties. For the function L(p), the inverse of the minimum singular value of the process steady state gain matrix and the sum of the real parts of the right-half plane eigenvalues of the openloop system linearization are used. A small value for the inverse

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

of the minimum singular value indicates that the steady state gain matrix is not close to being singular and that the plant easily moves between different steady states27,28 while small values for the real parts of the right-half plane eigenvalues ensure that no large control efforts are needed to achieve stability.29 Therefore, L(p) is defined as follows:

[∑

1/VS,min(p)

L(p) )

Re{RHPEi}

i

]

(16)

where VS,min(p) and RHPEi denote the minimum singular value of the process steady state gain matrix and the right-half plane eigenvalues of the open-loop system linearization at p, respectively. It must be remarked that more effective and sophisticated methods to estimate open-loop plant controllability are currently available and could be used to structure the problem (13-15). We refer the reader to Sakizlis et al.22 for a survey of these approaches. Finally, the function M(p) is defined as the sum of the distances of state and manipulated variables from the boundaries of the process constraint set C. This gives the following expression: n

M(p) )

∑ [|x - x

⊥ -1 max,i | Q

⊥ + |x - xmin,i | Q-1] +

i)1

∑ [|u -

u⊥max,i | Q-1

+ |u -

⊥ umin,i | Q-1]

(17)

i)1

⊥ ⊥ and uk,i denote the projections of x and u onto the where xk,i planes x ) xk,i and u ) uk,i of the state and parameter spaces, respectively. The function M(p), also known as inverse barrier function,30 takes positive values in the interior of the process constraint set C and becomes infinitely large as p approaches the border of C. Hence, replacing (17) in (13) has the effect of penalizing the selection of curves Γ corresponding to steady states ps ∈ Ω(Γ) close to the border of C. This approach resembles the barrier method implemented to handle nonlinear constraints when minimizing a given objective function.30 However, it is important to note that the spirit by which it is here applied is completely different. The function M(p) is, indeed, introduced in (13) to prevent that steady states ps ∈ Ω(Γ) close to the border of C are found while the fulfilment of process constraints is separately enforced by fulfilling (14). 4.2. Computational Issues. Solving problem (13-15) requires the finding of a parametrized curve of the reference variables domain minimizing the functional Φ. This gives an infinite dimensional optimization problem. Therefore, discretization is necessary to make the problem (13-15) numerically affordable. To this aim, the parametrization Γ(R) is replaced by the sequence ΓΝ ){σ1, ..., σN} and the optimal reference trajectory shaping problem is recast in the following form:

ΓN ) min ΦN(ΓN) ) ΓN

min ΓN

N

l

k)2

i

∑ (( ∑ w

1,iLi(pk))

+ w2M(pk))|Ω(σk) - Ω(σk-1)| Q (18)

Ω(σk) ∈ C L(Ω(σk)) e η

Once the problem (18-20) is solved, the elements of the optimal sequence ΓΝ ){σ1, ..., σN} can be interpolated to get an approximate solution Γ*(R) to (13-15). In this framework, the objective is to compute a parametrization Γ*(R) guaranteeing that the corresponding steady states Ω(Γ*(R)) exhibit, as R varies, characteristics similar to those observed for the steady states ΨN ) {Ω(σ1), ..., Ω(σN)}. The choice of the steady states in ΨN is, indeed, optimized by solving (18-20), whereas there is no control on the steady states Ω(Γ*(R)) with Γ*(R)∉{σ1, ..., σN}. However, solving (18-20) may lead to sequences ΓΝ ) {σ1, ..., σN} corresponding to steady states ΨN ) {Ω(σ1), ..., Ω(σN)} far from each other. In this case, the characteristics of Ω(Γ*(R)) may likely result, as R varies, significantly different compared to those observed for ΨN ) {Ω(σ1), ..., Ω(σN)}. In particular, steady states Ω(Γ*(R)) with Γ*(R)∉{σ1, ..., σN} may be found violating process constraints or corresponding to unsatisfactory controllability properties of the open-loop system linearization. This can be prevented by interpolating ΓΝ at each step while solving (18-20) and checking that Ω(Γ*(R)) fulfill the feasibility conditions (14-15) as R varies. Obviously, this significantly increases the computational resources needed to compute Γ*(R). Alternatively, the following constraint can be added to (18-20): -δ e Ω(σk) - Ω(σk-1) e δ,

m

∀k ∈ {1, ..., N} ∀k ∈ {1, ..., N}

(19) (20)

where the norm is defined by recourse to a weighting matrix Q. In particular, |p|Q ) pTQp, with Q positive definite.

9133

∀k ∈ {1, ..., N}

(21)

where δ ∈ Rn+m is a vector of bounds. This forces the steady states ΨN ){Ω(σ1), ..., Ω(σN)} to be close to each other enabling us to ensure feasibility of the interpolating branch Ψ(Γ*(R)) by adequately selecting δ. In particular, a feasible solution Γ*(R) can be obtained by solving (18-21) for decreasing values of the elements of δ until the obtained branch Ψ(Γ*(R)) is found to fulfill (14-15) as R varies. It must be here observed that reducing the elements of δ may require an increase in N in order to get a feasible solution to (18-21). Therefore, the following algorithm is followed to compute Γ*(R): 1. select initial guesses for δ and N 2. solve (18-21) 3. if a solution to (18-21) is found, interpolate the computed sequence ΓΝ to get Γ*(R) otherwise increase N and repeat step 2 4. if steady states Ω(Γ*(R)) not fulfilling (14-15) are found, reduce the elements of δ and repeat the previous steps until a feasible Ω(Γ*(R)) is found We note that an initial guess for δ can be obtained based on the physics of the process, namely, looking at the characteristic variations of state and manipulated variables. Then, an estimate of the minimum N ensuring feasibility of (18-21) can be computed as the minimum number of steady states ps,i ∈ Ω(S(σ0,σf)) fulfilling -δ < ps,i - ps,i-1 < δ and ps,N ) psf. An initial guess for N can be, therefore, derived by slightly increasing this value. The proposed algorithm enables one to get a feasible solution Γ*(R). Nevertheless, it must be remarked that the sequence ΓN solving (18-21) and, hence, its interpolation Γ*(R) approximate with increasing accuracy the curve Γ solving (13-15) when N is increased. Therefore, a rigorous way to get an interpolating function close to the solution of (13-15) is, after obtaining a feasible Γ*(R), to solve (18-21) for increasing N values until no significant changes in the observed solution occur. The required computational resources needed to perform this task can be reduced by employing the sequence ΓN solving (18-21)

9134

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

with N ) K to construct an initial sequence when solving (18-21) with N ) K + 1. 5. Bifurcation Analysis of the Gain-Scheduled Closed-Loop System Bifurcation analysis has been extensively applied to characterize the behavior of nonlinear systems and can be considered an established tool in the analysis of engineering processes. The studies conducted in this area have been traditionally aimed at improving the understanding of the qualitative dynamics of the process units. A representative example of such an approach is the seminal paper on the dynamics of an exothermic continuous tank reactor by Ray and his co-workers.31 Over the past few decades, the use of bifurcation analysis has been proposed to address synthesis problems in a rigorous manner. In this framework, two different approaches to the application of bifurcation analysis to control system design can be distinguished. The first approach assumes that a control structure has been preliminary selected and makes use of parametric continuation to identify regions of the controller parameter space characterized by unstable plant operation and/or undesired global dynamics as, for example, multistability.13,32-34 The second approach, also known as bifurcation control,35 incorporates information about the bifurcation characteristics of the uncontrolled plant into the synthesis of the control structure. In this context, control strategies capable of tailoring the bifurcation characteristics of nonlinear processes are formulated.14,36,37 In what follows, we discuss the application of bifurcation analysis and parametric continuation to compute a parametrized family of feedback gain matrices K(Γ(R)). The presented methodology enables one to account for the global stability features of the GS-closed-loop system and extends the approach proposed in Richardson et al.15 by providing a framework for the design of the control requirements that improve the global stability features of the GS-closed-loop system. To compute K(Γ(R)), feedback gain matrices enforcing a set of control requirements (e.g., settling time, risetime, etc.) about specified steady states of the curve Ω(Γ) are interpolated.38,39 With this approach, poor controller performance or even instability can occur due to hidden coupling terms.40,41 Moreover, control requirements at specified steady states are generally defined based on trial-and-error procedures aimed at finding a good compromise between closed loop performance and reliability of the control law. However, the effect of these choices on the global stability features of the controlled plant is generally ignored. To overcome these limitations, bifurcation analysis can be performed to investigate the effect of interpolating feedback gains and/or modifying controller parameters on the occurrence of multiple stable solution regimes of the GS-closed-loop system. To account for the dependence of closed-loop dynamics on controller parameters, we recast eqs 11 in the following form:

{

x˙(t) ) f(x, u(x(t), Γ(R), λ)) y(t) ) g(x, u(x(t), Γ(R), λ))

(22)

with u(x(t),Γ(R), λ) described as follows: u(x, Γ(R), λ) ) us(Γ(R)) + K(Γ(R), λ)(x - xs(Γ(R))) ) uff(R) + ufb(x, R, λ)

(23)

Here, λ denotes the vector of parameters which is necessary to fix to compute the feedback controller gains. For example, λ defines the desired poles of the closed-loop system linearization when applying a pole placement procedure or the elements of

Table 1. Reactor Parameter Values V (L) F (g/L) Cp (J/g · K) Cpc (J/g · K) U (J/min · K) E1/R (K) E2/R (K) k1 (min-1)

100 1000 0.239 0.25 18462 8355 3759 7.46 × 1011

k2 (min-1) ∆Η1 (J/mol) ∆Η2 (J/mol) Tf (K) Tcf (K) CA,in (mol/L) CB,in (mol/L)

16632 20596 3913 300 293 1 0

weighting matrices when making use of a linear quadratic regulation algorithm.29 Therefore, the effect of varying λ on the global stability characteristics of the steady states ps ∈ Ω(Γ) must be assessed. To this aim, the evolution of the bifurcations of the GS-closed-loop system (22) in the R-λ space is analyzed. In so doing, possible bifurcations responsible for the instability of the steady states ps ∈ Ω(Γ) and/or the occurrence of multistability are detected. Hence, λ values can be selected to guarantee that the steady states ps ∈ Ω(Γ) exhibit prescribed margins from the bifurcation boundaries. Note that the range of λ over which the analysis should be performed must be selected according to the specific control problem to guarantee satisfactory closed-loop behavior about the steady states ps ∈ Ω(Γ). 6. Simulation of the Start-up of a Continuous Stirred Jacketed Tank Reactor To prove the effectiveness of the proposed approach, we present in this section the results of its application to the problem of controlling the start-up of a CSTR where the series of two first order irreversible exothermic reactions A f B f C takes place. Depending on the reaction heat values and the activation energies, these reactors can exhibit multiple stable solution regimes and nonminimum-phase behavior.42 6.1. Closed-Loop System Design. Under standard modeling assumptions, reactor dynamics can be described by the following nonlinear system:

( )

ER1 dCA Q ) (CA,in - CA) - k1 exp C dt V RT A

( )

(24)

ER2 dCB Q ) (CB,in - CB) - k2 exp C + dt V RT B

( )

k1 exp -

ER1 C RT A

( ) ( )

∆H1k1 ER1 dT Q ) (Tin - T) exp C dt V FCp RT A ER2 ∆H2k2 U exp C (T - Tc) FCp RT B VFCp Qc dTc U ) (Tc,in - Tc) + (T - Tc) dt Vc VFCpc

(25)

(26)

(27)

where CA is the concentration of the species A, CB is the concentration of the species B, T is the reactor temperature, Tc is the coolant temperature, Q is the feed flow rate, and Qc is the coolant flow rate. We assume, in the following, that Q and Qc are the inputs to the system which can be modified. The other parameters appearing in (24-27) have the usual meaning and are kept at constant values reported in Table 1. The problem we address is to perform transition of (24)-(27) from ps0 to psf as described in Table- 2. ps0 corresponds to reaction extinction conditions while psf is characterized by high conversion of the intermediate product B. In particular, psf is

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009 Table 2. Process Constraints, Initial and Final Equilibrium

CA (moli/L) CB (moli/L) T (K) Tc (K) Q (L/min) Qc (L/min)

ps0

psf

lower bound

upper bound

0.89 0.11 305 300.1 70 51.19

0.19 0.57 354 329.1 100 51.2

1 1 400 400 110 70

0 0 293 293 40 30

( ) ( )

( )

xs3(σ) ) σ2

(36)

E1 σ h (σ) Rσ2 1 2 + (σ1 - x1,in) E2 h (σ)σ1∆H2(V/U) + σ2 (37) k2 exp Rσ2 1

xs6(σ) )

σ1UV(A - B) CD + EF

( )

A ) k1 exp

(28)

(39)

E2 h (σ) Rσ2 2

( )

B ) k2(σ1 - x1,in)∆H2 exp

( ) ( ) ( )

D ) exp

E1 h (σ) Rσ2 1

E ) FCpc exp

(30)

(32)

E1 h (σ) Rσ2 1

C ) k2∆H2FCpcVσ1(σ1 - x1,in)

(29)

(31)

E2 Rσ2

E1 (σ - x4,in)U - k1σ1h2(σ)V Rσ2 2

us1(σ) ) xs5(σ)

(40)

us2(σ) ) xs6(σ)

(41)

with h1(σ) and h2(σ) defined as follows: k1(σ1 - x1,in - x2,in) exp h1(σ) )

(34)

( )

k1σ1 exp

(33)

Equations 32-33 describe the dynamics of the filter. Here, τ1 and τ2 are time constants while u1 and u2 are the inputs to the filter and define the desired values of feed and coolant flow rates. In the following, we assume that τ1 ) τ2 ) 0.1. Therefore, we pick u1 and u2 as manipulated variables and apply the methodology described in the previous sections to implement the transition between ps0 and psf. The concentration of the species A and the reactor temperature are chosen as scheduling variables, that is σ ) (x1, x3). Therefore, eqs 28-33 are solved under steady state conditions with unknown variables (x2, x4, x5, x6, u1, u2). This gives the following map Ω(σ) ) (xs1(σ), xs2(σ), xs3(σ), xs4(σ), xs5(σ), xs6(σ), us1(σ),us2(σ)): xs1(σ) ) σ1

(38)

where

x6 dx4 U ) (x4,in - x4c) + (x - x4) ) g4(x, u) dt Vc VFCpc 3

dx6 1 ) - (x6 - u2) ) g6(x, u) dt τ2

E1 σ Rσ2 1 (σ1 - x1,in)

xs5(σ) )

F ) (σ1 - x1,in) exp

dx5 1 ) - (x5 - u1) ) g5(x, u) dt τ1

( )

-k1Vexp -

( )

x5 ∆H1k1 E1 dx3 ) (x3,in - x3) exp x dt V FCp Rx3 1 E2 ∆H2k2 U exp x (x - x4) ) g3(x, u) FCp Rx3 2 VFCp 3

(35)

( ) ( )

xs4(σ) )

x5 E1 dx1 ) (x1,in - x1) - k1 exp x ) g1(x, u) dt V Rx3 1

x5 E2 dx2 ) (x2,in - x2) - k2 exp x + dt V Rx3 2 E1 x ) g2(x, u) k1 exp Rx3 1

xs2(σ) ) -σ1h(σ)

-k1(V/U) exp -

unstable and coexists with a stable and an unstable steady state. In this framework, we assume that state and input variables are constrained within the set C described in Table 2. To avoid sudden variations of the control law, we assume the control inputs to the plant to be generated through a lowpass filter. Therefore, we pose x ≡ (CA, CB, T, Tc ,Q, Qc) and describe the filter and plant dynamics by means of the following nonlinear system:

( )

9135

( ) ( ) E2 Rσ2

E2 E1 - k2(σ1 - x1,in) exp Rσ2 Rσ2

h2(σ) ) (σ1 - x1,in)∆H1 + FCp(σ2 - x3,in)

(42) (43)

An optimal curve Γ of the scheduling variables domain can be now constructed following the approach illustrated in section 4. To this aim, the nonlinear program (18-21) is solved with parameters values reported in Table 3 and the spline toolbox of Matlab is used to interpolate the resulting sequence of steady states. The interpolating curve Γ and its projection Ω(Γ) onto the bifurcation manifold are described in Figure 4. It is important to note that the evolutions of the manipulated variables u1(R) and u2(R) shown in Figure 4e,f provide the Table 3. Parameters Used to Compute the Optimal Γ w1

w2

N

η

δ

(5, 6)

10

20

(2, 3)

(0.1, 0.1, 5, 5, 5, 5, 5)

9136

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

Figure 4. Parameterization Ω(Γ(R)) with Γ solving the optimal reference trajectory shaping problem.

feedforward control law uff(R) which needs to be implemented to achieve transition to psf. Besides this control law, feedback must be implemented to guarantee stability and desired output behavior about Ω(Γ(R)) as R varies. To this aim, a parametrized family of feedback gain matrices K(Γ(R),λ) is computed by applying the linear quadratic regulation (LQR) algorithm (Kailath, 1996). This gives the following GS-control law: u(x, R) ) us(x, Γ(R)) - K(Γ(R), λ)(x - xs(Γ(R))) 1 ) uff(x, Γ(R)) - BT(Ω(Γ(R)))P(Ω(Γ(R)), λ) × 2 (x - xs(Γ(R))) (44) where B ) Du[g(xs(Γ(R)),us(Γ(R))]| and P is the solution of the algebraic Riccati equation: ATP + PA - PBBTP + R(λ) ) 0

(45)

Here, A ) Dx[g(xs(Γ(R)),us(Γ(R))] and R(λ) is a positive definite matrix whose elements are parametrized by λ. In particular, it is assumed that λ is a scalar and R is the diagonal matrix defined by the vector V ) (5λ., 5λ., 5λ., 5λ., 1, 1), that is R ) diag(V). At this point, the PR-control scheme (10) can be applied to perform transition of the GS-closed-loop system to psf. With this objective, the control law (44) is implemented with λ ) 1, which is found to ensure satisfactory output behavior about Ω(Γ(R)) as R varies. Furthermore, an horizon interval Th ) 5 min and a period Τp ) 0.2 min are selected and a bisection algorithm is used to solve the nonlinear program (10). Results of the implementation of (10) are described in Figure 5. Here, the reactor temperature simulation response (Figure 5a) and the evolution of the reference signal imposed by the PR-controller (Figure 5b) are described. The reactor temperature value xs3(1) corresponding to psf is also shown in Figure 5a. It is apparent from Figure 5 that the reference control scheme (10) with the computed Γ fails to bring the GS-closed-loop system to psf. In particular, after an initial rise, the reactor temperature starts to decrease and an undesired stationary value is eventually reached.

Figure 5. Evolution of the GS-closed-loop system under PR-control with the optimal Γ and λ ) 1: (a) reactor temperature simulation response; (b) reference signal.

It is remarkable to note that the observed behavior occurs irrespective of the horizon interval Th used to implement (10). To prove it, it must be observed that the inversion in the pattern initially followed by the reactor temperature occurs at t = 8.2 min when the reference signal is switched between R ) 0.58 and R ) 1. This switch is performed, in accordance to (10), as it enables reaching the desired value of the reference signal R ) 1 and does not produce process constraint violation over the horizon interval t ∈ [8.2, 8.2 + Th]. However, process constraints are found to be fulfilled not only in the range t ∈ [8.2, 8.2 + Th] but ∀ t > 8.2 min. For this reason, the observed switch of the reference signal is performed, in accordance with (10), even if an infinite Th is used. In these conditions, undesired transitions can be prevented by modifying the design of the GS-closedloop system so that undesired solution branches coexisting with Ω(Γ) are removed. With this objective, the approach described in section 5 is applied. Figure 6a shows steady state solution regimes of the GSclosed-loop system as R varies. The reactor temperature is

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

9137

Figure 6. Bifurcation structure of the GS-closed-loop system: (a) solution diagram with R as bifurcation parameter (stable and unstable branches are denoted by solid and dashed line respectively); (b) projections of the saddle node bifurcations in the R-λ plane.

Figure 7. Evolution of the GS-closed-loop system under PR-control with the optimal Γ and λ ) 3: (a) reactor temperature simulation response; (b) reference signal.

chosen as a variable representative of the state of the system. It can be observed that the steady states Ω(Γ(R)) are stable over the entire range R ∈ [0,1]. Nevertheless, undesired solution branches are detected at large R values because of the presence of the saddle node bifurcation point S1. In particular, a stable low conversion and an unstable steady state coexisting with Ω(Γ(R)) are detected at R values larger than RS1. To examine the effects of modifying controller parameters on such coexistence, the evolution of the saddle node bifurcation point S1 is analyzed by varying λ in the range [1, 4] Projections of the saddle node bifurcation points onto the plane R-λ are shown in Figure 6b. It can be observed that no saddle node bifurcation is detected in the range R ∈ [0,1] when λ > 2.3. For these λ values, no undesired steady state coexisting with Ω(Γ) can be found. Therefore, the reference control scheme (10) is applied by implementing the control law (44) with λ ) 3, which is still found to ensure satisfactory output behavior about Ω(Γ(R)) as R varies. The reactor temperature simulation response observed in this case is described in Figure 7. It is apparent that the controlled plant moves to the desired steady state psf. 6.2. Performance Analysis. To analyze the performance of the proposed approach, results of the implementation of the reference control scheme (10) with the optimal Γ computed in the previous subsection and the segment Γ ) S(σ0, σf) are compared. In this framework, we preliminarily describe in Figure 8 the curve Ω(Γ) and Ω(S(σ0, σf)).

It is apparent that the profiles of plant outputs (Figure 8a-d) for the considered curves stay close to each other as R varies. In particular, temperature profiles are barely distinguishable. Moreover, output values far from the bounds of the constraint set are found in both cases. On the contrary, significant differences between the profiles of plant inputs found for Ω(Γ) and Ω(S(σ0,σf)) are observed (Figure 8e-f). Feed flow rate values us1 (Figure 8e) for Ω(S(σ0,σf)) are invariably larger than those found for Ω(Γ), while the contrary holds for the coolant flow rate us2 (Figure 8f). However, while feed flow rate values us1 far from the bounds of the constraint set are observed for both Ω(Γ) and Ω(S(σ0,σf), the us2 profile is, for Ω(S(σ0,σf)), close to its lower bound around R ) 0.2. This means that when R is close to such value, reduced variations of the coolant flow rate are allowed if S(σ0,σf) is used in (10). The effects of this on transition control are shown in Figure 9. Here, the reactor temperature simulation response observed when the reference control scheme (10) is implemented with Γ ) S(σ0, σf) (Figure 9a) and the evolution of the reference signal imposed by the reference controller are reported (Figure 9b). The same value λ ) 3 selected when implementing (10) with the Γ computed in the previous subsection is used in (44). This value guarantees uniqueness of the steady states ps ∈ Ω(S(σ0,σf)). With this setting, differences between the behaviors observed when S(σ0,σf) and Γ are used in (10) must be imputed to different characteristics of the corresponding steady states ps ∈ Ω(S(σ0,σf)) and ps ∈ Ω(Γ). Figure 9 shows that the transition period found when Γ ) S(σ0, σf) is used in (10) is about 30 min. This value is significantly larger than that found with the optimal Γ computed in the previous subsection (=13 min) (Figure 7). The observed difference is mainly due to the presence of steady states ps ∈ Ω(S(σ0,σf)) close to the border of the process constraint set. As shown in Figure 8, the coolant flow rate for ps ∈ Ω(S(σ0,σf)) approaches its lower bound when R = 0.2. In this range, the reference governor is forced to slowly vary R in order to guarantee the fulfilment of the constraints. This clearly appears in Figure 9b. It can be here observed that the reference controller takes a long time to move R away from the range around R ) 0.2. On the contrary, a fast growth in R is found once this region is crossed. 7. Summary and Future Directions A novel method for transition control of nonlinear processes was presented. In accordance with the proposed approach, gainscheduling is implemented to ensure desired output behavior over the operating region of interest while transition between the initial and the final steady state is achieved by predictive

9138

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

Figure 8. Comparison between the parametrizations Ω(Γ(R)) obtained with Γ solving the optimal reference trajectory shaping problem (dashed line) and Γ ) S(σ0,σf) (solid line).

Figure 9. Evolution of the GS-closed-loop system under PR-control with Γ ) S(σ0,σf) and λ ) 5: (a) reactor temperature simulation response; (b) reference signal.

reference control. In so doing, a feasible sequence of reference inputs is computed so that the closed-loop system moves along an optimal curve of the bifurcation manifold. This restricts the search of a feasible reference input sequence to a onedimensional subset of the reference variables domain, reducing the required computational resources. Moreover, nonlinearities responsible for transitions to undesired solution regimes can be, with this approach, promptly identified. It is, indeed, sufficient to trace the solution diagram of the closed-loop system by varying the reference signal according to the selected curve to detect multiple stable solution regimes and/or bifurcations causing instability of the selected steady states. In this framework, the problem of tailoring a curve of the bifurcation manifold and a family of feedback gain matrices ensuring desired output behavior and preventing the occurrence of undesired solution branches was addressed. The following two-steps procedure was proposed to tackle this problem: first, a curve of steady states is computed so that its points lie far

from the border of the process constraint set and correspond to satisfactory controllability properties of the uncontrolled plant (optimal reference trajectory shaping); then, bifurcation analysis of the closed-loop system is performed enabling the choice of the controller parameters values preventing transitions to undesired solution regimes (robust GS-closed-loop design). With this approach, an optimal curve of steady states is computed before designing the closed-loop system by only exploiting information about the uncontrolled plant. At this stage, steady states corresponding to satisfactory controllability of the openloop system linearization are selected. Fulfilling this requirement is, indeed, needed to ensure the existence of feedback controller gains enforcing desired closed-loop behavior. On the other hand, bifurcation analysis is performed enabling a thorough examination of the effect of varying feedback controller parameters on the global stability properties of the overall controlled plant. Therefore, feedback controller parameters can be computed ensuring satisfactory output behavior and ruling out the presence of undesired solution branches. The proposed method was validated on the problem of controlling the start-up of a jacketed continuous tank reactor. In this context, the approach proved to effectively handle the presence of multiple stable solution regimes enabling the removal of a low conversion solution branch responsible for plant-shut down. Also, the importance of adequately shaping a curve of the bifurcation manifold to guarantee fast transition to the final steady state was demonstrated by comparison with the classical reference governor approach.9 In particular, PR-control was implemented by restricting the reference inputs fed to the closed-loop system to (a) the segment connecting the reference inputs associated to the initial and the final steady state and (b) to the curve obtained by solving the optimal reference trajectory shaping problem. In the first case, the steady states spanned by the considered segment were found close to the border of the process constraint set resulting in a large transition period. On the contrary, fast transition was achieved when restricting the reference inputs to the computed optimal curve.

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

It is important to remark that the advantage of the proposed method, compared to alternative existing control approaches as, for example, model predictive control, mainly resides in its capability of fulfilling optimal control requirements and, at the same time, accounting for the presence of critical operability boundaries during the design of the closed-loop system. This approach is expected, in general, to ensure fast transition to the final steady state while preventing the occurrence of undesired phenomena as, for example, secondary reactions or thermal runaways. It is worth remarking that the proposed method was here applied to a problem where the steady states of the uncontrolled plant can be parametrized in terms of internal plant variables. However, situations can be found where such a parametrization cannot be used due, for example, to the occurrence of output and/or input multiplicity. It must be stressed that, in this case, the application of the proposed methodology results unchanged. Obviously, an alternative parametrization must be initially constructed. With this objective, an effective approach is to introduce, based on the idea of pseudo-arc-length,43 a curvilinear coordinates system onto the bifurcation manifold of the uncontrolled plant.16 Here, each steady state is identified by a set of curvilinear coordinates on the bifurcation manifold. This approach works disregarding the occurrence of input and/or output multiplicity providing the basis for a broad application of the proposed methodology. Nevertheless, the application of this parametrization scheme to solve the optimal reference trajectory shaping problem is not trivial and deserves to be thoroughly discussed. Therefore, we leave a detailed analysis of such alternative approach to the parametrization of plant steady states as the subject of a future work. It must be finally stressed that significant improvements of the proposed method are expected to be achieved by simultaneously computing the curve of steady states to be tracked and the feedback controller gains. In particular, not only the feedback controller gains but also the curve of steady states to be tracked should be selected to enforce desired global stability features of the closed-loop system. To this aim, nonlinear constraints can be formulated imposing that the distance of each steady state of the selected curve from bifurcation and performance boundaries do not exceed a prescribed value.26 This approach also enables the handling of process constraints which are implicitly defined by nonlinear equations in the state and input variables. However, its application to the formulated problem becomes cumbersome when a large number of steady states is needed to describe a curve connecting the initial and the final steady state. In this respect, we are currently investigating reduction techniques which allow the decomposition of the optimal reference trajectory shaping problem into a sequence of subproblems whose solution involves few steady states. Acknowledgment The authors kindly acknowledge professor Costin Sorin Bildea for the valuable suggestions provided during the preparation of the manuscript. Notation Acronyms GS ) gain-scheduled PR ) predictive reference RG ) reference governor RHPE ) right-half plane eigenvalue

9139

Notation C ) constraint set D ) reference variables domain Dx ) derivative with respect to x diag ) diagonal matrix defined by vector v E ) set of steady states of the uncontrolled plant L ) controllability function M ) barrier function p ) state-input vector r ) reference input vector S(a, b) ) segment of the reference variables domain connecting a and b Th ) time horizon Tp ) sampling time u ) input vector VS ) singular value x ) state vector y ) output vector Greek symbols R ) scalar reference input η ) vector of bounds on plant controllability δ ) vector of bounds on the distance between plant steady states λ ) vector of feedback controller parameters σ ) vector of scheduling variables Ω ) parametrization of the set of steady states Γ ) continuously differentiable curve of the reference variables domain Z(a, b) ) set of the continuously differentiable curves of the reference variables domain connecting a and b Subscripts 0 ) initial f ) final fb ) feedback ff ) feedforward max ) maximum min ) minimum N ) number of steady states p ) period s ) stationary Superscripts m ) dimension of the input vector n ) dimension of the state vector p ) dimension of the output vector l ) dimension of the controllability vector function Notation: Exothermic CSTR with Reactions A f B f C CA ) concentration of A CB ) concentration of B Cp ) specific heat (J/g · K) ER ) activation energy (J/mol) ∆H ) heat of reaction (J/mol) k ) kinetic constant (min-1) Q ) volumetric flow rate R ) gas constant (J/(mol K)) F ) density (g/L) T ) temperature (K) U ) heat-transfer coefficient (J/(min · K)) V ) reactor volume (L) Subscripts: Exothermic CSTR with Reactions A f B f C 1 ) reaction A f B 2 ) reaction B f C c ) coolant

9140

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

in ) inlet X ) reactant A or B

Literature Cited (1) McAuley, K., B.; MacGregor, J., F. Optimal grade transitions in gas phase polyethylene reactor. AIChE J. 1992, 38, 1564–1576. (2) Takeda, M.; Ray, W., H. Optimal grade transition strategies for multistage polyolefin reactors. AIChE J. 1999, 45, 1776–1793. (3) Zhang, S., X.; Read, N., K.; Ray, W., H. Runaway phenomena in low-density polyethylene autoclave reactors. AIChE J. 1996, 42, 2911– 2925. (4) Mancusi, E.; Merola, G.; Crescitelli, S.; Maffettone, P., L. Mulstistability and hysteresis in an industrial ammonia reactor. AIChE J. 2000, 46, 824–828. (5) Mayne, D., Q.; Rawlings, J., B.; Rao, C., V.; Skokaert, P., O., M. Constrained model predictive control: stability and optimality. Automatica 2000, 36, 789–814. (6) Chen, H.; Allgower, F. A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability. Automatica 1998, 34, 1205–1217. (7) Astrom K., J.; Wittenmark B. AdaptiVe Control; Addison: Wesley, MA, 1989. (8) Rugh, W. J,; Smamma, J., S. Research on gain-scheduling. Automatica 2000, 36, 1401–1425. (9) Bemporad, A. Reference governor for constrained nonlinear systems. IEEE Trans. Automat. Control 1998, 43, 415–419. (10) Gilbert, E.; Kolmanovsky, I. Nonlinear tracking control in the presence of state and control constraints: a generalized reference governor. Automatica 2003, 38, 2063–2073. (11) Chen, P., C.; Shamma, J., S. Gain scheduled l1-optimal control for boiler-turbine dynamics with actuator saturation. J. Process Control 2003, 14, 263–277. (12) Lauzze, K. C.; Chmielewski, D., J. Power control of a polymer electrolyte membrane fuel cell. Ind. Eng. Chem. Res. 2006, 45, 4661–4670. (13) di Bernardo, M. Bifurcation analysis for control system applications. In Bifurcation Control: Theory and Applications; Springer-Verlag, New York, 2003. (14) Wang, X.; di Bernardo, M.; Stoten, D., P.; Lowenberg, M., H.; Charles, G. Bifurcation tailoring via Newton flow-aided adaptive control. Int. J. Bifurcation Chaos 2003, 13, 677–684. (15) Richardson, T.; Lowenberg, M.; di Bernardo, M.; Charles, G. Design of a gain-scheduled flight control system using bifurcation analysis. J. Guidance Control Dyn. 2006, 29, 444–453. (16) Kwatny, H., G.; Chang, B., C. Constructing linear families from parameter dependent nonlinear dynamics. IEE Trans. Automat. Control 1998, 43, 1143–1147. (17) Gilbert, E.; Kolmanovsky, I. Set-point control of nonlinear systems with state and control constraints: A Lyapunov function reference governor approach. Proc. 38th Conf. Decis. Control 1999, 2507–2512. (18) Skogestad, S.; Morari, M. Effect of disturbance directions on closedloop performance. Ind. Eng. Chem. Res. 1987, 26, 2029–2035. (19) Skogestad, S.; Hovd, M.; Lundstrem, P. Towards integrating design and control: use of frequency-dependent tools for controllability analysis. Proceedings of the Process Systems Engineering, Montebello, Canada, August 1991. (20) Luyben, M., L.; Floudas, C. A. Analyzing the interaction of design and controls1. A multiobjective framework and application to a binary distillation column. Comput. Chem. Eng. 1994, 18, 933–969. (21) Seferlis, P.; Grievink, J. Process design and control screening based on economic and static controllability criteria. Comput. Chem. Eng. 2001, 25, 177–188.

(22) Sakizlis, V.; Perkins, J. D.; Pistikopoulos, E., N. Recent advances in optimization-based simultaneous process and control design. Comput. Chem. Eng. 2001, 28, 2069–2086. (23) Monnigmann, M.; Marquardt, W. Normal vectors on manifolds of critical points for parametric robustness of equilibrium solutions of ODE systems. J. Nonlin. Sci. 2002, 12, 85–112. (24) Monnigmann, M.; Marquardt, W. Steady state process optimization with guaranteed robust stability and robust feasibility. AIChE J. 2003, 49, 3110. (25) Monnigmann, M.; Marquardt, W. Steady-state process optimization with guaranteed robust stability and flexibility: Application to HDA reaction section. Ind. Eng. Chem. Res. 2005, 44, 2737–2753. (26) Grosch, R.; Monnigmann, M.; Marquardt, W. Integrated design and control for robust performance: Application to an MSMPR crystallizer. J. Process Control 2008, 18, 173–188. (27) Morari, M. Design of resilient processing plantssIII. A general framework for the assessment of dynamic resilience. Chem. Eng. Sci. 1983, 38, 1881–1891. (28) Perkins, J. D.; Wong, M., P. Assessing controllability of chemical plants. Chem. Eng. Res. DeV. 1985, 63, 358–362. (29) Kailath, T. Linear Systems; Prentice Hall: NJ, 1996. (30) Fletcher, R. Practical Methods of Optimization; Wiley: New York, 1994. (31) Uppal, A.; Ray, W., H.; Poore, A., B. On the dynamic behaviour of continuous stirred tank reactors. Chem. Eng. Sci. 1974, 29, 967–985. (32) Alhumaizi, K.; Elnashaie, S., E. Effect of control loop configuration on the bifurcation behaviour and gasoline yield of industrial fluid catalytic cracking (FCC) units. Math. Comput. Modell. 1997, 25, 37–56. (33) Hahn, J.; Monnigmann, M.; Marquardt, W. A method for robustness analysis of controlled nonlinear systems. Chem. Eng. Sci. 2004, 59, 4325– 4338. (34) Zavala-Tejeda, V.; Flores-Tlacuahuac, A.; Vivaldo-Limac, E. The bifurcation behavior of a polyurethane continuous stirred tank reactor. Chem. Eng. Sci. 2006, 61, 7368–7385. (35) Chen, G.; Moiola, J., L.; Wang, H. O. Bifurcation control: Theories, methods, and applications. Int. J. Bifurcation Chaos 2000, 10, 511–548. (36) Abed, E., H.; Fu, J., H. Local Feedback stabilization and bifurcation control I. Hopf bifurcation. Syst. Control Lett. 1986, 7, 11–17. (37) Wang, H. O.; Abed, E. H. Bifurcation control of a chaotic system. Automatica 1995, 31, 1213–1226. (38) Stilwell, D., J.; Rugh, W., J. Stability preserving interpolation methods for the synthesis of gain scheduled controllers. Automatica 2000, 36, 665–671. (39) Fernandez Anaya, G.; Flores-Tlacuahuac, A. Interpolated controllers for the robust transition control of a class of reactors. AIChE J. 2006, 52, 247–254. (40) Rugh, W. Analytical Framework for Gain Scheduling. IEEE Control Syst. Mag. 1991, 11, 79–84. (41) Lawrence, D.; Rugh, W. Gain scheduling dynamic linear controllers for a nonlinear plant. Automatica. 1995, 31, 381–390. (42) Gamboa-Torres, A. E.; Flores-Tlacuahuac, A. Effect of process design/operation on the steady state operability of a CSRT-reactions Af BfC. Chem. Eng. Res. Des. 2000, 78, 481–491. (43) Seydel, R. Practical Bifurcation and Stability Analysis: From Equilibrium to Chaos; Springer-Verlag: New York, 1994.

ReceiVed for reView January 29, 2009 ReVised manuscript receiVed July 23, 2009 Accepted August 8, 2009 IE9001553