Adaptive Switching Structure Detection for the Solution of Dynamic

Ind. Eng. Chem. Res. 2006, 45, 8083-8094

8083

Adaptive Switching Structure Detection for the Solution of Dynamic Optimization Problems Martin Schlegel† and Wolfgang Marquardt* Lehrstuhl fu¨r Prozesstechnik, RWTH Aachen UniVersity, D-52056 Aachen, Germany

Recent research in the field of dynamic optimization is striving for increased robustness and efficiency of the numerical solution methods. In addition, the physical interpretation of the optimal solution, with respect to the control switching structure, has been gaining increased interest, because it provides valuable insight that is useful for applying the optimization results on the process. In a previous publication, we have addressed this issue by proposing an automatic switching structure detection method. This paper presents an extension of this idea by combining the structure detection approach with a control profile adaptation technique. The resulting adaptive structure detection approach enables a faster numerical solution, which is more likely to contain the correct switching structure. This further enhances the numerical solution of dynamic optimization problems. The concept is illustrated by means of a small example problem and a large-scale industrial case study. 1. Introduction Dynamic optimization provides a rigorous mathematical framework for the optimization of batch process operation or transient phases in continuous process operation. However, despite its conceptual strength and broad applicability, dynamic optimization currently is a methodology still far from widespread industrial acceptance. It remains a challenge to obtain a highquality solution sufficiently fast, especially for problems that involve large-scale process models. Even the results for small problems, which can be obtained quite quickly with today’s methods and computers, do not always show a satisfactory solution quality. Often, an interpretation of the numerical results is required. In most cases, there is no means to solve a dynamic optimization problem analytically. Rather, approximative numerical solution techniques based on a discretization of the process variables are applied. The accuracy of a numerical solution and the solution speed are mainly governed by the type of discretization used. Primarily, this applies to the control variables. The true solution of a dynamic optimization problem consists of distinct intervals, the so-called “arcs”.1 The control variables are continuous and differentiable on each arc but can jump from one arc to the next at the “switching points”. For example, the optimal feed rate of a semi-batch reactor might stay at a given upper limit for a certain time period, until it suddenly jumps to a lower level to avoid the reactor temperature exceeding its upper bound. Later during the operation, the reactor holdup might hit a specified limit, so that the feed has to be switched off. Thus, an optimal control profile might show a rather nonuniform behavior over the optimization horizon, which is referred to as the “switching structure”. This fact still poses a challenge to all numerical solution methods. The quality of the solution is dependent on the chosen discretization for the control variables and can be insufficient if the discretization does not properly reflect the switching structure, which is typically not known before the problem has been solved. * To whom all correspondence should be addressed. Tel.: +49241-80-94668. Fax: +49-241-80-92326. E-mail: marquardt@ lpt.rwth-aachen.de. † Current address: BASF Aktiengesellschaft, 67056 Ludwigshafen, Germany.

However, this structure is actually one key result when it comes to an implementation of the optimal solution on a real process. Certainly, the size and complexity of the dynamic process model has a significant impact on the computational performance of any optimization algorithm. “Large-scale” dynamic optimization involves process models that contain many (i.e., several hundred or thousand) equations and variables.2,3 However, the computational effort required for the solution is not only governed by the model size, but also by the problem complexity, i.e., the number of control variables, path constraints, and the complexity of the switching structure. Therefore, the focus of recent research is on advanced algorithms for further extending the range of possible applications of dynamic optimization, in terms of model size and problem complexity. In this context, we proposed an adaptive strategy for the automated refinement of the control profile discretization in a previous paper.4 Starting from an initial coarse grid, the parametrization is successively refined based on a wavelet-based analysis of the optimal solution obtained in the preceding optimization step. This leads to nonuniform grids on which the optimization can be evaluated with less computational effort than that on a uniform mesh of comparable accuracy. Alternatively, it is possible to extract the switching structure from a numerical solution automatically, as we have shown in a recent paper.5 This structural information yields additional insight into the problem, which is particularly valuable when the results are actually applied on a real process. Furthermore, it can be used for a very efficient parametrization of the control profiles, because the decision variables used for the parametrization are specifically selected for a best fit to the switching structure. An exploitation of the switching structure has been attempted before in various works (Winderl and Bu¨skens,6 Szymkat and Korytowski,7 Kaya and Noakes8). However, in all of the previous work, the parametrization approach has been restricted to a special class of dynamic optimization problems, such as time-optimal control, resulting in bang-bang solutions. The present work is covering a very large class of nonlinear dynamic optimization problems instead. These two approaches can be nicely combined to a unified solution framework that comprises the advantages of both concepts. This “adaptive switching structure detection” framework is the main contribution of this paper. It will be introduced

10.1021/ie060496e CCC: $33.50 © 2006 American Chemical Society Published on Web 10/31/2006

8084

Ind. Eng. Chem. Res., Vol. 45, No. 24, 2006

and illustrated by means of a small-scale but still challenging example. To show the applicability to a realistic, large-scale problem, we further present a case study that involves an industrial polymerization process. The paper is organized as follows. In section 2, we will present some fundamentals on the problem formulation and the solution method used in this paper. Section 3 reviews the concepts of multiscale control profile adaptation and automatic structure detection. The adaptive switching structure detection is the topic of section 4. Section 5 describes the experiences made with applying the suggested method to the industrial case study problem. Finally, in section 6, we conclude the paper with a summary and future perspectives. 2. Preliminaries 2.1. Problem Formulation. The general formulation of a dynamic optimization (DO) problem used in this paper reads as follows:

min Φ ) ΦM(y(tf))

(DO)

u(t),tf

s.t. Bx3 ) f(x(t), z(t), u(t), t)

t∈I

(1a)

0 ) g(x(t), z(t), u(t), t)

t∈I

(1b)

0 ) x(t0) - x0 0 g h(y(t), u(t)) 0 g e(y(tf))

(1c) t∈I

(1d) (1e)

x(t) ∈ Rnx and z(t) ∈ Rnz are the vectors of the differential and algebraic state variables, respectively. The entire state vector is denoted as y ) [xT, zT]T. The dynamic process model is given by eqs 1a and 1b in semi-explicit form, with initial conditions given in eq 1c. The problem is defined on the time horizon I ) [t0, tf]. The time-dependent control variables u(t) ∈ Rnu and possibly the final time tf are decision variables for the optimization. The goal of the optimization is to find the optimal set of decision variables to minimize the objective function Φ. Without a loss of generality, we assume a Mayer-type objective function, i.e., a terminal cost criterion, ΦM(y(tf)). The search space for finding the optimum is restricted by constraints, which describe, for example, safety or other operational requirements to be fulfilled during the operation of a plant. We consider path constraints h(‚‚‚) and endpoint constraints e(‚‚‚). 2.2. Necessary Conditions of Optimality and Switching Structure. Optimal control theory allows a reformulation of a problem such as that described by the dynamic optimization equation (eq DO) as that of minimizing a Hamiltonian function, using Pontryagin’s Minimum Principle.9 For the sake of brevity, a detailed derivation of this reformulation shall not be given here; instead, we refer to, e.g., the work of Bryson and Ho.1 Pontryagin’s principle forms the basis for a class of solution methods for dynamic optimization problems, the so-called “indirect methods” (cf. section 2.3). Furthermore, it allows the derivation of the “necessary conditions of optimality” (NCO). These provide insight into possible solution structures of the optimal control profiles u(t), as described by Srinivasan et al.10 The behavior of a particular optimal control profile ui(t) on an arc k can be classified into four different possible types, namely, (1) ui,k(t) ) ui,min: the control is at its lower bound, (2) ui,k(t) ) ui,max: the control is at its upper bound,

(3) ui,k(t) ) ui,k,path(t): the control is determined by an active state path constraint, (4) ui,k(t) ) ui,k,sens(t): the control is not determined by any active constraint; it is sensitivity-seeking. In these possible types, ui,k denotes the profile of control variable ui(t) on the arc k. The optimal solution of a given dynamic optimization problem always has a distinct sequence of arcs, the “switching structure”. The classification into the different types of arcs is useful for the automated detection of the switching structure, as presented by Schlegel and Marquardt.5 Note that this method requires neither the derivation of a Hamiltonian function nor the formulation and solution of an adjoint system. Rather, only the classification of an optimal control trajectory into the aforementioned four possible arc categories is used. Strictly speaking, the application of optimal control theory to derive NCO and a control switching structure is restricted to problems where the dynamic constraints are present in form of an ordinary differential equation (ODE) model. It turns out that the derivation of an equivalent theory for DAE-constrained optimization problems has been considered by researchers just recently, at least for nonlinear DAEs. NCO explicitly derived for DAE-constrained problems have been derived by Roubı´cˇek and Vala´sˇek11 and de Pinho and Vinter.12 However, the derivations in both references only treat bounds on the control variables; neither state path constraints nor endpoint constraints are considered. Apparently, research is still required to derive general NCO for problems of the form (DO), in particular, for higher-index systems. For our work, we assume that the classification into four possible types of arcs also exists for the DAE-constrained problems under consideration. 2.3. Numerical Solution. Solution strategies for dynamic optimization problems of the form described above can be distinguished into three basic classes: dynamic programming methods, indirect methods, and direct methods, with specific advantages and disadvantages. A comparative overview has been given, for example, by Binder et al.13 and by Srinivasan et al.10 The reader is referred to these surveys and the literature cited therein for details. For large-scale dynamic optimization problems, such as those which arise in chemical engineering, direct methods have proven to be most successful.13 Direct methods transform the original problem into a nonlinear programming problem (NLP) via a suitable discretization of the state and control variables and the constraints. A key advantage is that they do not require an explicit treatment of the necessary conditions of optimality (and, hence, are straightforwardly applicable to problems constrained by large-scale DAE models). In this work, we make use of the direct sequential method for dynamic optimization (also known as single shooting method or control vector parametrization). Here, only the control variables u(t) are discretized and form the set of decision variables for the optimization problem. The dynamic process model described by eqs 1a and 1b is solved numerically in each iteration step of the NLP solver via a suitable numerical integration method. Thus, the discretization of the state variables y(t) is performed implicitly by the integration and is typically error-controlled. Because process models that stem from practical applications are often quite stiff, error control is an advantage to ensure good approximation of all state profiles. The sequential approach is particularly well-suited for the solution of large-scale problems, where the control variables can be parametrized by few decision variables. A reasonable parametrization of the control profiles is crucial for the performance of the sequential method, because the computa-

Ind. Eng. Chem. Res., Vol. 45, No. 24, 2006 8085

Figure 1. Function and corresponding time-scale plot.

tional cost increases with the number of decision variables, whereas the chances for obtaining a high solution accuracy of the optimal control profiles are generally higher for a fine discretization. This tradeoff is a major motivation for the adaptive methods described in sections 3 and 4. For the parametrization of the control profiles ui(t), often piecewise polynomial approximations are applied. Using a B-spline representation14 for this purpose, we can write the discretized control variables as nui

ui(t) )

uˆ i,j φ(m) ∑ j j)1

(2)

with nui denoting the number of parametrization functions for the control variable ui. Depending on the choice of the order m of the B-spline function φ(m), different approximation orders can be realized. Our solution framework uses piecewise constant (m ) 1) and piecewise linear (m ) 2) parametrizations. ∆ui denotes the set of discretization time points for each control variable ui(t). After discretization of the control variables, the dynamic optimization problem (DO) is reformulated as an NLP, with the discretization parameters uˆ ) [uˆ T1 , ..., uˆ nTu]T as decision variables. The path constraints are evaluated pointwise at all time points tj contained in a set ∆y, for example, on the unified nu mesh of all control variables ∆y ) ∪i)1 ∆ui. A typical method for the solution of NLP problems that are encountered in direct dynamic optimization methods is sequential quadratic programming (SQP) (e.g., Nocedal and Wright15). 2.4. Wavelets. One advantage of using a B-spline basis for control discretization is its relation to a certain class of wavelet functions. Wavelets are well-suited for a multiscale or multiresolution analysis of functions. A comprehensive description of wavelets and the theory of multiresolution methods is beyond the scope of this paper. The reader is referred to the extensive literature in this field (e.g., Dahmen16). As described by Schlegel et al.,4 such a wavelet analysis of control profiles can be used for the derivation of an automated adaptive refinement strategy for the control discretization. This concept will also be used in the adaptive switching structure detection framework to be described in section 4. The key idea of a multiresolution method is to develop a hierarchical representation of a function with a collection of coefficients, each of which provides some local information about the frequency of the function. In our context, this concept is applied to the control variables. By applying a wavelet transformation to the discretized vector uˆ , we obtain an equivalent representation encoded in the wavelet vector d. The magnitude of the wavelet coefficients dj,k is a measure for the local frequency content of the function. A graphical representation can be obtained using the time-scale plot, where the absolute value of every dj,k can be visualized by a particular color of a patch. Figure 1 shows a function and the time-scale plot of the

corresponding wavelet decomposition. Some of the patches are marked by their corresponding wavelet coefficients. Such a plot gives a useful visualization of the level of local detail contained in the analyzed function. A special feature of the wavelet transformation, the norm equivalence, is very useful for the derivation of adaptation techniques. It means that discarding small coefficients dj,k will cause only small changes in approximate representations of the original function (or the control profile in our case). Furthermore, functions that exhibit strong variations only locally, can be approximated very accurately by linear combinations of only relatively few wavelet basis functions corresponding to the significant coefficients dj,k. Transferred back into the time domain, this means that a control profile ui(t) will be discretized only at those points in I, which are necessary for a prespecified accuracy. In other words, the discretization mesh ∆ui is not necessarily equidistant. Furthermore, it can be different for every control variable ui(t). 3. Methods for Adaptation of the Control Profile Discretization As discussed previously, optimal control profiles exhibit a certain switching structure. On the other hand, direct numerical solution methods such as the sequential approach work with discretizations of those profiles. Hence, adaptation is the key concept to achieve accurate, yet computationally efficient, solutions. The problem of choosing optimal discretization grids for direct dynamic optimization problems has been addressed by several authors. In the context of direct simultaneous methods, various approaches have been suggested by Cuthrell and Biegler,17 von Stryk,18 Schwartz,19 Tanartkit and Biegler,20,21 Betts and Huffmann,22 and Binder.23 Similar concepts applied in the context of a sequential solution approach have been used by Vassiliadis et al.,24 Waldraff et al.,25 and Canto et al.26 In the past, we have developed two different methods for control profile adaptation for the sequential approach. In the following discussion, both will be summarized and illustrated by means of an example problem. 3.1. Multiscale Control Profile Adaptation. The adaptation method proposed in our previous paper4 uses the aforementioned properties of wavelets. To generate the desired problem-adapted control variable mesh, the optimization and a subsequent mesh refinement procedure are performed in a repetitive procedure, starting from initially coarsely discretized control profiles uˆ 0i (∆u0i). For this purpose, in each refinement step l, the is inspected by an a posteriori analysis. previous solution uˆ l-1 i Based on this analysis, a new discretization for the next optimization run is generated. The problem is then solved again on the improved discretization grid ∆ul i, where the interpolated old solution uˆ l-1 is used as an initial guess. i l The a posteriori analysis of the optimal solution uˆ /l i (∆ui) observed in refinement step l is implemented by means of the wavelet concept explained in section 2.4. The result of the analysis is a new discretization mesh ∆ul+1 , which is better i adapted to the solution structure and used as the starting point for the subsequent optimization run. The analysis involves two steps: a grid-point elimination and a grid-point insertion step. The idea of grid-point elimination is to remove unnecessary grid points in the representation of uˆ /l i , the optimal solution at refinement step l. This can be easily accomplished in the wavelet space using the norm equivalence (see section 2.4), i.e., neglecting sufficiently small wavelet coefficients. Analogously,

8086

Ind. Eng. Chem. Res., Vol. 45, No. 24, 2006 Table 1. Bioreactor Example: Computational Statistics for Adaptive Refinement Approach

Figure 2. Bioreactor example: Solution profiles and wavelet coefficients of the feed rate F for l ) 0 and l ) 7.

large wavelet coefficients are significant for an accurate representation of the profile and, therefore, are candidates for a refinement by grid-point insertion. The grid-point insertion algorithm identifies those wavelet coefficients that have the largest absolute values. Back transformation into the time domain then tells us where to insert additional grid points to the discretization. For the sake of brevity, we refer the reader to a previous paper4 for more details and, instead, introduce an example to illustrate the concept. 3.2. Example: Bioreactor (Part 1). This case study involves the optimal production of a secreted protein in a fed-batch reactor. It was originally formulated by Park and Ramirez27 and has been considered by several researchers. Here, the modified version as presented by Canto et al.28 is used. The dynamic model contains five differential and five algebraic equations and can be found in Appendix A. The objective is the maximization of the amount of secreted protein, satisfying bounds on the feed rate only. Therefore, the dynamic optimization problem is free of state path and endpoint constraints and reads as follows:

minΦ(tf) F(t)

s.t. model equations 3-12, Fmin e F(t) e Fmax with Fmin ) 0.0 L/h, Fmax ) 2.0 L/h, and tf ) 15 h. The example is quite challenging; numerical studies show that it requires an above-average number of NLP iterations. The optimal control trajectory has an interesting structure, with the single-control variable switching several times between the bounds and sensitivity-seeking arcs. If we solve the problem with the control profile adaptation approach, we obtain the optimal solution shown in Figure 2. (For some details about the software implementation, see Appendix B.) For brevity, only the refinement iterations l ) 0 and l ) 7 are shown. Note that a piecewise linear control vector parametrization has been used here. The resulting profile shows that, indeed, parametrization functions are inserted nonuniformly at those positions, where they are required for a better approximation of the optimal control profile. The objective function value successively improves over the iterations and converges to an end value, as

l

CPU time

accumulated CPU time

number of degrees of freedom

Φl

0 1 2 3 4 5 6 7

1.45 3.31 8.98 21.42 26.52 54.80 43.92 61.20

1.45 4.77 13.75 35.17 61.69 116.49 160.41 221.61

9 15 21 33 46 56 71 84

-32.01138 -32.44776 -32.59591 -32.67345 -32.68496 -32.69464 -32.69505 -32.69518

shown in Table 1. This table also reveals the computational statistics for this problem. For each refinement iteration l, the required computer processing time (CPU) time, in seconds (on a 2.2 GHz AMD Athlon PC), and the accumulated CPU time over all iterations, the number of decision variables (degrees of freedom) and the optimal objective function value are given. Other examples that illustrate this adaptation concept can be found elsewhere.4 3.3. Automatic Detection of the Switching Structure. Despite its advantages, the adaptive refinement method presented in the previous section exhibits potential for further improvement. For the aforementioned example, we observe an accumulation of grid points around discontinuities in the profile. Those are a result of wavelet trees, which are built by the refinement and are inherent to the method. However, for a representation of the true solution, most of these grid points would actually not be required. Obviously, this is related to the control switching structure, which has not been taken into account (such as in all direct numerical solution methods). However, there is a solution for this problem: as described by Schlegel and Marquardt,5 it is possible to automatically detect the control switching structure. This method involves three main steps: (i) solution of the problem (DO), (ii) detection of the arcs in the optimal control profiles, and (iii) reparametrization and solution as a multistage problem (MSP). The purpose of the reparametrization in step (iii) is to utilize the results of step (ii) to obtain a numerical solution with fewer decision variables but higher accuracy than that obtained in step (i). Each arc determined in step (ii) is treated as a separate stage (hence the term multistage problem), where the parametrization order (constant or linear) is chosen according to the type of arc: On arcs with active control variable constraints (umax, umin), those can be excluded from the set of degrees of freedom. Sensitivity-seeking and path-constrained arcs can be approximated reasonably well by piecewise linear functions, knowing that optimal profiles are always continuous on the arcs. Furthermore, the multistage approach includes the lengths of the stages as decision variables in the optimization. Therefore, the solution of this problem yields precise values for the switching points present in the control structure at the optimal point. For technical details about the particular steps of the algorithm, we refer to Schlegel and Marquardt.5 The final output of the solution structure detection method is the number of arcs L, as well as the starting times of the individual arcs τ ∈ RL. Furthermore, we obtain the information about the types of the arcs, which will be denoted by the Boolean vector A of size L, where each Ai will be one out of the set {umin, umax, upath, usens}. 3.4. Example: Bioreactor (Part 2). For illustration of this concept, we revisit the bioreactor example. If we solve the problem following the previously described approach, we obtain the optimal control profile depicted in Figure 3. The dashed line shows the optimal control profile reached after step (i). Here,


Figure 3. Bioreactor example: optimal control profile obtained by the structure detection approach (the dashed line represents the problem (DO), whereas the solid line represents the multistage problem (MSP)).

a discretization of 32 piecewise constant B-splines has been used. The analysis of this profile in step (ii) yields a switching structure with K ) 4 arcs and A ) {usens, umax, umin, usens}. The solution of the final step (iii), which solves the reparametrized multistage problem (MSP), is shown by a solid line. Note that the switching points τ between the arcs, now being decision variables, have been moved by the optimizer. Furthermore, the approximation order is now adapted to the type of the particular arcs: piecewise linear profiles (with three degrees of freedom) lead to a relatively good representation of the sensitivity-seeking arcs. Overall, the reparametrized problem contains only 10 decision variables and can be solved quickly. The computational statistics are shown in Table 2. 4. Adaptive Switching Structure Detection A key assumption for the structure detection method presented in Section 3.3 to work properly is that all relevant arcs of the control switching structure are already present in the singlestage solution. If this is not the case, then certain arcs might be missed. On the other hand, even when, for example, sensitivityseeking or path-constrained arcs are detected correctly, the initial reparametrization might be too coarse for a satisfactory solution quality. The concepts of adaptation (section 3.1) and reparametrization (section 3.3) provide two perspectives on the same problem, i.e., obtaining the proper discretization that best reflects the true optimal solution. Hence, it is self-evident to combine the concepts in some way to exploit the advantages of both. This is the topic of this section and the main contribution of this paper. 4.1. Adaptation of the Single-Stage Solution. The construction of the structure detection method is built on the information about the status of the control and state path constraints. Thereby, the qualitative information on the sequence of the arcs is most important. As long as the initial guess is sufficiently close to the optimum, the solution of the multistage reformulation will be able to move the switching points accordingly. Note that this also requires a good initialization of the switching times; see Sager29 for an example on multiple local minima in a multistage formulation. For illustration, let us again at the bioreactor example. If we apply the adaptive refinement method, here for a piecewise constant discretization starting with eight intervals, we obtain the solutions shown in Figure 4 for the first three iterations. In addition, the labeling of the arcs resulting from the structure detection are shown. It is not surprising that the detected structure Al changes over the refinement iterations. The solution in Figure 4a contains only one (sensitivity-seeking) arc over the entire time horizon. The constraints on the control variable are not active at any

point on the time horizon in the single-stage solution. Hence, it is impossible to catch the umax and umin arcs from this solution. The solution structures in Figure 4b and c are the same and consist of three arcs. Two conclusions can be drawn from this example. On one hand, from a structural perspective, a fine resolution in the single-stage solution does not provide any significant additional benefit, compared to a coarser solution as long as the switching structure remains the same and the switching points are sufficiently close to the optimum. On the other hand, a coarse single-stage solution is more likely to miss certain arcs in the solution structure. Reconsidering the adaptation technique presented in section 3.1, we can use its refinement properties in this context. If we intend to switch to the multistage solution anyway (step (3)), it is sufficient to stop at some early point of the refinement. Clearly, the question arises, what should be the stopping criterion in this case? It is most likely not possible to derive a theoretically verified criterion for this purpose. Nevertheless, a straightforward heuristic can be developed. The computational effort of the structure detection algorithm is negligible; hence, it can be executed after each refinement solution of the multiscale-based adaptation framework. It then is easy to monitor any qualitative change in the solution structure, which can be measured in terms of the number Kl and sequence Al of the arcs, but not in terms of concrete values for the switching times, etc. A simple heuristic is then to stop the refinement as soon as two subsequent iterations have the same qualitative switching information. Typically, the qualitative structure as well as the quantitative data obtained by this approach are well-suited for the initialization of the variables in the NLP arising from the multistage formulation. 4.2. Adaptation of the Multistage Solution. Besides adapting the single-stage profiles, it is also possible to use adaptive refinement on the multistage level. The primary goal of the reparametrization approach is to achieve a reasonable solution accuracy with a minimum number of decision variables. However, it might be also preferable to find solutions of higher accuracy with more decision variables at a higher computational expense. As we already know, adaptation is a well-suited concept to trade off accuracy and computational effort. If a solution is sought, which properly captures the switching points, but is coarsely discretized on path-constrained or sensitivityseeking arcs, then a single multistage solution is sufficient. However, if even higher accuracy on those inactive control arcs is desired, the adaptation can be performed separately on these arcs in the same way as we applied it on the entire time horizon in section 3.1. This is actually attractive for the following reason: we saw that a drawback of the adaptation method is the accumulation of unnecessary mesh points around the switching points, which is due to the formation of the wavelet trees. Now, in the multistage setting, this problem is resolved, because the adaptation is performed separately on the various arcs. Assuming a properly detected switching structure, there will be no switching time on a particular arc, rather the switching times are between the arcs and are taken care of by the multistage formulation. Hence, an unnecessary mesh point accumulation will not occur. Instead, the adaptation will be even more effective than applied in a straightforward manner to a single-stage problem. This idea also has a relation to an issue discussed previously, namely the possibility that certain arcs in the solution structure might be missed because of a too-coarse resolution of the singlestage control profiles. The adaptation of the single-stage solution

8088


Figure 4. Bioreactor example: adapted solution (piecewise constant) for (a) l ) 0, (b) l ) 1, and (c) l ) 2. Table 2. Bioreactor Example: Computational Statistics for Structure Detection Approach

dynamic optimization, DO multistage problem, MSP

CPU time



Φ*

7.97 1.98

7.97 9.95

32 10

-32.49540 -32.41390

can reduce this risk. However, an additional, complementary, measure is the adaptation of the multistage solution. The idea is to execute the structure detection procedure independently per stage after each refinement. This results in two levels of arcs. On the upper level, there are the arcs resulting from the analysis of the single-stage solution, which have been converted to the K stages of the multistage problem. On the lower level, each stage k is again analyzed for possible “sub-arcs”. This result allows an analysis of the multistage solution for a possible change of the switching structure. Two cases can be distinguished: (1) By construction, each stage k contains exactly one subarc. If, however, the solution of the multistage problem leads to a number of sub-arcs that is larger than one, this indicates a possible change of the switching structure. (2) Even if the number of sub-arcs per stage remains 1 after the solution, it still might be the case that the type of the subarc (and, thereby, the entire stage) has changed. The first case is most likely to occur during successive refinement of the multistage solution. For example, it could happen that on an arc, which has been detected and parametrized as sensitivity-seeking, during the refinement suddenly the control variable hits its lower or upper bound. An additional umin or umax arc, respectively, then is detected. Subsequently, the problem can be again split up into further stages at the corresponding activation and/or deactivation points. An example for the second case is that an arc switches to the type of a neighboring arc. 4.3. Adaptive Switching Structure Detection. A unified adaptation framework can now be established by linking adaptation, structure detection, and reparametrization along the lines developed in the previous sections. The technical procedure of this concept is summarized in form of Algorithm 1. Here, the user-specified lmax and κmax provide upper limits on the two refinement iterations. The refinement stops according to |(Φκ - Φκ-1)/Φκ| > Φ, when no further significant improvement of the objective function is possible (cf. Schlegel et al.4). This unified framework adds an additional level of adaptation. The pure single-stage adaptation adjusts the resolution of the control variable profiles, and a pure multistage solution refines the position of the switching points in a given sequence of arcs. The combined approach refines both the control resolution and the sequence of the detected arcs. It allows one to detect the switching structure and obtain an accurate control profile resolution, even for complicated cases. Details about the software implementation can be found in Appendix B. 4.4. Example: Bioreactor (Part 3). Figure 2 shows that the optimal solution profile of the feed rate F reveals a complex

Algorithm 1. Adaptive Switching Structure Detection Choose initial control profiles and grids: u0i (t), ∆u0i, i ) 1, ..., nu Choose maximum number of refinement iterations: lmax, κmax Step 1: Single-Stage Adaptation for l ) 0, 1, ... do /l Solve discretized optimization problem (NLPl) on ∆ul i f Φl; uˆ /l i , tf Analyze solution structure f Al, Kl if (Al * Al-1 (for l > 0) and l < lmax) then Grid point elimination, insertion f refined grids: ∆uil+1 Interpolate profiles: ul+1 :) interp(uli, ∆ul+1 ) i i else Go to Step 2 end if end for f optimal single-stage solution on adaptive grid: u/i , Φ*, A*, Kl Step 2: Multistage Adaptation K0 ) Kl Set up multistage problem (MSP0) for κ ) 0, 1, ..., do κ solve optimization problem (MSPκ) on ∆uκi,k f Φκ; uˆ κi,k, t/κ k , k ) 1, ..., K Analyze solution structure f Aκ, Kκ if (|(Φκ - Φκ-1)/Φκ| > Φ (for κ > 0) and κ < κmax) then Grid point elimination, insertion f refined grids: ∆uκ+1 i,k κ κ+1 κ Interpolate profiles: uκ+1 i,k :) interp(ui,k, ∆ui,k ), k ) 1, ..., K else Exit end if if Aκ * Aκ-1 (for κ > 0) then Set up multistage problem (MSPκ+1) end if end for f optimal multistage solution on adaptive grid: u/i ) uκi , Φ* ) Φκ, A* ) Aκ

switching structure, with some of the arcs being rather short in duration. Thus, this problem is an excellent example to illustrate the adaptive structure detection concept. The procedure is started with an adaptation of the single-stage solution departing from a discretization with eight equidistant intervals. As soon as a consistency in the switching structure is detected, it is switched to a multistage solution with adaptation. This point is reached at l ) 2, as shown in Figure 4. Therefore, Algorithm 1 switches to the multistage solution, according to the heuristic. The solution of the first, coarsely parametrized multistage problem (κ ) 0) leads to the optimal control profile shown in Figure 5. Again, a successive refinement of this solution is performed, with an inspection for a change in the switching structure after each iteration. This occurs first in iteration κ ) 3 (cf. Figure 6), because, here, an umax arc appears at the end of the first stage which originated from an arc initially proposed to be of the sensitivity-seeking type. The time-scale plot in Figure 6 also


Figure 5. Bioreactor example: multistage solution for κ ) 0.



shows that the refinement is now performed independently on each stage, but not on the entire time horizon. The solution, which now contains four arcs, is further refined, until, in iteration κ ) 5, a second structure change is detected. A small umax arc appears at the end of the horizon. The procedure stops at κ ) 6. The final solution is shown in Figure 7. This study shows that, indeed, the stopping criterion for the single-stage adaptation might lead to incorrect results, in the sense that still certain arcs may be missed. Nevertheless, the multistage adaptation is able to detect the two missed arcs, which are then included into the solution structure, overall leading to the correct result. The computational statistics are shown in Table 3. In comparison to the adaptive refinement results (Table 1), we see that much less computation time is needed with this approach to obtain a solution of comparable accuracy. Furthermore, we obtain the control switching structure, which is also more precise than the one detected by applying the pure structure detection method (section 3.4). 5. Industrial Polymerization Process In comparison to the previous example, which is a small problem from the literature, we now consider a real industrial polymerization process. For confidentiality reasons, only those details can be presented here which have been published previously in the article by Du¨nnebier et al.30 The flowsheet of this large-scale continuous polymerization process is shown in Figure 8. The exothermic polymerization

occurs in a continuously stirred tank reactor (CSTR) that has been equipped with an evaporative cooling system. The reaction is operated at an open-loop unstable steady-state operating point that corresponds to a medium level of conversion. The polymer is separated from unreacted monomer and solvent, and the recycled monomer and solvent are being fed back via a buffer tank. To achieve the desired properties, downstream processing and blending units transfer the melt into the final product. Crucial process parameters that affect the final product properties include the polymer viscosity, which is related to the molecular weight, and the polymer content at the reactor outlet, which is equal to the conversion in the reactor. The study presented by Du¨nnebier et al.30 considers an optimal change of the reactor load, e.g., from 50% to 100% or vice versa. In contrast, here, an optimal change of the polymer grade is considered. The task is to perform a change of the polymer grade from one specification (grade A) to another one (grade B) in minimum time while enforcing several operational constraints. To achieve this goal, there are three control variables available, which are also shown in Figure 8: the mass flow rates of the monomer (u1) and the catalyst (u2) feed streams, as well as the mass flow rate of the recycle stream (u3). Operational path constraints are enforced on the following quantities: (i) reactor outlet flow rate (y1), (ii) buffer tank holdup (y2), (iii) solid content in the reactor (y3), (iv) polymer molecular weight (y4), (v) reactor temperature (y5), and (vi) condenser outlet flow rate (y6). In addition, there are endpoint constraints on the solid content in the reactor (y3) and the polymer molecular weight (y4) (corresponding to the grade B specification), which are more stringent than those enforced on these variables during the operation. The dynamic process model is of large scale, involving 148 differential and 1949 algebraic variables and equations. A rigorous reactor model including an extensive scheme of reaction kinetics was available from previous design and operation studies. Because of the fact that it has been designed and implemented for a different purpose, the model includes several options and capabilities (such as different catalysts and solvents) that are not needed for the optimization studies presented in this work. Some of the numerical problems reported below are caused by this fact. The model of the reactor and buffer tank is a very rigorous one; however, the separator unit and the condenser have been simplified significantly, based on process insights. The reactor is modeled as a CSTR with heat and mass balance and complex polymerization kinetics for all species involved. The separator has been modeled as an ideal splitter, neglecting its dynamics, compared to the other buffers located in the system. The condenser model is a gray box model, using a steady-state mass balance combined with a second-order linear dynamic model estimated from process data. Because of the model size and a possibly complex control switching structure induced by the large number of control variables and constraints, this problem poses a challenge to any dynamic optimization method and is well-suited to assess the applicability of the methods described in this paper. 5.1. Multiscale Control Profile Adaptation. As a first benchmark, we apply the multiscale control adaptation technique described in Section 3.1 to this problem. The three control variables are initially discretized as piecewise constant functions, choosing eight equidistant intervals for each. The initial profiles take the values corresponding to the steady-state operation for the production of polymer grade A as constant values over the entire time horizon. The solutions after the iteration l ) 5 are shown in Figure 9. (For confidentiality reasons, both axes in

8090


Table 3. Bioreactor Example: Computational Statistics for the Unified Adaptation Approach l

κ

0 1 2 0 1 2 3 4 5 6

CPU time



Φl

1.11 2.72 4.78

1.11 3.83 8.61

8 14 22

-31.76824 -32.30083 -32.49409

1.41 2.38 5.37 5.31 11.53 20.14 6.48

10.02 12.39 17.77 23.08 34.61 54.75 61.23

9 13 19 25 28 39 52

Φk

-32.38856 -32.57897 -32.66648 -32.67582 -32.69151 -32.69479 -32.69514

Figure 8. Industrial polymerization process: process flowsheet.

Figure 9. Industrial polymerization example: adaptive solution of the three control variables for l ) 5.

all solution plots for this example, as well as the objective function value, are scaled.) We can see that the resolution of the control profile discretization has been adapted during the refinement process separately for each control variable. The finest resolution of the catalyst feed rate (u2) shows some chatter during a portion of the second half of the time horizon. This behavior is due to the presence of a sensitivity-seeking arc in that region. Because the sensitivity of the objective function, with respect to the control variable u2, is low in this region, it is difficult for the optimizer to distinguish the influence of control parameters that are adjacent to each other in the discretized profile. For the computational statistics, we refer to Table 4. 5.2. Automatic Detection of the Switching Structure. For comparison, we now apply the switching structure detection technique described in section 3.3 to this problem. The singlestage solution is discretized with 16 equidistant intervals for each control variable. In Figure 10, both the single-stage solution

Table 4. Computational Statistics l

CPU time



Φl

0 1 2 3 4 5a

313.00 608.53 1303.36 1510.70 2727.05 5849.41

313.00 921.53 2224.89 3735.60 6462.64 12312.05

25 31 47 62 78 114

0.857702 0.858290 0.854916 0.854738 0.854271 0.853757

a Iteration l ) 5 did not converge to optimality, because of line-search failure.

(dashed lines) and the resulting multistage solutions (solid lines) are shown, together with the detected switching structure. As it turns out, only two of the six constrained state variables are hitting their bounds in some portions of the time horizon: the first one is the reactor outlet flow rate (y1), and the second one is the buffer tank holdup (y2) (cf. Figure 11). An interpretation of the resulting switching structure can be given as follows: The solution consists of five arcs. The reactor


Figure 10. Industrial polymerization example: solution with structure detection (control variables).

Figure 11. Industrial polymerization example: solution with structure detection (key path-constrained state variables).

outlet flow rate is active during the first four arcs, and the buffer tank holdup is active during the last three arcs. In the first arc, u1 is located at the upper bound and u2 is located at the lower bound, and u3 is used to keep the first path constraint active. When u3 reaches its lower bound, u1 follows a trajectory to keep the first path constraint further active. At the next switching point, the buffer tank holdup (the second path constraint) hits the upper bound and becomes active. Now, all three control variables follow intermediate arcs. Because the sensitivities of u1 and u3, with respect to the two constraints, are of the same order of magnitude, both are used to keep the constraints along their upper bounds. Variable u2 has a very low sensitivity, with respect to the constraints; hence, it could be identified as a sensitivity-seeking control trajectory. This is also consistent with the observations stated in the previous subsection. At the next switching point, u2 jumps to its lower bound, whereas u1 and u2 stay on their path-constrained trajectories. In the final arc, u2 jumps to the maximum, u1 to the minimum, and u3 is used to keep the buffer tank holdup at the active bound, whereas the reactor outlet flow rate has left its bound. The length of the last arc is adjusted such that the two endpoint constraints are met, which can be seen in Figure 12.

5.3. Adaptive Switching Structure Detection. Finally, the adaptive structure detection technique described in Section 4 is applied. The complex switching structure of the three control variable profiles is already known (cf. Figure 9). For brevity, the intermediate refinement results for the iterations l ) 0, ..., 4 are not shown here, but they are reported by Schlegel.31 They show that the switching structure for the catalyst feed rate and the recycle stream change in almost every refinement iteration. This supports the stopping criterion of the single-stage adaptation based on a consistent switching structure, because apparently here, quite some refinement would be necessary to arrive at the correct structure. However, to illustrate the adaptation of the multistage solution for this problem, the single-stage adaptation is not treated here. Instead, we start with a coarse single-stage solution (with eight intervals per control variable) and then switch to the multistage mode with adaptation. Figures 13 and 14 show the optimal profiles for the three control variables after refinement iteration κ ) 0 and κ ) 4. We observe that the initial structure just comprises three stages. Already in iteration κ ) 0, the last arc in the profile for u2

8092


Figure 12. Industrial polymerization example: solution with structure detection (endpoint-constrained state variables).

Figure 13. Industrial polymerization example: solution with adaptive structure detection (κ ) 0).

Figure 14. Industrial polymerization example: solution with adaptive structure detection (κ ) 4).

jumps from the initially detected usens arc to a umax arc. After iteration κ ) 2, a umax arc is detected in the profile for u1. The solution of this problem involved some computational problems. For example, the sensitivity-seeking arc in the catalyst feed rate cannot be solved to a higher accuracy between refinement iterations 2 and 4, although more degrees of freedom are provided by the refinement. Apparently, the NLP solver has convergence problems. This is also the reason for the termination of the refinement in iteration κ ) 4. Nevertheless, these results show that the adaptive structure detection approach is applicable to large-scale problems. The computational statistics in Table 5 show that satisfactory results, in terms of computation time versus objective function value, can be achieved. Compared to the adaptive refinement results reported in Table 4, here, better objective function values are available after significantly shorter CPU times. Another important result is the fact that a coarse control switching structure (with three arcs) is already very close to

the true optimum. In other words, the several small arcs observed for this problem in sections 5.1 and 5.2 do not contribute significantly to an improvement of the objective function value. This has also been confirmed by Kadam et al.,32 who applied the concept of NCO tracking33 to the same polymerization process. 6. Conclusions Various research activities in the field of dynamic optimization in the process engineering community are striving for increased robustness and efficiency of the numerical solution methods. Improving the algorithms to gain further digits of accuracy in shorter computation time is certainly important, but just one side of the coin. As noted, for example, by Srinivasan et al.,10 the physical meaning of the optimization problem and its solution is of equal importance, especially when it comes to applications of the solution, e.g., in a process control context.

Ind. Eng. Chem. Res., Vol. 45, No. 24, 2006 8093 Table 5. Example 4: Computational Statistics l

κ

0 0 1 2 3 4a a

CPU time


number of degrees of freedom, dof

Φl

718.90

718.90

32

0.857700

377.25 455.15 90.34 269.46 816.23

1096.15 1551.30 1641.64 1911.10 2727.33

19 26 44 43 71

Φk 0.854245 0.854609 0.853249 0.853294 0.853210

Iteration κ ) 4 did not converge to optimality.

This aspect is captured by the control switching structure. It might be sufficient or even more important for practical applications to find the correct switching structure rather than just improving the objective function numerically without gaining additional insight. Both of these perspectives have been taken in our previous publications in form of an adaptive refinement concept and an automatic structure detection method, respectively. In this paper, it has been shown that a sophisticated combination of both approaches can further enhance the numerical solution of dynamic optimization problems. Combining adaptation and structure detection in an iterative manner leads to better approximations of the control profiles with fewer degrees of freedom. As an overall result, we can obtain efficient solutions that are more likely to contain the correct switching structure than by applying the two concepts separately. This has been illustrated by means of an example problem used throughout the paper. Furthermore, the adaptive structure detection framework adds another dimension of refinement, because it allows to refine the sequence of arcs from a very simple structure to a possibly rather complex one with various small arcs. This is a new perspective of adjustable solution quality that, so far, is not known from the literature, which can be useful, for example, in on-line applications such as the measurement-based implementation of optimal trajectories via NCO tracking.32,34 A significant portion of the paper has been spent on the description of an industrial example problem. The results presented for this example show that, despite numerical problems, also large-scale problems can be tackled by the proposed method. Current work strives for a consolidation of the algorithms and the underlying software framework to address even more-complex industrial dynamic optimization problems. Acknowledgment The authors gratefully acknowledge the fruitful discussions with G. Du¨nnebier and K.-U. Klatt (Bayer Technology Services GmbH, Leverkusen), D. Bonvin and B. Srinivasan (EPFL Lausanne), and J.V. Kadam (RWTH Aachen University). The software implementation was largely supported by K. Stockmann (RWTH Aachen University).

F (m - S) V

V˙ ) F

(7)

4.75µx 0.12 + µx

(8)

Se-5.01S 0.1 + S

(9)

21.88S (S + 0.4)(S + 62.5)

(10)

Θ(µx) ) fP(S) ) µx(S) )

(6)

Y(S) ) 58.75µx2 + 1.71

(11)

Φ ) -PMV

(12)

The state variables appearing in the aforementioned differential equations are the amount of secreted protein PM and the total protein amount PT (given in units of 1/L), the culture cell density X (given in units of g/L), the culture glucose concentration S (given in units of g/L), and the holdup volume V (given in liters). The glucose feed rate F (expressed in units of L/h) is the control variable, with m being the feed concentration. The variables Θ, fP, µx, and Y denote the secretion rate constant, the protein expression rate, the specific growth rate and the biomass-to-glucose yield, respectively. Appendix B. Implementation All numerical concepts presented in this paper have been implemented into the software tool DyOS.35 The implementation allows the optimization of any model compliant to the CapeOpen ESO interface standard,36 e.g., through the modeling and simulation package gPROMS.37 Adaptation, switching structure detection, reformulation, and solution of the multistage problems are conducted automatically. The example problems have been solved using SNOPT38 as a NLP solver and a modified version of LIMEX39 as a numerical integrator. Different options for solving the problem can be selected by the user. DyOS takes care of the adaptation(s) and also switches automatically from the single-stage to the multistage problem without the need of any further interaction. Literature Cited

Appendix A. Model of the Park-Ramirez Bioreactor

F P˙ M ) Θ(µx)(PT - PM) - PM V

S˙ ) - 7.3µx < (S)X +

(3)

P˙ T ) fP(S)X -

F P V T

(4)

X˙ ) µx(S)X -

F P V X

(5)

(1) Bryson, A.; Ho, Y.-C. Applied Optimal Control; Taylor and Francis: Bristol, PA, 1975. (2) Leineweber, D.; Bauer, I.; Bock, H.; Schlo¨der, J. An efficient multiple shooting based reduced SQP strategy for large-scale dynamic process optimization. Part 1: theoretical aspects. Comput. Chem. Eng. 2003, 27, 157-166. (3) Cervantes, A.; Biegler, L. Large-scale DAE optimization using a simultaneous NLP Formulation AIChE J. 1998, 44, 1038-1050. (4) Schlegel, M.; Stockmann, K.; Binder, T.; Marquardt, W. Dynamic optimization using adaptive control vector parametrization. Comput. Chem. Eng. 2005, 29, 1731-1751.

8094


(5) Schlegel, M.; Marquardt, W. Detection and exploitation of the control switching structure in the solution of dynamic optimization problems J. Process Control 2006, 16, 275-290. (6) Winderl, S.; Bu¨skens, C. Exploiting the special structure for realtime control of optimal control problems with linear controls: Numerical experience. Proc. Appl. Math. Mech. 2002, 1, 484-485. (7) Szymkat, M.; Korytowski, A. Method of monotone structural evolution for control and state constrained optimal control problems. In European Control Conference ECC 2003, Cambridge, U.K., 2003. (8) Kaya, C.; Noakes, J. Computational method for time-optimal switching control J. Opt. Theory Appl. 2003, 117, 69-92. (9) Pontryagin, L.; Boltyanskiy, V.; Gamkrelidze, R.; Mishchenko, Y. The Mathematical Theory of Optimal Processes; Wiley-Interscience: New York, 1962. (10) Srinivasan, B.; Palanki, S.; Bonvin, D. Dynamic optimization of batch processes I. Characterization of the nominal solution. Comput. Chem. Eng. 2003, 27, 1-26. (11) Roubı´cˇek, T.; Vala´sˇek, M. Optimal control of causal differentialalgebraic systems J. Math. Anal. Appl. 2002, 269, 616-641. (12) de Pinho, M.; Vinter, R. Necessary conditions for optimal control problems involving nonlinear differential algebraic equations J. Math. Anal. Appl. 1997, 212, 493-516. (13) Binder, T.; Blank, L.; Bock, H. G.; Bulirsch, R.; Dahmen, W.; Diehl, M.; Kronseder, T.; Marquardt, W.; Schlo¨der, J. P.; von Stryk, O. Introduction to model based optimization of chemical processes on moving horizons. In Online Optimization of Large Scale Systems; Gro¨tschel, M., Krumke, S. O., Rambau, J., Eds.; Springer-Verlag: Berlin, 2001. (14) De Boor, C. A Practical Guide to Splines; Springer-Verlag: New York, 1978. (15) Nocedal, J.; Wright, S. Numerical Optimization; Springer-Verlag: New York, 1999. (16) Dahmen, W. Wavelet and multiscale methods for operator equations. Acta Numerica 1997, 6, 55-228. (17) Cuthrell, J.; Biegler, L. On the optimization of differential-algebraic process systems. AIChE J. 1987, 33, 1257-1270. (18) von Stryk, O. Numerische Lo¨sung optimaler Steuerungsprobleme: Diskretisierung, Parameteroptimierung und Berechnung der adjungierten Variablen; VDI-Fortschrittsbericht, Reihe 8, Nr. 441; VDI-Verlag: Du¨sseldorf, 1995. (19) Schwartz, A. Theory and implementation of numerical methods based on Runge-Kutta integration for solving optimal control problems, Thesis, Department of Electrical Engineering and Computer Sciences, University of California at Berkeley, Berkeley, CA, 1996. (20) Tanartkit, P.; Biegler, L. A nested, simultanous approach for dynamic optimization problemssI. Comput. Chem. Eng. 1996, 20, 735741. (21) Tanartkit, P.; Biegler, L. A nested, simultanous approach for dynamic optimization problemssII: The outer problem. Comput. Chem. Eng. 1997, 21, 1365-1388. (22) Betts, J.; Huffmann, W. Mesh refinement in direct transcription methods for optimal control. Optim. Control Appl. Methods 1998, 19, 1-21. (23) Binder, T. AdaptiVe Multiscale Methods for the Solution of Dynamic Optimization Problems; VDI-Fortschrittsbericht, Reihe 8, Nr. 969; VDIVerlag: Du¨sseldorf, 2002.

(24) Vassiliadis, V.; Sargent, R.; Pantelides, C. Solution of a class of multistage dynamic optimization problems: 1. Problems without path constraints. Ind. Eng. Chem. Res. 1994, 33, 2111-2122. (25) Waldraff, W.; King, R.; Gilles, E. D. Optimal feeding strategies by adaptive mesh selection for fed-batch bioprocesses. Bioprocess Eng. 1997, 17, 221-227. (26) Canto, E.; Banga, J.; Alonso, A.; Vassiliadis, V. Efficient optimal control of bioprocesses using second-order information. Ind. Eng. Chem. Res. 2000, 39, 4287-4295. (27) Park, S.; Ramirez, W. Optimal production of secreted protein in fed-batch reactors AIChE J. 1988, 34, 1550-1558. (28) Canto, E.; Banga, J.; Alonso, A.; Vassiliadis, V. Dynamic optimization of chemical and biochemical processes using restricted second-order information. Comput. Chem. Eng. 2001, 25, 539-546. (29) Sager, S. Numerical methods for mixed-integer optimal control problems; Der andere Verlag: To¨nning, Lu¨beck, Marburg, 2005. (30) Du¨nnebier, G.; van Hessem, D.; Kadam, J.; Klatt, K.-U.; Schlegel, M. Optimization and Control of Polymerization Processes. Chem. Eng. Technol. 2005, 28, 575-580. (31) Schlegel, M. AdaptiVe discretization methods for the efficient solution of dynamic optimization problems; VDI-Fortschrittsbericht, Reihe 3, Nr. 829; VDI-Verlag: Du¨sseldorf, 2005. (32) Kadam, J.; Srinivasan, B.; Bonvin, D.; Marquardt, W. Optimal grade transitions in industrial polymerization processes via NCO tracking. Submitted to AIChE J., 2006. (33) Srinivasan, B.; Palanki, S.; Bonvin, D. Dynamic optimization of batch processes II. Role of measurements in handling uncertainty. Comput. Chem. Eng. 2003, 27, 27-44. (34) Kadam, J.; Schlegel, M.; Srinivasan, B.; Bonvin, D.; Marquardt, W. Dynamic optimization in the presence of uncertainty: From off-line nominal solution to measurement-based implementation. J. Process Control, in press. (35) DyOS User Manual, Release 2.1, Lehrstuhl fu¨r Prozesstechnik, RWTH Aachen University, Aachen, Germany, 2002. (36) Keeping, B.; Pantelides, C. Numerical Solvers Open Interface Specification Draft, Technical Report CO-NUMR-EL-03, CPSE, Imperial College, London, 2000. (37) gPROMS User Guide (Release 2.1.1), Process Systems Enterprise, Ltd., London, 2002. (38) Gill, P.; Murray, W.; Saunders, M. SNOPT: An SQP algorithm for large-scale constrained optimization, Technical Report, Stanford University, Stanford, CA, 1998. (39) Schlegel, M.; Marquardt, W.; Ehrig, R.; Nowak, U. Sensitivity analysis of linearly-implicit differential-algebraic systems by one-step extrapolation. Appl. Numer. Math. 2004, 48, 83-102.

ReceiVed for reView April 20, 2006 ReVised manuscript receiVed August 28, 2006 Accepted August 29, 2006 IE060496E

Adaptive Switching Structure Detection for the Solution of Dynamic

Recommend Documents