Optimal Campaign Continuous Manufacturing - ACS Publications

Article pubs.acs.org/IECR

Optimal Campaign Continuous Manufacturing Ali M. Sahlodin and Paul I. Barton* Process Systems Engineering Laboratory, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, United States ABSTRACT: Campaign continuous manufacturing (CM), characterized by relatively short operational windows such as a few weeks, is being explored as an alternative to traditional batch-wise manufacturing in the pharmaceutical industry. However, optimal operation in campaign CM can be challenging because of the significance of startup and shutdown phases, which can negatively affect on-specification production and plant economy. In this paper, the effectiveness and computational tractability of several known optimization approaches when applied to campaign CM are investigated. Inspired by the turnpike property in optimal control, a new approach is then proposed, which aims to maximize on-specification production explicitly rather than minimizing the startup/shutdown times as commonly adopted in high-volume industries. A main contribution in the new approach is that the resulting optimization formulation is guaranteed to be differentiable, despite the underlying hybrid dynamic system. Thus, it can be solved reliably using gradient-based algorithms. Case studies are presented to demonstrate the effectiveness of the proposed approach.

1. INTRODUCTION Many chemical processes today are designed to operate in a continuous manner, which is characterized by a steady flow of feeds and products into and out of the plant, respectively. Oil refineries and petrochemical facilities are typical examples of continuous manufacturing (CM) processes. On the other hand, many chemical processes are run in a batch-wise manner, in which there is no steady flow through the unit operations, and each unit is often disconnected from the rest of the plant while it is in operation. CM offers several advantages over batch-wise manufacturing, including lower capital and operating costs, improved controllability and product consistency, and lower environmental footprint.1−3 Despite these benefits, some products, such as the majority of specialty chemicals and pharmaceuticals, are manufactured in batch processes.4 The tradition of batch-wise manufacturing in these sectors is attributed to the relatively small production scales5 and high profit margins2 that would demotivate further improvements offered by CM. Moreover, isolated operations in batch-wise manufacturing have the advantage of preventing upstream disturbances from propagating through the plant,6 and enable quality checks at certain points in the process before proceeding to the next steps. This is particularly important for the pharmaceutical industries where stringent quality specifications must be satisfied. Nonetheless, the situation is changing now as increasing competition, new regulations, and advancements in control systems push for a transition from batch-wise to CM.7−9 For the CM alternative to be successful, some operational challenges that are specific to continuous processes must be addressed. Most importantly, product specifications should be met during the operation, instead of only at the final time, as in a batch process. This is typically difficult or impossible to achieve, especially in the startup and shutdown phases, in which the process undergoes major transients. Thus, a common approach to optimizing the operation consists of minimizing startup and shutdown times so that the process is run at the desired steady state for the longest possible period. In doing so, © 2015 American Chemical Society

the product quality specifications are typically ignored to achieve a fast transition. While this approach can be suitable for processes with long operation times (e.g., oil refining), it may not be so for so-called campaign CM, in which the whole operation will take no more than a few weeks, and the transient phases take up a significant part of the period. Therefore, specialized approaches will be needed to ensure optimal operation in such conditions. In this paper, the problem of optimal operation in campaign CM is addressed. Instead of minimizing the transition times, explicit maximization of on-specification production is considered. However, it is found that the resulting dynamic optimization formulation may not exhibit some properties required by standard optimization techniques. Specifically, a hybrid dynamic system may result,10 which is not in general differentiable with respect to optimization variables. Consequently, it may not be possible to solve the optimization problem reliably using gradientbased methods. Supported by the turnpike property in optimal control, and based on the sensitivity theories for hybrid systems, a new formulation is proposed, which preserves differentiability of the problem. The remainder of the paper is organized as follows. In the next section, Conventional Problem Formulations for optimal campaign CM are reviewed, and their computational aspects are discussed. The proposed formulation is presented in section 3. Numerical case studies are given in section 4. Finally, section 5 closes the paper with some concluding remarks.

2. CONVENTIONAL PROBLEM FORMULATIONS In this section, two conventional problem formulations for optimal campaign CM are presented, and computational aspects of each formulation are discussed. Note that the goal in these formulations is produce as much on-specification Received: Revised: Accepted: Published: 11344

April 12, 2015 October 8, 2015 October 13, 2015 October 13, 2015 DOI: 10.1021/acs.iecr.5b01376 Ind. Eng. Chem. Res. 2015, 54, 11344−11359

Article

Industrial & Engineering Chemistry Research product as possible, which is taken as an indication of the profit in this work. 2.1. Maximization of Steady-State Duration. Assuming there exists a steady state that results in on-specification production, one approach to optimal operation is to find a steady state which yields the maximum on-specification production rate and then maximize the time the process is run at that steady state. This is equivalent to minimizing (i) the time it takes for the process to transition from its initial state to the optimal steady state (i.e., startup time) and (ii) the time it takes for the process to transition from the optimal steady state to the final shutdown state (i.e., shutdown time). The optimal steady state can be obtained from the following static optimization problem: max

u ss ∈ U , x ss ∈ X

s.t.

J(x ss, uss)

x|̈ t > tss =

Similarly, it can be shown that all higher-order time derivatives are also zero, and thus, the system has no momentum to depart from the steady state after tss. Now, a third optimization problem is formulated to minimize the shutdown time. Here, it is desired to move the process from the optimal steady state to a final state in the shortest time possible. This goal can be achieved by solving the following optimal control problem: min

u ∈ tss. In this way, the time derivatives of the state variables become zero (0 = f(x*ss , u*ss )), and the state profiles are locked at their optimal steady-state values. It can be shown that the state profiles will not depart from steady state after tss simply by obtaining the second-order time derivatives as 11345

DOI: 10.1021/acs.iecr.5b01376 Ind. Eng. Chem. Res. 2015, 54, 11344−11359

Article

Industrial & Engineering Chemistry Research variables too. Moreover, a similar change of variables as above can be implemented for each time stage. 2.2. Maximization of on-Specification Production. The quality constraints ensuring on-specification production are not enforced in the startup and shutdown phases of the formulation presented in subsection 2.1. As a result, on-specification production is only guaranteed during the steady-state phase, which may be too short depending on the process dynamics and the overall operation time. This motivates other formulations that address on-specification production by explicitly taking the quality specifications into account. To do so, two approaches will be considered, namely, a so-called quality-asobjective (QaO) approach and a quality-as-constraint (QaC) approach. These approaches and their practical implications are discussed in the following subsections. 2.2.1. Quality-as-Objective Approach. In this approach, on-specification production over the entire time horizon is maximized as follows:

∫ u∈< 0

max s.t.

tf

J (̃ x(t , u), u(t )) dt

x(̇ t , u) = f(x(t , u), u(t )),

s(x(t , u), u(t )) ≤ 0,

Figure 1. Sequence of modes visited by the trajectory of a hybrid system.

optimization variables.14 Provided the functions defining each individual mode are differentiable, Galán et al.14 established sufficient conditions for differentiability of the solution of a hybrid dynamic system. For simplicity, consider a constant control u(t) ≔ p̂, ∀t ∈ [0,tf]. According to the sufficient conditions, in order for the solution of the hybrid system to be differentiable with respect to p̂, small variations in p̂ must not change the sequence of modes visited along the solution trajectory. For example, let the sequence of modes given by q(x(t, p̂), p̂) be Mode 1 → Mode 2 → Mode 1, as shown in the top plot of Figure 2. The sufficient condition requires the modes

(P4)

∀ t ∈ (0, tf ]

∀ t ∈ [0, tf ]

(P4.1) (P4.2)

h(x(tf , u), u(tf )) ≤ 0

(P4.3)

x(0, u) = x 0

(P4.4)

with the integral objective function giving the total onspecification product produced over [0, tf], and J ̃ being the rate of on-specification production (in units of product per unit time), defined at each point in time as J ̃(x(t , u), u(t )) ⎧ J(x(t , u), u(t )), if q(x(t , u), u(t )) ≤ 0, =⎨ otherwise ⎩ 0, ⎪

⎪

Figure 2. Change in the sequence of modes visited by a perturbation in the parameter.

(2)

By embedding the quality constraints in the objective function, any portion of the product that violates the quality constraints will not be considered. This way, the optimizer would naturally drive the operation toward maximum on-specification production. Computational Aspects. The objective function can be treated as an auxiliary state variable with the time derivative xȧ (t , u) = J (̃ x(t , u), u(t ))

to be switched in the same order for q(x(t, p̂ + ϵ), p̂ + ϵ) for all |ϵ| ≤ δ with δ > 0, although the times at which these modes switch could naturally change. In case a sufficiently large perturbation in p̂ changes the sequence of modes visited (as shown in the bottom plot of Figure 2), the derivatives of the solution of the hybrid system with respect to p̂ may not exist at some intermediate value of p̂. As a result, the QaO approach, although being an ideal formulation, could potentially lead to a nonsmooth (or even discontinuous) optimization problem that cannot be solved reliably using conventional derivative-based algorithms. 2.2.2. Quality-as-Constraint Approach. The problem of potential nondifferentiability in the QaO approach would be resolved if the objective function is defined as the total production (on-specification and off-specification) and the quality specifications participate as explicit constraints instead. One way of doing this would be through enforcing the quality constraints over the entire time horizon, as below:

(3)

Because of the discrete behavior in eq 1, the dynamics in eq 3 are discontinuous, which gives rise to a hybrid system.10 Hybrid systems can be seen as a sequence of switching dynamic models with switches triggered according to certain discrete events. In ̃ ̃ eqs 2 to 3, J(x(t, u), u(t)) = J(x(t, u), u(t)) and J(x(t, u), u(t)) = 0 are the alternate dynamics or modes, and switches are triggered according to the discrete event defined by the quality constraints. These modes are illustrated schematically in Figure 1 for a hypothetical scenario. The system starts and continues in Mode 1 for as long as the quality constraints remain satisfied. Then, it switches to Mode 2 due to violation of the constraints, and finally switches back to Mode 1 when the constraints are again satisfied. Hybrid systems pose major challenges regarding the existence of sensitivity information that is required for derivativebased optimization algorithms. Specifically, the solution of the hybrid system may not be differentiable with respect to

∫ u∈< 0

max s.t.

tf

J(x(t , u), u(t )) dt

x(̇ t , u) = f(x(t , u), u(t )),

s(x(t , u), u(t )) ≤ 0, 11346

(P5)

∀ t ∈ (0, tf ]

∀ t ∈ [0, tf ]

(P5.1) (P5.2)


Article

Industrial & Engineering Chemistry Research q(x(t , u), u(t )) ≤ 0,

∀ t ∈ [0, tf ]

On the basis of the above arguments, it is possible to define an interval containing the turnpike, over which all the production is on-specification. This interval is called optimal on-specif ication interval, and denoted by [ton, toff] ⊂ [0,tf], with ton being the time at which the production becomes on-specification, and toff being the time at which the production becomes off-specification again. The on-specification and total productions are equal over [ton, toff], and can be used interchangeably. Moreover, the quality constraints are guaranteed to hold over [ton, toff]. As a result, eq P4 can be approximated by replacing the on-specification production J ̃ over [0,tf] with the total production J over [ton, toff] in its objective function, and explicitly enforcing the quality constraints over [ton, toff]. Since the optimal on-specification interval itself is not known a priori, it is included in the optimization variables. The new formulation can then be presented as follows:

(P5.3)

h(x(tf , u), u(tf )) ≤ 0

(P5.4)

x(0, u) = x 0

(P5.5)

A major problem with eq P5 is that it may end up being infeasible because the quality constraints in a complex process are likely to be violated during major transients such as startup and shutdown. Even if it is feasible, the requirement to meet the quality constraints at all times could be conservative compared to what could be achieved through eq P4. For example, consider a worst case in which the only way to meet the impurity thresholds would be to stop the production.

3. PROPOSED FORMULATION In the following, a new formulation is proposed that addresses the above issues while also avoiding potential nondifferentiability of the QaO formulation. It is inspired by turnpike theory,15 which characterizes the optimal evolution of a dynamic system in optimal control problems. It was initially investigated for optimal control problems with convex cost functions (in a minimization case). In such a particular case and provided a sufficiently long operation time, it can be established that the optimal operation is formed of a transient operation, followed by an approximately steady operation in an ϵ-neighborhood of the optimal steady state, followed by another transient operation.16 The approximate steady operation is often called the turnpike, as depicted in Figure 3. The extension of the turnpike

max

u ∈ F̂inB , which means the process is first run in the on-specification mode, and then switches to the off-specification mode. Consequently, the sequence of modes changes when exceeding the point F̂inB , thereby violating a sufficient condition for differentiability described in subsection 2.2.1. 4.1.3. Formulations III and IV. Formulation III solves the following optimization problem:

solver raised derivative-related errors and failed to converge for this problem, for example, with four uniform time stages and flow rate profiles initialized at the midpoint of their lower and upper bounds. To further investigate the situation, the objective function is evaluated for a large number of constant-over-time FinA and FinB , and the results are plotted in Figure 8. The left plot shows the total on-specification production as a function of both inlet flow rates, while the right plot shows a cross section of the 3D graph at a nominal FinA along with the total on- and off-specification production as a reference. It is easy to see, especially in the right plot, that the objective function is discontinuous over the domain of the optimization variables. From the right plot, it is observed that the total production increases monotonically by using a higher flow rate of B. On the other hand, the on-specif ication production stops following this trend and jumps to a very low level after F̂inB ≈ 1.17 × 10−3. From a process viewpoint, the trend is explained by the fact that an excess of B favors the side reaction producing the undesired species I, thereby raising CI above its threshold and making the product off-specification. A more precise explanation is also possible by looking into the sequence of modes for different values of FinB . As shown in Figure 9, the impurity level stays within the threshold during the entire operation for FinB = 1.13 × 10−3< F̂inB , which means the process is run entirely

max in in

FA , FB

∫0

tf

F out(t )C P(t ) dt

s.t. dynamic model (10) V (tf ) − V0 ≤ 0 C I(t ) ≤ 0.1,

∀ t ∈ [0, tf ]

initial conditions (11) FAin(t ) ∈ [0, 0.01], FBin(t ) ∈ [0.002, 0.01], a.e. in [0, tf ] (P12)

It turns out that eq P12 is infeasible for any number of time stages considered in the profile discretization. This means it is not possible to bring the volume back to its initial value at t = tf and still keep the CI level within the threshold. Unlike the A inlet, the B inlet cannot be closed fully during shutdown, thereby making the excess of B and high concentration of I inevitable. Now, the proposed formulation IV is employed as follows: max

FAin , FBin , τ1, τ2 , τ3∈ [0, t f ]

∫1

2

F out(t ′)C P(t ′) dt ′

s.t. transformed dynamic model, ∀t ′ ∈ (0,3] V (3) − V0 ≤ 0 C I(t ) ≤ 0.1,

∀ t ∈ [1, 2]

τ1 + τ2 + τ3 = tf initial conditions (11) FAin(t ) ∈ [0, 0.01], FBin(t ) ∈ [0.002, 0.01], a.e. in [0, 3] (P13)

Figure 9. Profile of the impurity concentration at a constant FinA = 0.001 and for different values of FinB ; the remaining operation interval has been truncated in the interest of clarity.

where the transformed dynamic model is obtained from the transformation of eq 10 into the t′ domain (see eq P7.1). Notice 11352


Article

Industrial & Engineering Chemistry Research

Figure 10. Optimal profiles for eq 10 based on Formulation IV. The units for the concentrations and volume are mol L−1 and L, respectively.

the two formulations are plotted in Figure 11. For Formulation IV, the on-specification product increases continuously except at the end of the process. As a result, only 4% of the product is off-specification. On the other hand, Formulation I leads to 23% off-specification product, as can be seen from the right plot. Notice that most of the off-specification is produced during the startup. 4.1.4. Shorter Operation. It is expected that for a shorter operation time, Formulation IV is even more advantageous over Formulation I as the steady state becomes less dominant. To verify this, eq 10 is revisited with tf = 25 min. This change has no effect on the results obtained from Formulation I, namely the minimum startup/shutdown times as well as the optimal transient profiles. For Formulation IV, the optimal onspecification interval is obtained as [ton * , toff * ] = [ 0, 22.3], which corresponds to (τ1*, τ2*, τ3*) = (0, 22.3, 2.7) in the transformed model. Also, the optimal control and state profiles are given in Figure 12, where the control profiles are discretized in the same manner as in the previous case. Note also that the optimization variables are initialized with the same numbers as for eq P13. Here, the control profiles reach their optimal steady-state values at t ≈ 7.5 min. However, the turnpike is only reached at t ≈ 16 when the state variables also reach their optimal steady state. The process moves away from the turnpike upon lowering the flow rates at t = 22 and starting the shutdown phase. The total and on-specification productions from the two formulations are depicted in Figure 13. Formulation I results only in 0.17 mol on-specification product out of a total production of 0.33 mol (i.e., 49% offspecification). On the other hand, the off-specification for Formulation IV is only 9%. In comparison with Formulation I, Formulation IV improves the on-specification production by 69%, which is much more than the improvement of 23% achieved in the base case with tf = 50. This further emphasizes the advantage of the proposed formulation for processes with short operation times. For convenience, the comparison of different formulations applied to the CSTR example are summarized in Table 2. 4.2. Van de Vusse Reactor. This example is adapted from the jacketed nonisothermal Van de Vusse reactor in Rothfuss et al.,28 in which the following reactions take place:

that the same variable symbols are used in the transformed formulation for the sake of notational simplicity. The control profiles are parametrized using the special discretization scheme described in subsection 3.1.1, with three uniform time stages in each of the intervals [0, 1], [1, 2], and [2, 3]. The initial guesses of (τ1, τ2, τ3) = (0, 20, 5) are used to start the optimization. Also, both flow rate profiles are initialized at their optimal steady-state values as obtained previously. The optimal control and state profiles in the original t domain are shown in Figure 10. It is interesting to notice that for much of the operation interval, the process stays close to the optimal steady state obtained earlier in subsection 4.1.1. This is a manifestation of the turnpike emerging in the optimal solution. Note, however, that information about the optimal steady state or a turnpike was not included in the problem formulation. Although it may be hard to confirm visually, the optimal steady state is reached around t ≈ 25 min, that is, much later than in the case of Formulation I. This is quite expected as a minimum startup time is not the objective anymore. Therefore, the B flow rate need not be pushed to its highest value in order to fill up the reactor quickly (unlike in Formulation I; see Figure 6). Thus, the excess of B in the reactor remains under control, so does the impurity concentration during the entire startup. The shutdown is triggered at t ≈ 47.3 min by lowering the inlet flows from their optimal steady state, which is the same as that obtained from Formulation I. The optimal onspecification interval is also obtained as [t*on, t*off] = [0,47.3], which corresponds to (τ*1 , τ*2 , τ*3 ) = (0,47.3,2.7) in the transformed model. The optimal operation yields 0.65 mol on-specification product, which is 23% more than what would be attained from Formulation I. The total and on-specification productions from

k1

k1

A→B→C k2

2A → D

where B is the desired product, and C and D are undesired byproducts. The process is depicted in Figure 14, which has

Figure 11. Total and on-specification productions for Formulation IV (left) and Formulation I (right) in eq 10. 11353


Article


Figure 12. Optimal profiles for eq 10 based on Formulation IV. The units for the concentrations and volume are mol L−1 and L, respectively.

CȦ (t ) = −k1(T (t ))CA(t ) − k 2(T (t ))CA2(t ) C in − CA(t ) in + A F (t ), V (t ) Ċ B(t ) = k1(T (t ))(CA(t ) − C B(t )) − ĊC(t ) = k1(T (t ))C B(t ) − Figure 13. Total and on-specification productions for Formulation IV (left) and Formulation I (right) in eq 10 with tf = 25 min.

Ċ D(t ) =

formulation I II

on-spec. product

0.53 (mol) not converged III infeasible IV (proposed) 0.65 (mol)

on-spec. product

improvement

N/A

0.17 (mol) N/A

N/A

N/A 23%

N/A 0.29 (mol)

N/A 69%

C (t )F in(t ) 1 k 2(T (t ))CA2(t ) − D , 2 V (t ) T in − T (t ) , V (t )

Tċ (t ) = β(T (t ) − Tc(t )) − γPc(t ), V̇ (t ) = F in(t ) − F out(t ),

tf = 25 min

improvement

CC(t )F in(t ) , V (t )

Ṫ(t ) = h(t ) + α(Tc(t ) − T (t )) + F in(t )

Table 2. Comparison of the on-Specification Production for Different Formulations Applied to Eq 10. The Improvement Is Calculated with Respect to Formulation I tf = 50 min

C B(t )F in(t ) , V (t )

F out(t ) ≡ η V (t ) , h(t ) ≡ −δ(k1(T (t ))(CA(t )ΔHAB + C B(t )ΔHBC) + k 2(T (t ))CA(t )2 ΔHAD), ⎛ ⎞ −E1 k1(T (t )) ≡ k 0,1exp⎜ ⎟, ⎝ T (t ) + 273.15 ⎠ ⎛ ⎞ −E2 k 2(T (t )) ≡ k 0,2exp⎜ ⎟ ⎝ T (t ) + 273.15 ⎠ (12)

with the initial conditions (CA(0), C B(0), CC(0), C D(0)) = (0, 0, 0, 0) mol m−3, (T (0), Tc(0)) = (108, 107.7)°C, V (0) = 10−3 m 3 (13)

The model parameters are given in Table 3. Note that the dependence of the normalized heat transfer coefficients α and β on the volume V are ignored for simplicity. This assumption is particularly valid for α if the jacket surrounds the side area of the reactor only, as illustrated in Figure 14. The operation takes tf = 0.8 h. The product is considered on-specification at each time t if CC(t) ≤ 500 mol m−3 and CD(t) ≤ 350 mol m−3, and off-specification otherwise. Moreover, safety reasons stipulate that the reactor temperature T must not exceed 111 °C at any

Figure 14. Van de Vusse reactor considered in subsection 4.2.

been modified to accommodate varying volume, and is modeled as follows: 11354


Article

Industrial & Engineering Chemistry Research Table 3. Model Parameters for Eq 12 1.29 × 1012 9.04 × 106 9758.3 8560 4.2 −11 −41.85

k0,1 k0,2 E1 E2 ΔHAB ΔHBC ΔHAD

h−1 m3 (mol·h)−1 K K kJ mol−1 kJ mol−1 kJ mol−1

α β γ δ η Tin CinA

30.828 86.688 0.1 3.52 × 10−4 30 104.9 5.10 × 103

h−1 h−1 K kJ−1 m3 K kJ−1 m3/2 h−1 °C mol m−3

Figure 15. Optimal profiles for the startup time minimization of eq 12. The simulation has been extended past the vertical lines (t = t*ss ) to validate the steady state. See the text for the units of the quantities.

obtained as J* = 39373 mol h−1, which is attained at Fin, ss * = 40 m3 h−1, P*c,ss = 2873 kJ h−1, C*ss = (2922,984,500,347) mol m−3, (Tss*,Tc,ss * ) = (110,107) °C, and Vss* = 1.78 m3. Now, the minimum startup time is obtained from

time. Also, the reactor volume V must be 0.01 m3 or less at the end of the operation. The optimization variables are taken as the flow rate Fin and the cooling power Pc, which can vary within [0, 40] m3 h−1 and [0,4000] kJ h−1, respectively. In the sequel, the results of applying Formulations I and IV are presented. 4.2.1. Formulation I. First, the steady-state optimization is formulated as

min

F in , Pc , tss ∈ [0, t f ]

s.t. dynamic model (12) x(tss) − x *ss = 0

J ≔ FssoutC B,ss

max

Fssin , Pc ,ss , Css, Vss , Tss , Tc ,ss

T (t ) − 111 ≤ 0,

s.t. steady‐state model, Fssin ∈ [0, 40],

CC,ss ∈ [0, 500],

F in(t ) ∈ [0, 40], Pc(t ) ∈ [0, 4000], a.e. in [0, tf ]

C B,ss ∈ [0, 2000],

(P15)

where x ≔ (C, T, Tc, V) represents the vector of the state variables, and xss* is the vector of the corresponding optimal steady-state values. The time transformation described in subsection 2.1.1 with a discretization of the control profiles over 20 uniform time stages is used. Also, the optimization variables are initialized as 4 for the startup time tss and as the midpoint of the lower and upper bounds for the control profiles. The minimum time to the optimal steady state is computed as t*ss = 0.64 h. The optimal startup profiles for the control and state variables are shown in Figure 15. The steady state at tss* is validated by continuing the simulation beyond tss*, while setting

C D,ss ∈ [0, 350],

Tss ∈ [0, 111], Tc,ss ∈ [0, 115],

∀ t ∈ [0, tss]

initial conditions (13)

Pc,ss ∈ [0, 4000],

CA,ss ∈ [0, 3000],

Vss ∈ [0, 2],

tss

(P14)

where Css = (CA,ss,CB,ss,CC,ss,CD,ss), and the steady-state model refers to eq 12 with the time derivatives set to zero. Here again, the upper bounds on CA,ss, CB,ss, Tc,ss, and Vss are considered only because the global solver requires them to proceed. At convergence, the solution is checked to make sure these bounds are not active. The maximum steady-state production rate is 11355


Article


Figure 16. Optimal profiles for the shutdown time minimization of eq 12. The simulation includes the preceding steady-state phase t ≤ tf − tsd * as a reference. See the text for the units of the quantities.

Fin and Pc to their optimal steady-state values. For the entire startup, the flow rate is optimally operated at its maximum value, which is also its optimal steady-state value. This allows for a quick filling of the reactor. On the other hand, the cooling power undergoes moderate fluctuations. It is increased initially in order to prevent the reactor temperature from exceeding the specified threshold 111 °C. Regarding the impurity concentrations, while CD remains below the on-specification limit 350 mol m−3, CC exceeds its on-specification limit of 500 mol m−3 for t ∈ [0.32,0.54], leading to off-specification production during this interval (see Figure 19 and related discussions). Finally, the minimum shutdown time is obtained from min

F in , Pc , tsd ∈ [0, t f ]

on-specification product. This means 41% of the product ends up off-specification given a total production of 27731 mol. 4.2.2. Formulation IV. The proposed formulation reads:

∫ F , P , τ , τ , τ ∈ [0, t ] 1 in

c 1 2 3

2

F out(t ′)C B(t ′) dt ′

f

s.t. transformed dynamic model, ∀t ′ ∈ (0,3] V (3) − 10−2 ≤ 0 T (t ′) ≤ 111,

tsd

∀ t ′ ∈ [0, 3]

CC(t ′) ≤ 500,

∀ t ′ ∈ [1, 2]

C D(t ′) ≤ 350,

∀ t ′ ∈ [1, 2]

τ1 + τ2 + τ3 = tf

s.t. dynamic model (12)

initial conditions (13)

V (tsd) − 10−2 ≤ 0 T (t ) − 111 ≤ 0,

max

F in(t ′) ∈ [0, 40], Pc(t ′) ∈ [0, 4000], a.e. in [0, 3]

∀ t ∈ [0, tsd]

(P17)

x(0) − x *ss = 0

where the transformed dynamic model is obtained from the transformation of eq 12 into the t′ domain (see eq P7.1)). The control profiles are parametrized using the special discretization scheme described in subsection 3.1.1, with two uniform time stages in each of the intervals [0,1], [1,2], and [2,3]. The initial guesses of (τ1, τ2, τ3) = (0,0.7,0.1) are used to start the optimization. Also, the Fin and Pc profiles are initialized at 40 and 2600, respectively. The optimal control and state profiles in the t′ domain are shown in Figure 17. The interval [1,2] corresponds to the optimal on-specification interval. Notice that the state profiles in the first interval [0,1] are constant. This is explained by the optimal values for τ, which are obtained as (τ*1 ,τ*2 ,τ*3 ) = (0,0.72,0.08) . Since τ*1 = 0, the right-hand side of the transformed dynamic model is zero for t′ ∈ [0, 1] (see eq P7.1), and thus, the state variables remain unchanged. This is equivalent to the first off-specification interval having a zero

F in(t ) ∈ [0, 40], Pc(t ) ∈ [0, 4000], a.e. in [0, tf ] (P16)

Similar to the previous case, the time transformation described in subsection 2.1.1 is used, and then, the control profiles are discretized over 10 uniform time stages. The optimization variables are initialized as in eq P15. The minimum shutdown time is obtained as tsd * = 0.085 h. The optimal shutdown profiles are shown in Figure 16, where part of the preceding steady-state phase is also included as a reference. The shutdown begins at t = 0.715 h, which results in a short steady-state duration of approximately 0.08 h. The quality constraints are violated during the entire shutdown. The concatenation of the obtained startup, steady-state, and shutdown profiles yields 16241 mol 11356


Article


Figure 17. Optimal profiles for eq 12 based on Formulation IV in the t′ domain. See the text for the units of the quantities.

Figure 18. Optimal profiles for eq 12 based on Formulation IV. See the text for the units of the quantities.

stays on it until t = 0.72. Notice the cooling power is increased at t = 0.72 in order to maintain the reactor temperature within the threshold during the shutdown phase. However, the impurity concentrations rise well beyond their limits in this phase. An on-specification production of 25377 mol is obtained in this case (92% of the total production), which is 56% more than what would be produced from using Formulation I. The total and on-specification productions from both formulations are plotted in Figure 19. Interestingly, Formulation IV results

duration in the t domain, rendering the control and state profiles in t′ ∈ [ 0,1] physically meaningless. From eq 7, the optimal onspecification interval is calculated as [ton, toff] = [0, 0.72]. The optimal control and state profiles in the original time domain are given in Figure 18. The Fin profile reaches its optimal steady state at t = 0, while the Pc profile reaches a neighborhood of its optimal steady state at t = 0.36. From the state profiles, the process begins its turnpike at t ≈ 0.40, with the approximation accounting for the process dynamics, and 11357


Article


(2) Benyahia, B.; Lakerveld, R.; Barton, P. I. A Plant-Wide Dynamic Model of a Continuous Pharmaceutical Process. Ind. Eng. Chem. Res. 2012, 51, 15393−15412. (3) Mascia, S.; Heider, P. L.; Zhang, H.; Lakerveld, R.; Benyahia, B.; Barton, P. I.; Braatz, R. D.; Cooney, C. L.; Evans, J. M. B.; Jamison, T. F.; Jensen, K. F.; Myerson, A. S.; Trout, B. L. End-to-End Continuous Manufacturing of Pharmaceuticals: Integrated Synthesis, Purification, and Final Dosage Formation. Angew. Chem., Int. Ed. 2013, 52, 12359− 12363. (4) Leuenberger, H. Eur. J. Pharm. Biopharm. 2001, 52, 289−296. (5) Goršek, A.; Glavič, P.; Senčar, P. Optimal process design for specialty products. Comput. Chem. Eng. 1992, 16 (Suppl. 1), S321− S328. (6) Zhang, H.; Lakerveld, R.; Heider, P. L.; Tao, M.; Su, M.; Testa, C. J.; D'Antonio, A. N.; Barton, P. I.; Braatz, R. D.; Trout, B. L.; Myerson, A. S.; Jensen, K. F.; Evans, J. M. B. Application of Continuous Crystallization in an Integrated Continuous Pharmaceutical Pilot Plant. Cryst. Growth Des. 2014, 14, 2148−2157. (7) Lakerveld, R.; Benyahia, B.; Braatz, R. D.; Barton, P. I. Modelbased design of a plantwide control strategy for a continuous pharmaceutical plant. AIChE J. 2013, 59, 3671−3685. (8) Lakerveld, R.; Benyahia, B.; Heider, P. L.; Zhang, H.; Wolfe, A.; Testa, C. J.; Ogden, S.; Hersey, D. R.; Mascia, S.; Evans, J. M. B.; Braatz, R. D.; Barton, P. I. The Application of an Automated Control Strategy for an Integrated Continuous Pharmaceutical Pilot Plant. Org. Process Res. Dev. 2015, 19, 1088−1100. (9) Heider, P. L.; Born, S. C.; Basak, S.; Benyahia, B.; Lakerveld, R.; Zhang, H.; Hogan, R.; Buchbinder, L.; Wolfe, A.; Mascia, S.; Evans, J. M. B.; Jamison, T. F.; Jensen, K. F. Development of a Multi-Step Synthesis and Workup Sequence for an Integrated, Continuous Manufacturing Process of a Pharmaceutical. Org. Process Res. Dev. 2014, 18, 402−409. (10) Barton, P. I.; Lee, C. K. Modeling, Simulation, Sensitivity Analysis, and Optimization of Hybrid Systems. ACM Trans. Model. Comput. Simul. 2002, 12, 256−289. (11) Teo, K.; Goh, C.; Wong, K. A Unified Computational Approach to Optimal Control Problems; Pitman Monographs and Surveys in Pure and Applied Mathematics; Longman Scientific and Technical: 1991. (12) Chachuat, B. Nonlinear and Dynamic Optimization: From Theory to PracticeIC-32: Spring Term 2009; Polycopiés de l’EPFL: 2009. (13) Vassiliadis, V.; Sargent, R.; Pantelides, C. Solution of a class of multistage dynamic optimization problems. 2. Problems with path constraints. Ind. Eng. Chem. Res. 1994, 33, 2123−2133. (14) Galán, S.; Feehery, W. F.; Barton, P. I. Parametric sensitivity functions for hybrid discrete/continuous systems. Appl. Numer. Math. 1999, 31, 17−47. (15) McKenzie, L. Turnpike theory. Econometrica 1976, 44, 841− 865. (16) Trélat, E.; Zuazua, E. The turnpike property in finitedimensional nonlinear optimal control. J. Differ. Equ. 2014, 258, 81− 114. (17) Zaslavski, A. J. Turnpike Properties in the Calculus of Variations and Optimal Control; Nonconvex Optimization and Its Applications; Springer: 2006; Vol. 80. (18) Zaslavski, A. J. Turnpike Phenomenon and Infinite Horizon Optimal Control; Springer Optimization and Its Applications; Springer: 2014; Vol. 99. (19) Faulwasser, T.; Korda, M.; Jones, C. N.; Bonvin, D. IEEE 53rd Annual Conference on Turnpike and dissipativity properties in dynamic real-time optimization and economic MPC; Decision and Control (CDC); 2014; 2734−2739; DOI: 10.1109/ CDC.2014.7039808 (20) Carlson, D.; Haurie, A.; Leizarowitz, A. Infinite Horizon Optimal Control: Deterministic and Stochastic Systems; Springer-Verlag: 1991. (21) Lee, H.; Teo, K.; Rehbock, V.; Jennings, L. Control parametrization enhancing technique for optimal discrete-valued control problems. Automatica 1999, 35, 1401−1407. (22) Khan, K.; Yunt, M.; Barton, P. Some New Results on Sensitivity Analysis of Hybrid Systems. J. Nonlinear Syst. Appl. 2011, 2, 64−72.

Figure 19. Total and on-specification productions for Formulation IV (left) and Formulation I (right) in eq 12.

in a steady production of on-specification product, whereas Formulation I leads to intermittent periods of off-specification production during the startup phase. Intermittent formation of off-specification product in a continuous process is highly undesirable because of difficulties associated with isolating such product and separating it from the continuous production line.

5. CONCLUSIONS A new formulation for optimal operation in campaign CM is proposed. Owing to the typically short operation times, the proposed formulation emphasizes maximizing on-specification production explicitly during the entire operation rather than minimizing the startup and shutdown times. However, it is shown that such an objective function leads to a hybrid dynamic system, which is nondifferentiable or even discontinuous in general since the sequence of modes visited can change with different values of the optimization variables. Therefore, a number of reformulations are presented to fix the sequence of modes during optimization. First, the operation window is split into three variable-length intervals, where only the middle interval results in on-specification product. The middle interval contains the steady-state turnpike if it exists for the system under study. Moreover, a time transformation together with a special control discretization scheme is employed to facilitate numerical solution while maintaining differentiability of the optimization problem. It is shown that the proposed formulation yields significantly improved results in terms of on-specification production compared to the conventional formulations, especially when the operation window is short as will often be the case in continuous pharmaceutical industries.

■

AUTHOR INFORMATION

Corresponding Author

*Phone: +1 617 253 6526. Fax: +1 617 258 5042. E-mail: pib@ mit.edu. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS Financial support from Novartis-MIT Center for Continuous Manufacturing is gratefully acknowledged. REFERENCES

(1) Schaber, S. D.; Gerogiorgis, D. I.; Ramachandran, R.; Evans, J. M. B.; Barton, P. I.; Trout, B. L. Economic Analysis of Integrated Continuous and Batch Pharmaceutical Manufacturing: Ind. Eng. Chem. Res. 2011, 50, 10083−10092. 11358


Article

Industrial & Engineering Chemistry Research (23) Wächter, A.; Biegler, L. T. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Math. Program. 2006, 106, 25−57. (24) Tolsma, J. E.; Barton, P. I. DAEPACK: An Open Modeling Environment for Legacy Models. Ind. Eng. Chem. Res. 2000, 39, 1826− 1839. (25) Tolsma, J. E.; Barton, P. I. Hidden Discontinuities and Parametric Sensitivity Calculations. SIAM J. Sci. Comput. 2002, 23, 1861−1874. (26) Tawarmalani, M.; Sahinidis, N. V. polyhedral branch-and-cut approach to global optimization. Math. Program. 2005, 103, 225−249. (27) GAMS−A User’s Guide; GAMS Software GmbH: 2014; http:// www.gams.com/help/topic/gams.doc/userguides/GAMSUsersGuide. pdf. (28) Rothfuss, R.; Rudolph, J.; Zeitz, M. Flatness based control of a nonlinear chemical reactor model. Automatica 1996, 32, 1433−1439.

11359


Optimal Campaign Continuous Manufacturing - ACS Publications

Recommend Documents