445 1
Ind. Eng. Chem. Res. 1995,34, 4451-4462
Flexibility Analysis of Dynamic Systems Veniamin D. Dimitriadis and Efstratios N. Pistikopoulos* Centre for Process Systems Engineering, Department of Chemical Engineering, Imperial College, London SW7 2BY, UK
This paper presents a unified approach for the quantification of feasibility and flexibility of systems that operate dynamically under the effect of time-varying uncertainty. Based on the steady state flexibility analysis framework, the formulation for both the dynamic feasibility and dynamic flexibility analysis problems is developed. It is shown that when the assumption of given uncertainty profile holds, both problems reduce to dynamic optimization problems that can be solved using currently available techniques. For the case when this assumption is not valid, a discretization scheme for the daerential equations is combined with a n active constraint strategy to transform both problems into mixed integer nonlinear programming problems. A number of examples are used to illustrate the applicability of the proposed framework.
Introduction
equalities and inequalities of the following form:
Chemical plants are complex dynamic systems that have to operate succesfully throughout time in the presence of uncertainties. These uncertainties are, in general, also dynamic in nature and correspond to variations in either external variables (such as quality and quantity of feedstocks, product demands, etc.) or internal process parameters (such as heaumass transfer coefficients, kinetic constants, etc.). Under these circumstances, flexibility, i.e., the ability to maintain feasible operation over a range of uncertain conditions, is clearly a characteristic of key importance and has to be explicitly incorporated into the design objectives along with any profit-based terms. This can only be done if a measure exists on which the comparison of different design alternatives with respect to flexibility of operation can be based. The determination of such a measure establishes the flexibility analysis problem (Grossmann and Straub, 1991). The flexibility analysis problem generally consists of two tasks which are complementary t o each other. (i) The first task is determining if a given design can feasibly operate over the range of uncertainty considered. This problem is known as the feasibilityproblem or flexibility test problem. (ii) The second task is calculating a measure to quantify the ability of the design to operate in the presence of uncertainty. This is known as the flexibility index problem and is usually tackled by establishing the maximum parameter range over which the design can operate feasibly. During the past decade, a lot of research work has focused on providing the necessary theoretical background for the flexibility analysis of chemical processes and incorporating flexibility objectives into the corresponding design procedures. Most of this work has concentrated on processes that operate at steady state. Such processes can be described by sets of algebraic
* Author to whom all correspondence should be addressed. Telephone: f 4 4 171 5946620. FAX: +44 171 5946606. E-mail:
[email protected].
zEZ={ZIZ
L
I Z I Z
U
}
(4)
where dim{h} = dim{x}. In this description, d is the vector of the design variables, x is the vector of the state variables, z is the vector of the control variables and 8 is the vector of the uncertain parameters. The design variables correspond to decisions that are taken a t the design stage ( e g . , process configuration, equipment sizes, etc.),and therefore their values remain unchanged during the operation of the process. The control variables are the degrees of freedom that are available during process operation and can be adjusted for different values of the uncertainty vector 8 in order to drive the system state to desirable points. The state variables can be implicitly eliminated between (1)and (2),leading to the following simplified process description:
Halemane and Grossmann (1983) showed that the feasibility problem for processes that are described by (5) is equivalent to the following max-min-max optimization problem:
where J is the index set for the inequalities. If x(d) I 0, the design is feasible in T. If, on the other hand, x(d) > 0, the solution will identify a critical point 8c E T for which the greatest violation in the constraints occurs. Swaney and Grossmann (1985a) extended this work and proposed a scalar flexibility index, F, to quantijj, the ability of such processes to operate at points other than the nominal point of operation. The determination of F requires the solution of the following
0888-5885/95/2634-4451$09.00/00 1995 American Chemical Society
4452 Ind. Eng. Chem. Res., Vol. 34, No. 12, 1995
optimization problem:
(FI)
F=max6 subject to
max min maxfj(d,z,B)I0 BcT(6) ZEZ j d
6 2 0
T(6)= {$I@ - 6A8-
SN + 6A8')
:I/
I8 I
Qualitatively, F represents the largest scaled deviation that the design can accommodate, while remaining feasible. Deviations are measured from a given nominal value of the uncertain parameter vector, @, and are scaled according to given expected deviations from this value, he- and A8+. Swaney and Grossmann (1985b) also showed that, under certain convexity assumptions, critical points that limit feasibility and/or flexibility lie on the vertices of the uncertainty space. Grossmann and Floudas (1987) exploited the fact that sets of active constraints are responsible for limiting the flexibility of a design and posed the problem of determining these active constraints as a mixed integer linearhonlinear programming problem. Their formulation does not assume that the critical points are vertices of the uncertainty space. Later, Pistikopoulos and Mazzuchi (1990) introduced the stochastic flexibility index, SF,t o quantify the flexibility of processes with stochastic parameters. The stochastic flexibility index is defined as the integral of the multivariate joint probability distribution function of the uncertain parameters over the feasible region of operation of the process. As was noted before, the common assumption behind the previous approaches has been that the process operates at steady state for any value of the uncertain parameter vector. Although this assumption is true for most of the operating life of a continuous process, there are certain situations in which the results of the steady state flexibility analysis are of little help. Such situations include startup or shutdown of a process, transitions from one operating point to another, actions taken t o eliminate or minimize the effect of disturbances, and so on (Grossmann and Morari, 1983). Furthermore, the steady state analysis cannot be applied directly to batch processes because of their inherently dynamic nature. Morari (1983) emphasized the need for a systematic method to assess the dynamic operability characteristics of chemical processes and incorporate such considerations a t the design stage. He proposed a general framework for the evaluation of the dynamic resilience of process plants that relies on the fact that certain system-related factors pose a limit on the achievable closed-loop behavior of a process. Methods for estimating the effect that these factors have on the dynamic resilience of the system were presented by Morari (19831, Holt and Morari (1985a,b), and Skogestad and Morari (1987). They are restricted to process models that are linearized around a specified steady state operating point and can be used to provide a basis for comparison at the early design stages. Soroush and Kravaris (1993a,b) addressed the problem of flexibility of operation for batch reactors. For the purposes of their analysis, flexibility is defined as the ability of a reactor to operate according to a predetermined optimal way in the presence of parametric uncertainty. However, this definition is rather conservative as it does not account for the fact that process operation may be adjusted to compensate for
Figure 1. Storage tank dynamic system.
changes in the values of the uncertain parameters. Furthermore, their analysis does not provide a quantitative measure of flexibility and cannot be directly applied to the case of time-dependent uncertainty. The effect of uncertainty on the dynamic behavior of chemical processes was also studied by Walsh and Perkins (19941,with particular reference to wastewater neutralization processes. The identification of critical points for a given design is done using the heuristic that such points involve most uncertain parameters lying at their upper or lower bound. Once the possible critical points are identified, the performance of the system for each one of them is evaluated by integration. More recently, dynamic optimization techniques were also used to tackle the feasibility problem of dynamic systems. White et al. (1994)presented an approach for the evaluation of the switchability of a proposed design, i.e., its ability to perform satisfactorily when moving between different operating points. Dimitriadis et al. (1994)considered the feasibility problem from the safety verification point of view. Given a mathematical model for the process, their formulation can be used t o determine if there is a combination of events, from within the uncertainty space considered, that can lead the system to any one of a number of predefined dangerous operating regions. All these approaches share the limiting assumption that the control scheme is completely determined before the analysis step and hence the system model has no degrees of freedom. This paper addresses the problem of determining the feasibility and flexibility of operation for dynamic systems, i.e., for systems that are described by sets of differential and algebraic equations and are subject to time-varying uncertainty. Based on the steady state approach by Swaney and Grossmann (1985a), a formulation for the dynamic feasibility problem is presented and a dynamic flexibility index is defined. It is shown that, when the variation of the uncertain parameters over time is known or when the critical points that limit feasibility are vertices of the time-varying uncertainty space, both the feasibility problem and the flexibility index problem transform into dynamic optimization problems that can be solved using established techniques. Finally, when neither of the previous assumptions is valid and the system is subject to time-varying uncertainty with unknown time profile, an active constraint strategy is proposed for the identification of the critical uncertainty profile that relies on the explicit discretization of the time horizon. A number of example problems are used to illustrate the applicability of the proposed approach.
Dynamic Feasibility/FlexibilityAnalysis To motivate the problem to be addressed in this work, consider the simple storage tank dynamic system depicted in Figure l. The volume of liquid in the tank,
Ind. Eng. Chem. Res., Vol. 34, No. 12, 1995 4463 Table 1. Mathematical Model for Motivating Example
(ii)The second category is point constraints that have to be satisfied a t certain time instances within the period of operation:
gjl”int(d,x(tk),z(tk),8(tk),tk) I0, 12 = 1, ...,NP V ( t ) ,is the variable that defines the state of the system at any time. The input flow rate, Fo(t),is the uncertain parameter and can vary with time, taking any value between a lower and an upper bound. The output flow rate, F(t),is a linear function of the volume of liquid in the tank. The proportionality parameter, a@),either can vary with time in a predetermined manner or can be regarded as a control variable that can be used to compensate against changes in the input flow rate. The first case occurs when the system is uncontrolled or when the type of control action is predetermined and incorporated in the system model a t the design stage. The second case occurs when a is an extra degree of freedom in the system model and its behavior with time remains to be determined. The dynamic model of the system is shown in Table 1. For the system to operate normally, the volume of liquid in the tank should lie between certain limits. This requirement defines the inequality constraints that determine the feasibility of the system a t any time:
P(t)IV(t)IP(t) It is clear that steady state flexibility analysis tools cannot be applied t o problems of this kind, i.e., to processes that have an inherently dynamic nature and/ or are subject to dynamic uncertainty and feasibility conditions. In the sequel, a formulation for the feasibility problem will be proposed that can be used t o establish whether the system can feasibly operate throughout a specified time horizon, tolerating any possible profile of the uncertain parameter, Fo(t). Furthermore, a dynamic flexibility index will be defined in order to quantify the ability of the system to dynamically tolerate timedependent variations in the value of the uncertain parameter. The physical performance of a dynamic chemical process can generally be described by a set of ordinary differential and algebraic equations of the following form:
h(d,t(t),x(t),z(t),e(t),t) = 0 ; x(0) = xo
Two common types of point constraints are initialtime and end-time point constraints, i.e.,constraints that have to be satisfied a t the beginning and at the end of the period of operation. For example, in a batch reactor, the requirement that the h a l conversion of the reactant should be greater than a lower bound is an end-time point constraint. Obviously, constraints 7 determine the feasibility of the process, i.e., its ability to operate successfully and accomplish its stated objectives. Equations 6 and 7 are very similar to eqs 1and 2 that are used by Swaney and Grossmann (1985a) to describe the behavior of steady state processes, the only difference being that the variables describing the operation of the process now depend on the time variable, t. This implies that the feasible region of operation of .the process and the uncertainty space change with time. In the steady state case, for fured design d and accounting for the action of the control variables z, a feasible region can be defined in the 8-space, constrained by combinations of the fi constraints. The feasibility of this combination of design and control variables is then determined by whether or not this feasible region contains the given uncertainty space or not. In the dynamic case, however, the definition of the feasible region and the uncertainty space additionally depend on the actual time considered, because z and 8 are now functions of time. This is schematically illustrated in Figure 2, where it can be seen that a system can be feasible a t t = tl and infeasible at t = tz. In order t o model this, a time horizon, H , is defined which corresponds t o the maximum time for which the feasibility of the system is of interest. Then, the dynamic feasibility problem for this system may be stated as follows. Establish if, for every possible profile of the uncertain vector @ t ) , there is at least one profile of the control vector z ( t )such that the process feasibility constraints are satisfied at any time t E [O,K1. This statement is mathematically equivalent t o the following max-min-max optimization problem (Halemane and Grossmann, 1983):
(6)
where dim{h} = dim{x}. In the previous equation, d is the vector of the design variables that remains constant at the operation stage and x ( t ) , zW, and e ( t ) are the time-dependent vectors corresponding to the state variables, control variables, and uncertain parameters of the process, respectively. A number of constraints are always imposed on the operation of the process. These constraints can be classified into two categories: (i) The first is path constraints that have to be satisfied throughout the period of operation:
gpath(d,x(t),z(t),e(t),t) I0
(7b)
(7a)
For instance, in the motivating example, the requirement that the volume of liquid in the tank should lie between certain bounds during operation is a path constraint.
subject to
h(d,t(t),x(t),z(t>,e(t),t) = 0; x(0) = xo ~(= t ){e(t)leL(t) Ie(t)Ie%)} Z(t)= {z(t)lzL(t) 5 z(t) 5 zU(t)}
where J is the index set for the feasibility constraints (for notational simplicity, both path and point constraints have been combined into one feasibility constraints vector, g). Problem DFTl is the dynamic feasibility problem or dynamic flexibility test problem. If x(d) I0, the proposed design is dynamically feasible in T ( t ) . Otherwise, there exists a t least one possible profile, 8 ( t ) ,for which no corresponding control action, z(t),can be found to keep the system feasible at all times within the horizon considered. Problem DFTl can be transformed into a bilevel optimization problem by
4454 Ind. Eng. Chem. Res., Vol. 34, No. 12, 1995
a time horizon, H , the dynamic flexibility index evaluation problem can be formulated as follows:
DF(d) = max 6
(DFI1)
subject to X(d) = max min rnax g,(d,x(t>,z(t),e(t),t) I0 B(t)sT(d,t)z(tl&(t) jd,telO,K
+
~ ( 6 , t=) {e(t)leN(t)- sae-(t)Ie(t)IeN(t)
I
dAO+(t)) Z(t) = {z(t)lzL(t)Iz(t) IzU(t)) Qualitatively, the dynamic flexibility index, DF, represents the largest scaled deviation of the uncertain parameter profile that the design can tolerate while remaining feasible within the horizon considered. Deviations are measured from a nominal profile, eN(t),and are scaled using the expected deviation profiles Ae-(t) and A@%). Normally,
@(t) - A W ( t ) = &(t) 6%)
i
I
* ai
Feasible region boundary - Uncertainty space boundary ~LLLL
Figure 2. Variation of the feasible region of operation with time. (a)Design is feasible at t = tl for all values of 01 and 82. (b) Design is infeasible at t = t 2 for some values of 01 and 02.
introducing the feasibility function Y(d,O(t)):
+ AO+(t)= 6%)
Both the flexibility test problem DF"2 and the flexibility index problem DFIl are two-stage, semi-infinite, dynamic optimization problems with an infinite number of decision variables. In the next section, problem DFT2 will be simplified by assuming that the uncertainty profile, e(t),is given. Problem DFIl will also be simplified by assuming that the critical direction, A W ) , from is known. Both the nominal uncertainty profile, eN(t), assumptions will be relaxed in the last section of the paper.
Feasibility/FlexibilityAnalysis with Given Uncertainty Profile subject to
Y!(d,e(t))=
min
u
u,z(t)€Z(t),t€Iom
The dynamic feasibility problem DFT2 can be greatly simplified by assuming that the variation of the uncertainty vector with time, e@),is given. This is the case when either the exact profile is known or the critical points that limit feasibility are assumed to lie on the vertices of the time-varying uncertainty space. Under this assumption, the outer optimization vanishes and problem DF"2 reduces to a dynamic optimization problem of the following form:
x(d)=
min
u
(DFTS)
u,z(t)eZ(t),tdO,HI
As was mentioned before, the feasibility problem determines the ability of a system to operate feasibly, withstanding variations of the uncertain parameters within a given space of possible values. On the other hand, the flexibility index evaluation problem rigorously quantifies this ability, by establishing the maximum uncertainty range over which feasibility can be maintained. In the case of dynamic processes, however, the feasible region of operation, and consequently the value of the flexibility index, change with time. By defining
g(d,x(t),z(t),e(t),t)Iu
Z(t>= (z(t)lzL(t)5 z(t) IzU(t)} The dynamic flexibility index evaluation problem can also be transformed into a dynamic optimization problem by assuming that the direction, A&(t), from the nominal profile, W ) ,on which the critical point that limits flexibility lies is either known or is assumed to
Ind. Eng. Chem. Res., Vol. 34, No. 12, 1995 4455 be one of the vertex directions of the time-varying uncertainty space. Then, problem DFIl takes the following form:
max
DF(d) =
6
(DFIB)
I 40
I 20
G,z(t)eZ(t),t~[O,Hl
subject to
h(d,%(t>,x(t>,z(t),e(t),t) = 0; x(0)= xo g(d,x(t>,z(t>,W,t> I0
+
O ( t ) = 0%) 6AeC(t)
100
0 80
0 60
620 0 40
Z(t)= {z(t)lzL(t)Iz(t) IzU(t)} Practical solution algorithms for dynamic optimization problems usually fall into one of the following categories: (i) The first is full discretization algorithms, where all the equations in the problem are discretized with respect to all variables. The resulting problem is a large-scale nonlinear programming problem. (ii) The second category is control parameterization algorithms, where only the control vector is parameterized using a finite set of parameters. The result is again a nonlinear programming problem, but here the optimization is coupled with direct integration that is used to derive the gradients of the objective function and constraints with respect to the parameters. In this work, dynamic optimization problems were solved using the DAEOPT optimal control code (Vassiliadis et al., 1994a,b). DAEOPT uses the control parameterization approach. It partitions the time horizon into a specified number of control intervals of variable length and utilizes low order polynomial functions to approximate the control variables in each interval. Equality and inequality path constraints are dealt with efficiently by using point constraints as well as forcing an integral measure of their violation to zero. Revisiting the motivating example of Figure 1, assume that the following data are given:
vo = 1.0 P(t)= 0.9,
V(t) = 1.1
The uncertain parameter F&) varies according t o the following equation:
0.05, t = 0
0.05, t
=- 5
First, consider the case where the proportionality parameter a is fxed ( a = 0.05); i.e., no control over the behavior of the system exists. The solution of the dynamic feasibility problem DFT3 for a time horizon H = 60 is x(d) = 0.342 > 0. The evolution of the system state with time is shown in Figure 3, where it can be seen that the system becomes infeasible a t approximately t = 1. Now assume that a can be used as a control variable with
aL(t)= 0.00, au(t>= 0.14 Additionally, assume that a can follow a piecewise linear behavior over 10 control intervals. In this case
0 20
Lr
1 i
0 00
Figure 3. Feasibility test for motivating example: uncontrolled case.
problem DFT3 results in x(d) = -0.064 (Figure 4). This shows that the system is now dynamically feasible, due to the fact that the control variable can be used to compensate succesfully for the uncertain parameter variation. Although the obvious choice of control variable for this system is a, in practical situations it may be unrealistic to assume that the installed control system can impose such a perfect behavior on the valve constant. If a certain type of controller, eg., a PI controller, is installed on the process, the controller equation has to be added to the process model:
r(t>= V(t)- v It can be seen that the degree of freedom represented by a is now taken up by the control loop. However, the system is not completely determined until the values of the tuning parameters, K, and t,are fxed. In fact these parameters are now the degrees of freedom (control variables), and the dynamic feasibility formulation can be used to determine if there are values for Kc and z such that the system is feasible for the uncertain parameter variation considered. The solution t o problem DFT3 with 0.0 IK,I5.0 1.0 It I30.0
is x(d) = -0.013 and the corresponding values for the tuning parameters are K, = 0.922 and z = 30.0. The behavior of the system is shown in Figure 5. It should be noted that the only criterion used for the selection of K,and z is feasibility of the system. In practice, once feasibility is guaranteed, one can use another objective function for the problem, eg., an integral error objective function, in order to tune the controller not only for feasibility but also for optimal system response. Consider again the uncontrolled case of the motivating example problem and assume that the nominal
4456 Ind. Eng. Chem. Res., Vol. 34,No. 12, 1995
,
I
I
I
100
IW
V
I
0 80
0 60
!
0 80
0 60
o 40 0 40
0 20
0 20
FO 0 00
Tune
000
000
20.00
ow
20 00
40
w
0 14
0 14 0 17
0 13 0 I2 0 I2 11 I I
0II
0 IO
0 io 0 09
0 00 0 ox 0 Ob 0 07
0 07 0 o(, 0 06
0 os
Time
0 OT
Tune 0 00
20 00
40.00
60.00
Figure 4. Feasibility test for motivating example: perfect control case.
profile and possible deviations of the uncertain parameter are the following:
6%) = 0.05, AO'(t) = 0.15, AO-(t) = 0.00 Also assume that the possible critical directions correspond t o vertex directions from this nominal point. For the positive direction (AO(t) = 0.1) the solution to problem DFI2 is DF+ = 0.054, whereas for the negative direction (AO(t) = -0.05) the solution is DF- = 0.106 (the time horizon used was H = 60 for both cases). The dynamic flexibility index for the system is the smallest of these values, z.e., DF = DF+ = 0.054. This result shows that the maximum variation in the uncertain parameter that the system can tolerate within the time horizon considered is (0.054)(0.1) = 0.0054 in the positive direction. If the same analysis is repeated for the controlled case the results are DF+ = 1.042, DF- = 1.000, and DF = 1.000. It is clear that the flexibility of the system increases drastically when the control variable becomes available. Clearly, the results of both the feasibility and flexibility index problems depend on the length of the time horizon as well as the number of control intervals used
40 00
60 00
Figure 6. Feasibility test for motivating example: PI control case.
for the solution of the corresponding dynamic optimization problems. By solving the feasibility and/or the flexibility index problems for different values of H, we can construct a profile of either feasibility or flexibility over time. In Figure 6, the evolution of the dynamic flexibility with time is shown for the controlled case of the motivating example problem. It can be seen that the dynamic flexibility index for the negative direction remains constant at its maximum possible value, DF= 1.000. The dynamic flexibility index for the positive direction has a big initial value and decreases with time, asymptotically reaching a value of DF+ = 1.042. This is logical if we consider that, as the length of the time horizon increases, the demand posed on the control system increases accordingly and therefore the flexibility of the process decreases. Two example problems will be considered below t o illustrate the applicability of formulations DFT3 and DFI2. The first problem is concerned with the operation of a batch reactor in which the heat of reaction is an uncertain parameter. The second problem examines the operation of a heat exchanger network under fouling. Example 1. A highly exothermic, first-order reaction (A B) is carried out in a jacketed batch reactor. The mathematical model that is used to describe the system is given in Table 2. The heat of reaction, AHr, is not
-
Ind. Eng. Chem. Res., Vol. 34,No. 12,1995 4457 Dynanuc Flexibility Index
170
1.50
Cooling water flowrate (kg/s)
H
-I
1
1
1l l
05 ,00
........................................
O0
i
1
1 Time(hr)
0 00
0 50
IW
I so
2 00
Figure 7. Cooling water flow rate vs time at critical point for example problem 1. Reactant holdup (mol)
20000
'4000 I2000
1 \
::I 8000
i
L-
i -i
,
,
,
0 50
I 00
I50
-
2000
0
(ii) Due to safety considerations, the temperature in the reactor must not exceed a maximum of 400 K at any time during operation, i.e.,
T(t) I400 K, 0 It IH (iii)The temperature of the reacting mixture at the end of the operation must not exceed 308 K, so that discharge and cleaning operations can take place, Le.,
T(t=H) I308 K Given that the duration of each batch cycle is 2 h, it is desirable t o determine if the reactor can accommodate the whole range of uncertainty while operating feasibly according to the aforementioned constraints. Note that the uncertainty considered in this example is parametric and hence time-invariant. However, steady state flexibility analysis techniques cannot deal with the problem directly because the process is inherently dynamic. The dynamic flexibility index evaluation problem DFI2 is solved for the two vertex directions from the nominal point of operation. The time horizon is divided into six intervals of variable length and the profile of
000
2 00
Figure 8. Reactant holdup vs time at critical point for example problem 1.
the cooling water flow rate, Fcw, is assumed to be piecewise constant over these intervals. Along the positive vertex direction, it is found that the system can handle 78% of the expected uncertainty range (DF+ = 0.7788, = 81 394 J/mol). Similarly, along the negative vertex direction, DF- = 7.2384 and = 41 308 J/mol. The dynamic flexibility index of the system is therefore DF = 0.7788,and the critical point of operation corresponds to a value of M,T't= 81 394 Jlmol. The required profile of the cooling water flow rate for the system to operate feasibly a t this critical point is shown in Figure 7. In this case the process constraints on final conversion and final reactor temperature are marginally satisfied as shown in Figures 8 and 9. Example 2. The heat exchanger network depicted in Figure 10 is considered in this example. It is taken from Kotjabasakis and Linnhoff (19871,and it involves two hot streams and two cold streams. Stream data and
e+
4458 Ind. Eng. Chem. Res., Vol. 34, No. 12, 1995 H1
Reactor temperature (K)
!
3Y0 W
H2 2UOOC
l55OC
h
60 kWiK
M)
kWK
380 00
370.U
I
L-
160 00 350 00
340 00
130.00
-
w
-
32n
i
\< oooc
(115°C
Figure 11. Retrofitted heat exchanger network for example problem 2.
-
I I
I
I
n 50
o no
Io0
1
I
150
2 00
Time (hr)
Figure 9. Temperature vs time at critical point for example problem 1. Hl
n2
s I lS"C
( ,
90°C
Figure 10. Existing heat exchanger network for example problem 2.
Table 3. Stream Data for Example Problem 2 stream TIN("C) Tom ("C) FCp (kWiK) H1 155 590 60 H2 200 5115 60
c1 c2
steam
20 20 500
2175 5 175 500
20 50
Table 4. Exchanger Data for Example Problem 2 existing retrofitted exchanger UN (kW/m2K) area (m2) area (m2) 1 2 3 4
0.225 0.120
0.150 0.360
65 305 600 175
170 305 600 175
exchanger data for the network are given in Tables 3 and 4, respectively. The network can operate feasibly when the heat transfer coefficients have the nominal values shown in Table 4. However, due to exchanger fouling, these values decrease with time according t o an asymptotic law of the following form:
U = [RI + RX1 The values of the constants RI,Rf,and t for each exchanger are given in Table 5. Given that exchanger
Table 5. Fouling Data for Example Problem 2 exchanger
RI (m2KkW)
Rf(m2KkW)
t (months)
1
4.44 8.33 6.67 2.78
0.45 5.32 0.48 0.23
1.48 4.27 0.62 0.39
2 3 4
maintenance takes place every 6 months of operation, it is desirable to determine if the network can operate feasibly throughout this period. It is assumed that the available control variables are the steam heat load in the heater and a bypass flow rate for each of the heat exchangers. The time horizon is divided into six control intervals, and each of the aforementioned control variables is taken to be piecewise constant over these intervals. The result of the dynamic feasibility problem DFT3 for this case is X(d) = 12.513 > 0, implying that the network cannot operate feasibly under the effect of the uncertainty considered. A retrofit attempt to increase the flexibility of the existing network is shown in Figure 11 (Papalexandri and Pistikopoulos, 1993). Note that, in the retrofitted network, the area of exchanger 1 has been increased and only two bypass streams are used. The solution of the dynamic feasibility problem is now x(d) = -1.444 < 0, thus showing that the new network can indeed operate feasibly. Once the feasibility of the proposed design is determined, it is possible to set the infeasibility variable, u , equal to zero and solve problem DFT3 again with an economic objective function, e g . , minimization of the total utility consumption over the 6 months of operation, in order to additionally accomplish optimality of operation. The solution indicates a total utility consumption of 7.11 x lo6kW h, and Figure 12 is the profile of energy use over time. As expected, the utility consumption increases with time as fouling gets worse.
An Active Constraint Strategy for Feasibility/ Flexibility Analysis of Dynamic Systems As mentioned before, the dynamic feasibility analysis problem DFT2 is a two-stage, semi-infinite, dynamic optimization problem with an infinite number of decision variables. This problem is very similar to the steady state feasibility analysis problem FT that was addressed by Grossmann and Floudas (1987). The solution procedure they proposed exploits the fact that
Ind. Eng. Chem. Res., Val. 34, No. 12, 1995 4459
6 and 7 can be transformed into the following algebraic form:
Utility consumption (kW)
h;,(d,x,z,O,t,a) = 0 dm,(d,x,z,8,t)I0 where m E M,n E N , i E I, and j E J are the index sets for the finite elements, collocation points, differential equations, and feasibility inequalities, respectively. Note that the time-dependent state, control, and uncertainty vectors have been parameterized and expressed using a finite number of coefficients. These coefficients are the elements of the time-invariant vectors x, z, and 8, respectively. t and a are the time locations of the collocation points and the lengths of the finite elements, respectively. Given this transformation, the inner optimization of the dynamic feasibility analysis problem, DFTB, is rewritten as follows: 0.00
2 00
4 00
min u
6.00
Figure 12. Profile of energy use over time for example problem 2.
subject to
the feasibility of opeiation of a steady state system, for a given value of the uncertainty vector, 8, is characterized by a number of limiting or active constraints. These constraints are a subset of the system feasibility constraints, f i , of dimensionality n, 1, where n, is the number of available control variables. The set of active constraints and the optimal values of the control variables for each possible realization of 8 are determined by expressing the solution of the inner optimization problem as a square system of equations and introducing binary variables, yj, to denote if feasibility constraint f i is active or not. The result is that the twostage optimization problem is reduced to a single-stage mixed integer linearlnonlinear programming problem. There are two main difficulties in extending the aforementioned ideas to the solution of the dynamic feasibility analysis problem: (i) In the dynamic case, the system is described by a set of differential and algebraic equations, as opposed t o the set of algebraic equations of the steady state problem. (ii)When the two-stage problem is reduced to a singlestage equivalent, the profiles of the control variables and the uncertain parameters over time, z ( t ) and 8(t), respectively, represent an infinite number of decision variables. Both problems can be overcome by converting the differential equations t o algebraic residual equations, thus reducing the dynamic optimization problem to a nonlinear programming problem. An accurate and stable method to accomplish this reduction is by using orthogonal collocation on finite elements. According to this technique, which has been used extensively in the literature for the solution of dynamic optimization problems ( e g . , Biegler, 1984; Cuthrell and Biegler, 19871,the profiles of the state and control variables are approximated in each of a number of finite elements by Lagrange polynomials. The approximation is based on a number of collocation points that are usually the roots of a class of orthogonal polynomials. Continuity of the state profiles is enforced a t the beginning of each element, and the element lengths can be optimization variables. The details of the method are briefly discussed in the Appendix, where it is also shown how eqs
+
hAn(d,x,z,8,t,a)= 0
m E M , nEN, iEI, j E J The solution of this optimization problem can be determined by solving a system of equations, derived from the Kuhn-Tucker optimality conditions and the fact that the number of active constraints a t the solution must be equal to n, 1, where n, = dim{z} is the number of unknown coefficients that are used to parameterize the control profiles over time. In order to express the second requirement, a set of binary variables, y',,, is introduced to indicate if the feasibility constraint dinis active or not. Under this definition the Kuhn-Tucker optimally conditions are the following:
+
hin(d,x,z,8,t,a)= 0
dmn - U(1- f i n )
I0
4460 Ind. Eng. Chem. Res., Vol. 34, No. 12, 1995
Y i n E { 0 , 1 > , Ai,, mEM, nEN,
din
2
0
iEI, j E J
Ai,
where din are slack variables and p i , and are the Lagrange multipliers for the equalities and inequalities, respectively. The last two constraints imply that if a feasibility inequality is active then its corresponding slack variable is set to zero, whereas if it is inactive then its correspondingLagrange multiplier is set to zero. The previous set of equations can be used to determine the optimal values of z and a for a given realization of the uncertainty 8 and can be incorporated as a set of constraints in the outer optimization problem of DFT2 to give the following mixed integer nonlinear programming problem:
max u
the maximum deviation from feasible operation under the best achievable parameterized control action. The generalized dynamic flexibility index evaluation problem is formulated as a mixed integer nonlinear programming problem using exactly the same concepts and is presented below: min 6 subject to
(DFI3)
hk,(d,x,z,B,t,a) = o
din + din(d,x,z,8,t)= 0
(DFT4)
hLn(d,x,z,8,t,a)= 0
subject to
din + g',,(d,x,z,B,t) ah;, pi,-+ m,n,i
-
u =0
admn
Ain-m,nj
8x
-0
C Yi,,= dim{z} + 1
C Yin= dim{z} + 1
Y'n
11,
ALn, dmn 2
0,
620
m E M , nEN, i E I , j c J
where the elements of the vector 8 are appropriately bounded to reflect the uncertainty space considered. This is the generalized formulation for the dynamic feasibility analysis problem. Note that, in this formulation no specific knowledge about the uncertainty profile is assumed except that the time-dependent uncertainty is parameterized' to make the number of decision variables finite. For example, the uncertainty may be (i) piecewise constant over the finite elements considered, (ii) piecewise linear (with or without continuity) over the finite elements considered, or (iii) any function of prespecified form in which one or more parameters remain to be determined, eg., a sinusoidal function with variable amplitude and/or frequency. The global solution of the problem will identify the parameterized uncertainty profile over time that causes
Consider again the motivating example of Figure 1. If the dynamic equation that describes the behavior of the system is discretized over 10 finite elements of variable length, each containing 2 collocation points, the evolution of the system state over time, V(t),is approximated by a continuous, piecewise quadratic function. Furthermore, the assumption of known uncertainty profile is relaxed and both the uncertainty profile, F&), and the control variable profile, a(t), are now represented by piecewise linear functions over the finite elements considered, with 0.00 Ia(t)I0.14 0.00 IF J t ) I0.20
The dynamic feasibility problem for this simple system is a mixed integer nonlinear programming problem involving 333 continuous variables, 143 binary variables, and 423 constraints. This problem was solved via the Generalized Benders Decomposition procedure (Geoffrion, 1972), using as a starting point the uncertainty profile of Figure 3. The solution of the problem is x(d) = -0.029, implying that the system is dynamically feasible for the critical realization of the uncertainty. The critical uncertainty profile together with the
Ind. Eng. Chem. Res., Vol. 34, No. 12, 1995 4461 .....
Element 1
110
Element m
, ,
...
Element M
1
I
Io0 0 90
*
L -
a,"
0 15
7
A
%+I T= I
5=0
Collocation point n
Figure 14. Orthogonal collocation on finite elements. 0 14
a
Figure 13. Feasibility test for motivating example: critical profiles.
corresponding state and control variable profiles are shown in Figure 13. The convexity conditions that must hold for global solution in both DFT4 and DF13 are the same as the ones presented by Grossmann and Floudas (1987). However, it must be noted that these conditions will usually not hold, especially if the finite element lengths are treated as optimization variables to be determined. Therefore, unless a global optimization algorithm is used, it is very difficult t o rigorously prove the global optimality of the solution in the general case.
Conclusions A systematic framework has been proposed in this paper for the assessment of feasibility and flexibility of dynamic systems. The dynamic feasibility and flexibility problems were formulated as two-stage, semiinfinite, dynamic optimization problems with an infinite number of decision variables. If the uncertainty profile over time is given, both problems reduce t o dynamic optimization problems that can be solved using currently available techniques. This was illustrated using two example problems. The same reduction is possible when the system model is completely determined, i.e., if no control variables are available (Dimitriadis et al., 1994). For the case when neither of these assumptions can be made, a discretization scheme for the differential equations is combined with an active constraint strategy to transform both problems into mixed integer nonlinear programming problems. The proposed approach extends the steady state flexibility analysis framework t o incorporate the time element and has important applications in the areas of operability analysis, control structure selection, and controller tuning for dynamic chemical processes (Mohideen et al., 1995). However, a number of computational issues remain to be resolved. The size of the optimization problems is currently rather big, even for small problems. This is mainly due to the fact that the constraints have to be written over all discretization points. More effective discretization techniques combined with decomposition schemes where the state variables and the corresponding equations are handled separately (e&., Logsdon and Biegler, 1992) should be worth exploring.
Nomenclature d = vector of design variables x = vector of state variables z = vector of control variables t = vector of time locations for collocation points u = infeasibility variable s = slack variables t = time i = index set for equality constraints j = index set for inequality constraints m = index set of finite elements n = index set for collocation points a = vector of finite element lengths 6 = flexibility index variable 0 = vector of uncertain parameters 1,p = Lagrange multipliers
Appendix Orthogonal collocation techniques apply a polynomial approximation to the differential equation and require satisfaction of the approximating equation at a number of discrete collocation points which are usually the zeros of a class of orthogonal polynomials. In addition, orthogonal collocation on finite elements divides the time horizon into a prespecified number of finite elements and uses a separate approximation for each element. Continuity of the state profiles at the beginning of the elements is enforced. The situation is depicted in Figure 14. The time horizon is divided into a number of finite elements m = 1,...,M. The state, control, and uncertainty profiles x(t), z(t),and 0(t),respectively,are approximated within each element by Lagrange polynomials, based on a number of collocation points n = 1,..., N. The relative location of the collocation points in each element is the same and corresponds to the roots of a class of orthogonal polynomials. The approximating equations within element m are the following: M
M
7
-
7.
M
Note that the state profile is approximated using N +' 1 collocation points because the initial point of each element is also taken into account in the approximation.
4462 Ind. Eng. Chem. Res., Vol. 34, No. 12, 1995
The time derivative of the state profile is easily derived from the approximation and is the following:
Note that, due to the Lagrange condition:
where d,,
is the Kronecker delta. Therefore,
If the above approximations are written for every element and every collocation point and substituted into the differential equation, the result is an algebraic residual equation in terms of the unknown coefficients xmn= x, zmn z , and 8,, = 8. The following equations must also be added to enforce continuity of the state profile a t the beginning of each element:
Additional equations can be added, e&., t o control the error of the approximation, t o bound the element lengths, etc. (Cuthrell and Biegler, 19871,but in any case the result is that the initial system of differential and algebraic equations is transformed into a system of algebraic equations that depend on the collocation coefficients x, z, and 8, as well as on the element durations a,.
Literature Cited Biegler, L. T. Solution of Dynamic Optimization Problems by Successive Quadratic Programming and Orthogonal Collocation. Comput. Chem. Eng. 1984,8,243. Cuthrell, J. E.; Biegler, L. T. On the Optimization of Differential Algebraic Process Systems. AIChE J . 1987,33, 1257. Dimitriadis, V. D.; Shah, N.; Pantelides, C. C. Model-Based Safety Verification of Discrete/Continuous Chemical Processes. Presented at the AIChE Annual Meeting, San Fransisco, 1994. Geoffrion, A. M. Generalized Benders Decomposition. J . Optim. Theory Appl. 1972,10,237. Grossmann, I. E.; Morari, M. Operability, Resiliency and Flexibility-Process Design Objectives for a Changing World. In Proceedings Second International Conference Foundations of Computer Aided Process Design, Snowmass, CO, 1983. Grossmann, I. E.; Floudas, C. A. Active Constraint Strategy for Flexibility Analysis in Chemical Processes. Comput. Chem. Eng. 1987,11, 675. Grossmann, I. E.; Straub, D. A. Recent Developments in the Evaluation and Optimization of Flexible Chemical Processes. In Proceedings Computer Oriented Process Engineering, Barcelona, 1991.
Halemane, K. P.; Grossman, I. E. Optimal Process Design under Uncertainty. M C h E J . 2983,29,425. Holt, B. R.; Morari, M. Design of Resilient Processing Plants V-The Effect of Deadtime on Dynamic Resilience. Chem. Eng. Sci. 1985a,40,1229. Holt, B. R.; Morari, M. Design of Resilient Processing Plants VI-The Effect of Right Half Plane Zeros on Dynamic Resilience. Chem. Eng. Sci. 1985b,40,59. Kotjabasakis, E.; Linnhoff, B. Better System Design Reduces Heat Exchanger Fouling Costs. Oil Gas J . 1987,85. Logsdon, J. S.; Biegler, L. T. Decomposition Strategies for LargeScale Dynamic Optimization Problems. Chem. Eng. Sci. 1992, 47,851. Mohideen, M. J.;Perkins, J. D.; Pistikopoulos, E. N. Incorporation of Flexibility and Controllability in Optimal Design of Distillation Columns. To be presented at the IChemE Annual Meeting, Edinburgh, 1995. Morari, M. Design of Resilient Processing Plants 111-A General Framework for the Assessment of Dynamic Resilience. Chern. Eng. Sci. 1983,38,1881. Papalexandri, K. P.; Pistikopoulos, E. N. An MINLP Retrofit Approach for Improving the Flexibility of Heat Exchanger Networks. Ann. Oper. Res. 1993,42,119. Pistikopoulos, E. N.; Mazzuchi, T. A. A Novel Flexibility Analysis Approach for Processes with Stochastic Parameters. Comp. Chem. Eng. 1990,14,991. Skogestad, S.; Morari, M. Design of Resilient Processing Plants E-The Effect of Model Uncertainty on Dynamic Resilience. Chem. Eng. Sci. 1987,42,1765. Soroush, M.; Kravaris, C. Optimal Design and Operation of Batch Reactors. 1. Theoretical Framework. Ind. Eng. Chem. Res. 1993a,32,866. Soroush, M.; Kravaris, C. Optimal Design and Operation of Batch Reactors. 2. A Case Study. Ind. Eng. Chem. Res. 199313,32, 882. Swaney, R. E.; Grossmann, I. E. An Index for Operational Flexibility in Chemical Process Design. Part I: Formulation and Theory. M C h E J . 1985a,31,621. Swaney, R. E.; Grossmann, I. E. An Index for Operational Flexibility in Chemical Process Design. Part 11: Computational Algorithms. A K h E J . 1985b,31,631. Vassiliadis, V. S.;Sargent, R. W. H.; Pantelides, C. C. Solution of a Class of Multistage Dynamic Optimization Problems. 1. Problems without Path Constraints. Ind. Eng. Chem. Res. 1994a,33,2111. Vassiliadis, V. S.; Sargent, R. W. H.; Pantelides, C. C. Solution of a Class of Multistage Dynamic Optimization Problems. 2. Problems with Path Constraints. Ind. Eng. Chem. Res. 1994b, 33,2123. Walsh, S.; Perkins, J. D. Application of Integrated Process and Control System Design to Waste Water Neutralization. Comput. Chem. Eng. 1994,18,S183. White, V.; Perkins, J. D.; Espie, D. M. Switchability Analysis. In Proceedings IFAC Workshop on Integration of Process Design and Control, Baltimore, 1994. Received for review December 9,1994 Revised manuscript received May 30, 1995 Accepted August 11, 1995@
IE940730D Abstract published in Advance ACS Abstracts, November 15, 1995. @
