Optimal Dynamic Continuous Manufacturing of Pharmaceuticals with

Article pubs.acs.org/IECR

Cite This: Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Optimal Dynamic Continuous Manufacturing of Pharmaceuticals with Recycle Michael Patrascu†,‡ and Paul I. Barton*,† †

Downloaded via GUILFORD COLG on July 19, 2019 at 01:32:12 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

Process Systems Engineering Laboratory, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, United States

ABSTRACT: Continuous manufacturing of pharmaceuticals is challenged by short campaigns due to the low volume of some products. This calls for optimization of the entire dynamic campaign to achieve the highest overall yields and productivity and to minimize waste. In this paper, we investigate how recycles could be incorporated into the design of an effective end-to-end continuous manufacturing process using nonsmooth dynamic simulation and optimization. We explore two configurations for recycle that use an organic solvent nanofiltration membrane to concentrate and recycle the mother liquor back to either the crystallizer or the reactor and find the optimal dynamic operation given a short manufacturing horizon. The campaign performance is investigated in terms of the overall yield and the on-specification productivity. It is demonstrated that recycling back to the crystallizer does increase the yield for the same overall productivity; however, recycling to the reactor significantly improves both the yield and the productivity.

■

INTRODUCTION In the pursuit of more reliable, cost-effective, and sustainable processes, the pharmaceutical industry is considering to switch to continuous flow processes from the traditional batch-wise manufacturing. Continuous manufacturing (CM) allows one to integrate intensified processes, such as microreactors and microseparators, creating a more efficient, cost-effective, and flexible process that is much less challenging to scale up. These units and manufacturing techniques also have a smaller footprint and hold the potential to reduce wastes and improve overall yields.1−3 Continuous manufactuing of pharmaceuticals is likely to be characterized by short campaigns (a few days to weeks), making the relatively long and unproductive start-up and shutdown periods more significant compared to other industries. This is one of the most challenging problems CM of pharmaceuticals is facing to make such technologies economically attractive. Therefor, optimal campaigns are of interest, in which the process is optimized dynamically over the full time horizon,4,5 as opposed to searching only an optimal steady state. © XXXX American Chemical Society

At the Novartis-MIT Center for Continuous Manufacturing, a pilot plant was developed and operated, performing integrated synthesis, purification, formulation, and tableting in an entirely continuous mode.6 The manufacturing approach, according to the traditional paradigm of continuous processes, was to operate the process at a steady state.7,8 The steady state of each operating unit was set based on their individual performance.9−12 The integrated design, short campaign time, and productivity constraints were not taken into account when determining the optimal operating procedures. Furthermore, recycle streams were eventually not implemented in the demonstration process, probably due to technical difficulties. Computational tools are extremely important to improve yields, productivity, energy consumption, and other key Special Issue: Dominique Bonvin Festschrift Received: Revised: Accepted: Published: A

January 31, 2019 June 13, 2019 June 18, 2019 June 18, 2019 DOI: 10.1021/acs.iecr.9b00646 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 1. Process flow diagram of the CM plant highlighting the recycle loop in red. Control valves used in the optimization, u1−u6, are highlighted, mem as well as the three-way valve, < R1 , directly. mem , which directs the recycle to either Cr2 or to R1; u6 controls the recycle ratio, R

process attributes, especially where low margins and critical operating constraints exist. Recently, various efforts have been made to study the dynamics of individual process units, e.g., crystallization, granulation, and blending.13−16 These studies focused on the time to steady state, unit operations interactions, the effect of process upsets, and process optimization. We recently studied the optimal performance of a fully integrated end-to-end continuous pharmaceutical plant over the full campaign, rather than focusing on the steady state alone.5,17 Here, we build upon this model to study the potential feasibility and benefits of using organic solvent nanofiltration (OSN) membranes to recycle the concentrated mother liquor and potentially improve the yield and productivity for a short campaign. Obviously, effective implementation of process systems engineering tools (i.e., modeling, simulation, optimization, and control) is based on fundamental development and quantification of kinetic and thermodynamic parameters related to the process units. When such data do not exist, such as during the development of an entirely new chemical route, most of the early development work will address that, ideally using design of experiments (DOE)18 and quality by design (QbD)19 methods. This was the approach taken at the Novartis-MIT Center for CM. Another example for such work is the scale-up of continuous flow systems by Lonza,20 where implementation strategies for production campaigns of up to tons of pharmaceutical chemicals are discussed, based on Lonza projects. The parametric uncertainty should be addressed during the implementation of the in-line control system.21,22

Recycling is largely avoided from pharmaceuticals manufacturing processes. A few papers have explored the use of recycle loops in process design, both theoretically and experimentally; however, most were focused on the steadystate performance. Griffin et al.23 have proposed a continuous crystallization process to produce active pharmaceutical ingredients (APIs) with a small mean size and narrow distribution by recycling the large crystals from a classifier to a dissolver along with the mother liquor from a filter. They have used kinetic models developed for aspirin as a case study to explore the stability and the regions of sustained limit cycles of the crystallization process and have demonstrated that with simple changes to the crystallizer operation a product of a small mean size and a narrow distribution can be attained directly. The effect of recycle ratio on both crystal purity and process yield was evaluated by Alvarez et al.24 in a three-stage mixed-suspension, mixed-product removal (MSMPR) system. An evaporator was used in the recycle line to remove the solvent. They showed that the yield of the process can be significantly improved by increasing the recycle ratio while the crystal purity decreases. Benyahia et al.25 studied a theoretical recycle from the crystallization step back to the reactor using evaporation to remove the solvent using a dynamic model and simulations, focusing on the steady-state performance, again showing how yields can be improved. Recycling to a mixer in a direct compaction process was studied by Boukouvala et al.26 They have demonstrated the benefit of recycle to be lower sensitivity of the final API concentration in the tablets to feed rate fluctuations. They pointed out that the trade-off is the increase in the cost of the system and the slow response of the B

DOI: 10.1021/acs.iecr.9b00646 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research system to online control actions. Recently, Wang et al.27 showed that the attainable size and design space can be enlarged significantly by extending conventional crystallization with membranes, assuming total rejection of the solutes (product and impurity) by the membrane. In particular, larger crystals or shorter residence times can be achieved. Membranes allow solvent removal and provide an additional method for supersaturation generation, thereby increase the process flexibility. A sequence of batch cooling crystallization experiments was carried out to investigate the effect of mother liquor recycle on the buildup of impurities in a work by Keshavarz et al.28 The results showed that the rate of impurity buildup decreases with an increasing number of cycles until the amount of the impurity remains constant. The number of cycles required to reach the steady state increases when larger fractions of mother liquor were reused. They showed how to optimize recycle conditions to enhance process efficiency by reducing product and solvent waste. It is likely that the longer transients introduced by recycle loops will negatively impact short manufacturing campaigns. Also, due to the trade-off between yield and product quality (impurities buildup), it is not clear whether a recycle operation will show improvement over a campaign without recycle that has been optimized for yield and is producing a product with the quality already at the specification level. It is thus the goal of this work to demonstrate the benefits obtained by such recycle configurations in continuous pharmaceuticals manufacturing. We consider two recycle loop alternatives: recycling the concentrated mother liquor to the crystallizers for further separation and recycling it to the reactor for further conversion. Both of these alternatives may result in undesired increased levels of impurities in the final product. In the following section, we describe the chemical process explored, with the two recycle alternatives, followed by a detailed mathematical formulation of the process model, using nonsmooth differential-algebraic equations (DAEs).29 Thereafter, we investigate the simulation results for the two configurations, based on a naive approach, which includes simulations of specific operating approaches. In the main part of the paper, the dynamic optimization formulation and the optimization results are presented and discussed, with emphasis on the improved performance predicted versus the naive simulations. The benefits of the two recycle configurations are compared to the optimal case without recycle.

water (S2), in which C2 and cat1 are dissolved to form the aqeuous phase. LLE1 separates the two phases into pur1 (the aqueous phase) and an organic phase, which continues downstream to purify C3, the intermediate compound, in crytallizers Cr1 and Cr2.11 Supersaturation is achieved both by cooling and by adding heptane as an antisolvent, AS, at a fixed ratio of 1:3 (AS to S1). The solids are separated in WF1 by a continuous wash and filtration unit. The mother liquor is split between a direct purge stream and a membrane (Mem), which concentrates the solution and recycles it back. The permeate from the membrane (the solvent) is then mixed with the purge stream (pur2). The filter cake passes into a well-mixed vessel, D1. Ethyl acetate is added to D1 to dissolve the crystals and bring the concentration of the intermediate down to the required level. The slurry from D1 is mixed with aqueous HCl in M3 to synthesize C4 in PFR R2. This is achieved by removing the Boc protecting group at room temperature, Reaction 2a. Other side reactions occur as well which produce the impurity I2, as stated in Reaction Set (2) C3 → C4 + Boc

(2a)

C1 → I 2 + Boc

(2b)

C4 → I 2 + C2

(2c)

NaOH is then added to neutralized the acid. C4 is separated together with the organic phase in LLE2 and diluted with ethyl acetate in M5. The aqueous phase is separated and purged in pur3. The API, aliskiren hemifumarate, is produced in a twostage reactive crystallization step in crystallizers Cr3 and Cr4. Fumaric acid (FA) is added to Cr3.10 Separation of the API crystals takes place in WF2, which are subsequently diluted in dilution tank D2. Silicon dioxide is added in M6 to increase the flowability of the powder. The slurry mixture is fed to a drum dryer (DD) to remove most of the solvent and to a tubular dryer, equipped with a screw, (SD) to complete the drying process. Eventually, the API product powder is blended in the extruder Ex with 6000 Da polyethylene glycol (PEG). A product is considered on-spec if it satisfies these critical quality attributes (CQA): (1) API concentration higher than 34 wt %, (2) impurity I2 concentration up to 0.3 wt %, (3) total impurities concentration up to 0.5 wt %, and (4) solvent S1 concentration up to 0.5 wt %.

■

■

NONSMOOTH DYNAMIC MATHEMATICAL MODEL We model the dynamic system by a semiexplicit generalized index-1 nonsmooth DAE system,29,30 represented by the general form of equation set 3. In this system of equations, f and g are continuous, not necessarily smooth, functions (i.e., not differentiable everywhere). A more comprehensive discussion of nonsmooth DAEs can be found in Stechlinski et al.29 Here, x are the differential variables, and y are the algebraic ones, with u being the control functions. The system is composed of 2116 equations, including 880 differential equations (eq 3a). Equations 3b are the algebraic equations. Notice that the nonsmooth functions g and f do not depend on time explicitly.

MANUFACTURING PROCESS DESCRIPTION In this paper, we investigate a case study inspired by a pilot plant, which was designed and operated at MIT,6 illustrated in Figure 1, with the main modification being the recycle loops. The main reagent, C1, the second reactant, C2 (which is considered as the solvent) and the catalyst, cat1, are fed to the first mixer (M1) at a fixed molar ratio of 1:10:1 of C1:C2:cat1. In one alternative, this stream is mixed with the concentrated mother liquor from the first crystallizers train (Cr2, see below), Figure 1. This stream enters a PFR reactor (R1), where the intermediate C3 and the impurity I1 are produced according to Reaction Set (1)12 C1 + C2 ↔ C3

(1a)

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

(3a)

C3 → I1

(1b)

0 = g(x(u , t ), y(u , t ), u(t ))

(3b)

x(u , 0) = x 0

(3c)

In mixer M2 two solvents are added: ethyl acetate (S1), which dissolves C1, C3, and I1 to form the organic phase, and C


Article

Industrial & Engineering Chemistry Research The full details of the equations are given elsewhere,17 and only the equations relevant for the recycle and crystallizers and given here. The model is implemented and simulated using the Jacobian software (RES Group Inc.31). For optimization we use DAEPACK32 with IPOPT.33 Mixer M1. The mass and component balances on the mixer, assuming that blending is rapid, are F C2(t ) = 10F C1(t )MWC2 /MWC1 F cat1(t ) = F C1(t )MWcat1/MWC1

V R1ρl

dwiR1 ,j dt

3 R1 R1 R1 + V R1 ∑ siR1 , l rl , j (t )MWi , i = 1 ,···, nc , j = 1 ,···, n l=1

(11)

wiR1 ,0 (t )

(4)

=

(5)

Y (t ) =

C1 C1 C2 C2 cat1 cat1 FinR1(t )wiR1 ,in(t ) = F (t )wi + F (t )wi + F (t )wi

r1R1(t , z) = k1R1wCR11 (t , z)ρl /MWC1

(8a)

r2R1(t , z) = k 2R1wCR13 (t , z)ρl /MWC3

(8b)

z)ρl /MWC3

ρl

k=1

FinCr1(t ) = F raf1(t )

(15)

raf1 Cr1 wiCr1 ,in (t ) = wi (t ), i = 1 ,···, nc

(16)

+ (1 − εwCr1(t ))wiCr1, s(t )), i = 1 ,···, ncCr1

(17)

MiCr1(t ) = M Cr1(εwCr1(t )wiCr1, l(t ) + (1 − εwCr1(t ))wiCr1, s(t )), i = 1 ,···, ncCr1

(8c)

(18)

ρmCr1(t ) = εCr1(t )ρl + (1 − εCr1(t ))ρs

(19)

εwCr1(t ) = εCr1(t )ρl /ρmCr1(t )

(20)

εCr1(t ) = 1 − kvμ3 (t )

(21)

M Cr1(t ) = ρmCr1(t )V Cr1(t )

(22)

ncCr1

∑ wiCr1,l(t ) = 1 i=1

(23)

where and are the mass flow rates in and out, respectively, εCr1 and εCr1 are the liquid mass fraction and w liquid volume fraction, respectively, which are calculated based on the crystal size distribution, eq 35. Here, ρl and ρs are the constant densities of the mother liquor and the solids, respectively, kv is the shape factor of the crystal, and μ3 is the third moment of the crystal size distribution (CSD). This model allows the holdup volumes to change with no upper limit, enabling a model-based design of the crystallizer vessels size. A distribution coefficient approach is used to calculate the solid and liquid phases mass fractions, wCr1,s and wCr1,l , i i respectively. Here, we assume that the only impurity inclusions in the intermediate crystals, C3, are I1 and C1, and that the FCr1 in

(9)

3

∑ siR1,k MWri kR1(t , z)

(14)

dMiCr1 Cr1 Cr1 Cr1, l (t ) = FinCr1(t )wiCr1 (t ) ,in (t ) − Fout (t )(εw (t )wi dt

R1 ∂wiR1 F R1(t ) ∂wi (t , z ) = − (t , z ) ∂t A ∂z

+

wCR11,in(t )MWC3

Using the mixed-suspension mixed-product removal (MSMPR) model, the dynamic mass and component balances are

where rR1 is the reaction rate of reaction i, z is the axial i coordinate of the reactor, ρl is the constant density of the liquid, and kR1 i are the rate constants. Numerical values for all the parameters used in the model are given in ref 17. Due to the large excess of C2 in the feed, the rate eq 8a is assumed to be first order.34 The reactor is assumed to be always full. The mass and component balances are R1 FinR1(t ) = Fout (t ) ≡ F R1(t )

(wCR13,out(t ) − wCR13,in(t ))MWC1

Crystallizers 1 and 2. Purification of the intermediate C3 is achieved in two consecutive crytallizers with antisolvent (AS) added to Cr2.6,8,11 The feed to Cr1 is the raffinate from LLE1

(7)

where Fi are the mass flow rates, nR1 c represent the number of components in R1, wji are the mass fractions of the raw materials (assumed pure), and MWi is the molecular weight of species i. In this unit the number of species, nR1 c , is 8. The flow rate of C1 in the feed, FC1, is directly manipulated (FC1(t) = u1(t)). The catalyst ratio in the R1 feed is based on results from previous studies.12,34 < R1 mem is the three-way valve controlling the recycle stream from the membrane to the reactor or the crystallizer. Its value is either 1 or 0, depending on which recycle configuration is tested: 1 for the reactor recycle and 0 for the crystallizer recycle. Reactor R1. The first reaction is conducted in a reactor tube, 2.7 L in volume, where plug flow is assumed. The intermediate, C3, is produced according to the reversible Reaction 1a, and the irreversible Reaction 1b produces the impurity I1. The reaction rates are calculated by eq 812

z) =

(13)

The instantaneous yield of the reactor, Y (t), is calculated by R1

k 3R1wCR13 (t ,

(12)

R1

(6)

r3R1(t ,

wiR1 ,in(t )

wiR1 (t ) = wiR1 ,out(t ) , nR1

mem FinR1(t ) = F C1(t ) + F C2(t ) + F cat1(t ) + < R1 memFout (t )

mem mem R1 + < R1 memFout (t )wi ,out (t ), i = 1 ,···, nc

R1 (t ) = −n R1F R1(t )(wiR1 , j (t ) − wi , j − 1(t ))

(10)

where A is the cross section area, sR1 i,k is the stoichiometric coefficient of component i in reaction k, and rR1 k is the rate of reaction k. Using backward finite differences to approximate the spatial derivative leads to the well-known nR1 CSTRs in series approximation D

FCr1 out


Article

Industrial & Engineering Chemistry Research dNnCr1 Cr1

inclusion of the impurities is not kinetically limited. The constant distribution coefficient is DCCr1 s l l s wICr1, (t )wCCr1, (t ) = DC Cr1wICr1, (t )wCCr1, (t ) 1 3 1 3

(24)

s l l s wCCr1, (t )wCCr1, (t ) = DC Cr1wCCr1, (t )wCCr1, (t ) 1 3 1 3

(25)

s wCCr1, (t ) 3

(26)

+

s wICr1, (t ) 1

+

s wCCr1, (t ) 1

=1

wiCr1, s(t ) = 0, ∀ i ∉ {C1, C3, I1}

dt

Cr1 Cr1 + Q inCr1(t )nnCr1 Cr1 (t ) − Q out (t )n Cr1 (t ) ,in n

μ3Cr1(t )

(27)

∂(Vn) ∂(Vn) (t , z) + (G(t ) + D(t )) (t , z ) ∂z ∂t (28a)

(Vn)(t , 0) = B(t )V (t )

(28b)

(Vn)(t , +∞) = 0

(28c)

(Vn)(0, z) = n0(z)

(28d)

=

∫0

nCr1

+∞

3

L n(L , t )dL ≈

∑ (njCr1(t )Lj3ΔL) j=1

(35)

Lj = (j − 0.5)ΔLCr1 , j = 1 ,···, nCr1

(36)

Cr1

GCr1(t ) = max(kGCr1(t )sCr1(t )|sCr1(t )|nG

Cr1

BCr1(t ) = GCr1(t )

(29)

where the saturation mass fraction, is assumed to be dependent on the antisolvent concentration as follows (in both Cr1 and Cr2)

dN jCr1 dt

kCr1 B

−

j = 2, ..., n

−1

(37b)

DCr1(t ) Cr1 n1 (t ) 2ΔL

(38a) (38b)

nCr1 D ,

nCr1 G ,

nCr1 B

where is a constant, and and are empirical rate orders. The rate coefficient of dissolution is assumed to be 10 times higher than the growth rate coefficient, eq 38b. It is only used when the supersaturation s is negative to calculate the dissolution rate, DCr1(t), by the nonsmooth eq 37b (if s is positive, DCr1(t) will be zero, and practically is omitted from the population balance). In a similar fashion, the growth rate is only considered when s is positive, eq 37a. The absolute function is required to avoid raising a negative number to a positive, noninteger power. Cr1 FCr1 out is controlled by a valve position, u , using the following equation

(31)

Cr1 Fout (t ) = uCr1(t )CvCr1

Cr1 Cr1 Cr1 + DCr1(t )(N jCr1 + 1(t ) − N j (t ))) + Q in (t )nj ,in (t ) Cr1

+

kDCr1 = 10kGCr1

(30)

(t ) = −ΔL−1(GCr1(t )(N jCr1(t ) − N jCr1 − 1(t ))

Cr1 Q out (t )njCr1(t ),

− nGCr1

)

Cr1 Cr1 kGCr1 = kGCr1 + 273)) ,0 exp( − kG ,1 /(T

There are two potential dynamic regimes: formation-growth of crystals and dissolution-death of crystals. Both should be taken into account when modeling the dynamics of crystallizers. Changes in operating conditions, such as feed flow rate or composition, may lead to switching between these regimes. We use an upwind finite volume method to discretize eq 28, leading to the following DAE set for Cr1 N jCr1(t ) = V Cr1(t )njCr1(t ), j = 1 ,···, nCr1

kGCr1(t )

Cr1

|sCr1(t )|nB

−1

(37a)

Here, min, max, and the absolute function (|·|) are all examples Cr1 of nonsmooth functions, and kCr1 G and kD are the temperaturedependent rate constants calculated by the empirical expressions

wsat C3 ,

Cr2, l wCsat3 (t ) = 0.0318exp( −2.422wAS (t ))

kBCr1(t )

, 0)

(37c)

wCl 3(t ) − wCsat3 (t ) wCsat3 (t )

−1

DCr1(t ) = min(0, kDCr1(t )sCr1(t )|sCr1(t )|nD

where Qin(t) = Fin(t)/ρl and Qout(t) = Fout(t)/ρm(t) are the volumetric flow rates in and out, respectively, and n is the population density of the crystals. Here, G, D, and B are the growth, dissolution, and nucleation (birth) rates of crystals. These rates depend on the supersaturation, s, defined by eq 29 s (t ) =

(34)

Here, it is assumed that no crystals are present in the inlet to Cr1 the first crystallizer; thus, nCr1 is the j,in (t) = 0 for all j. Also, nj discretized population density of crystals of size L ∈ [jΔLCr1, (j + 1)ΔLCr1], j = 0, ..., nCr1, nCr1 is the total number of such discrete size intervals, and ΔLCr1 is a constant. The boundary condition z → + ∞, eq 28c, is replaced with the assumption that NnCr1 Cr1 +1 = 0 at all times. This is accomplished by choosing ΔLnCr1 which is large enough. The following set of nonsmooth algebraic equations are necessary to explicitly distinguish between the growth and dissolution regimes, taking the correct upwind direction into account

The crystals mass at any given time is calculated from a dynamic population balance, assuming that the growth rate is independent of size (McCabe law). The general population balance takes the following form

= Q in(t )nin(t , z) − Q out(t )n(t , z)

(t ) = −ΔL−1( −GCr1(t )NnCr1 (t ) − DCr1(t )NnCr1 Cr1 Cr1(t )) −1

V Cr1(t ) − Vmin |V Cr1(t ) − Vmin| + ϵ

(39)

It is used to determine the optimal dynamic residence time. It is also important to allow draining the vessel during shutdown and filling it during start-up from empty. Here, CCr1 v is the flow coefficient of the valve, VCr1 is the liquid volume in the vessel, Vmin is a small liquid volume below which flow is not permitted, and ϵ is a small positive number making the function Lipschitz continuous.

(32)

dN1Cr1 (t ) = −ΔL−1(GCr1(t )N1Cr1(t ) + DCr1(t )N2Cr1(t )) dt Cr1 Cr1 + BCr1(t )V Cr1(t ) + Q inCr1(t )n1,in (t ) − Q out (t )n1Cr1(t )

(33) E


Article


Figure 2. Effect of recycle to Cr2 shows increasing yield at the expense of impurity (I1) buildup: (a) no recycle, (b) recycle ratio of 0.5, and (c) recycle ratio of 0.9. per wAS (t ) + wSper2 (t ) = 1 2

Cr2 is modeled using the same form of equations as above. The antisolvent flow rate, FAS, is set proportional to the inlet flow rate to Cr16 F AS(t ) = 0.53FinLLE1(t )

(40)

Cr1 mem FinCr2(t ) = Fout (t ) + F AS(t ) + (1 − < R1 mem)Fout (t )

(41)

Here, θmem is the stage cut which is set to a reasonable constant value of 90%. Also, Fper is the permeate mass flow rate (which includes the solvents only), and Rmem is the split fraction determining what fraction of FWF1 goes to the membrane separation and could be controlled by valve u6. Here, we manipulate the split fraction Rmem directly for simplicity. In practice, the stage cut will depend on the membrane area, feed flow rate, and feed composition; however, that is assumed to be a secondary, more detailed, design step which is out of the scope of this paper.

The crystal size distribution in the feed is assumed to be identical to one in Cr1 outlet. The feed is mixed with the recycle stream (in the case that < R1 mem = 0 ), which may contain some AS Cr1 njCr2 ,in (t ) = nj (t ), ∀ j

■

̈ APPROACH) SIMULATION RESULTS (NAIVE In this section, we present simulation results for the two recycle configurations for specific sets of the manipulated variable values (see below), with emphasis on the impurities evolution with time. The purpose of this section is to motivate the use of rigorous optimization techniques, as opposed to relying on simulations of discrete cases, even if these are directed by engineering intuition and judgment. The only parameter changing between the various cases is the recycle ratio, Rmem. In these simulations, we maintain the feed flow rate of the main reactant to the plant (C1) constant at FFeed C1 = 0.2 kg/h, corresponding to the optimal steady-state residence time chosen for the case without recycle in previous studies.6,17 Similarly, in order to focus on the recycle effects only as much as possible, the residence time in the crystallizers, τCrj, is controlled at a set point of 4 h by manipulating the outlet flow rate from each crystallizer using proportional feedback control of the form8

(42)

FinCr2(t )wiCr2 ,in (t ) Cr1 = Fout (t )(εwCr1(t )wiCr1, l(t ) + (1 − εwCr1(t ))wiCr1, s(t )) mem mem Cr2 + (1 − < R1 mem)Fout (t )wi ,out (t ), i = 1 ,···, nc

(43)

All the parameter values for the crystallizers can be found in the Appendix of ref 17. OSN Membrane. Assuming pseudo steady state, the OSN membrane feed and purge mass and component balances are determined by Finmem(t ) = Rmem(t )F WF1(t )

(44)

WF1 wimem (t ) ,in (t ) = wi

(45)

F pur2(t ) = F per(t ) + (1 − Rmem(t ))FWF1(t )

(46)

F pur2(t )wipur2(t ) = F perwiper + (1 − R

mem

(t ))F WF1wiWF1(t ),

i=

1 ,···, ncCr2

eτCrj(t ) = τ Crj(t ) − 4, τ Crj(t )

(47)

Crj = V Crj(t )/(Q out (t ) + ϵ), j ∈ {1, 2, 3, 4}

Assuming the membrane is only permeable to the solvent and antisolvent and that their ratio is maintained in the permeate, the mass and component balances on the membrane are F per(t ) = θ memFinmem(t )

(48)

mem Fout (t ) = Finmem(t ) − F per(t )

(49)

mem mem F per(t )wiper(t ) = Finmem(t )wimem ,in (t ) − Fout wi ,out

(50)

wiper(t ) = 0, i ≠ S1, AS

(51)

per mem wAS (t )/wSper2 (t ) = wAS,in (t )/wSmem 2 1(t ),in

(52)

(53)

(54)

uCrj(t )CvCrj = mid(0.1, 5, 1 + K PCrjeτCrj(t )), j ∈ {1, 2, 3, 4}

(55)

The mid function is used to take the operating limits of the valves into account. In the next section, on optimization, these valve positions (uCrj(t)) will be left as degrees of freedom in the optimization formulation. Recycling to Crystallizer Cr2. It is well known that recycling the mother liquor, after concentration, back to the crystallizers can significantly improve the yield but also the concentration of impurities. This is demonstrated in Figure 2. F


Article


Figure 3. Effect of recycle ratio when recycling the concentrated mother liquor back to Cr2 on the impurity concentrations in the final product vs the campaign time: (a) Rmem = 0.1, (b) Rmem = 0.5, and (c) Rmem = 0.9. The spec values are represented by the straight dashed black lines.

Figure 4. Effect of recycle to Reactor R1 on R1 yield and flow rate of I1 out of R1: (a) no recycle, (b) recycle ratio of 0.5, and (c) recycle ratio of 0.9.

Figure 5. Effect of recycle ratio when recycling the concentrated mother liquor back to R1 on the impurity concentrations in the final product vs the campaign time: (a) Rmem = 0, (b) Rmem = 0.5, and (c) Rmem = 0.9. Spec values are represented by the straight dashed black lines.

than zero; however, it is not possible to know a priori how and when to manipulate the manipulated variables during the campaign. Recycling to Reactor R1. In this case, the recycle is directed to the first reactor, where unreacted C1 can convert to the intermediate C3. This approach was demonstrated by Benyahia et al.,25 where it was claimed that higher yields are achieved by recycling (purging less C1 and C3) but at the expense of higher impurity levels, such that a minimal purge ratio is required in order to obtain on-spec production. Obviously, increased yield is obtained by recycling the main reactant back to convert further in R1; however, because C3 will be recycled as well, the concentration of impurity I1 may increase (see Reaction Set 1). Figure 4 demonstrates that increasing the recycle ratio leads to increased impurity production. The impurity I1 flow rate more than doubles at steady state, from less than 0.1 g/min

As the recycle ratio increases from 0 to 0.9, the yield of the crystallizer increases from 0.85 to 0.98 but at the expense of an increase in impurity concentration in the crystals from 0.6 to 1.6 wt %. Another thing to note is the slower response of the impurity buildup; at a recycle ratio of 0.9, the impurity concentration reaches steady state after 60 h, while without recycle it only takes about 20 h. This obviously impacts all the downstream units and the final product. The resulting impurity profiles of the final product are depicted in Figure 3. It appears that even a recycle ratio of 0.1 leads to a violation of the spec value for I2, such that none of the product is considered on spec. It is important to remember that we only changed the recycle ratio and maintained the controlled residence times (the set points) the same in all cases simulated. It is plausible, as we demonstrate next, that simultaneously changing the residence times can result in an on-spec product even with a recycle ratio greater G


Article


Figure 6. Effect of recycle ratio on the total and on-spec production vs the campaign time: (a) Rmem = 0, (b) Rmem = 0.5, and (c) Rmem = 0.9.

Figure 7. Effect of recycle ratio on the on-spec production for a 200 h long campaign: (a) recycle to Reactor R1, < R1 mem = 1 and (b) recycle to the crystallizer, < R1 mem = 0 .

course leads to very low overall yields. However, for lower recycle ratios, as long as the impurities are maintained below the spec value, the overall productivity and the on-spec productivity increase compared to the “no recycle” case from ∼65 kg on-spec product with no recycle to ∼80 kg with a recycle ratio of 0.5, Figure 6. Figure 7 summarizes the results of the recycle study in terms of productivity and overall yield (for the exact definition of the overall yield, see eq 57). Clearly, there is a great benefit to recycle the concentrated mother liquor back to R1. The onspec productivity increases from 65 kg at 50% yield with no recycle to 84 kg at 65% yield with a recycle ratio of 0.7. Nonetheless, there is still some 6 kg of wasted off-spec product at Rmem = 0.7. Furthermore, increasing the recycle ratio to 0.9 results in a significant off-spec product, such that the yield drops to ∼10%. In the case of recycling the concentrated mother liquor back to the crystallizer, the results are much less promising, and even with a small recycle ratio of 0.1 (remember, all other manipulated variables are kept unchanged), all the product becomes off-spec. Running this kind of study, by simulating various cases (similar to what was done in ref 25), can help identify better solutions, but it is practically impossible to find an optimal solution or to know what better performance can be achieved by dynamically varying the variables one has control over. It is also practically impossible, only by simulations, to determine exactly when is the optimal time to switch from one mode of operation to another (e.g., start the “shutdown” phase, etc.). A more practical way to find superior performance is to solve a

with no recycle to more than 0.2 g/min when setting a recycle ratio of 0.9. The steady-state yield decreases only slightly with increasing recycle ratio (the residence time shifts from its optimal value) because the impurity concentration is still small, on the order of 1 wt %. Nonetheless, given the low spec values of impurities in the final product, this could be a significant increase, as we show below. Figures 5 and 6 show the impact of recycle on the impurity concentrations in the final product. With no recycle, the steady-state concentration of impurity I2 is just below the spec level of 0.3 wt %. There is a transient period where both specs are violated, up to ∼50 h (Figure 5(a)), which reduces the onspec productivity to ∼70 kg from a total productivity of ∼80 kg (Figure 6(a)). Increasing the recycle ratio has a positive effect on the concentration of I2 (reduces), but a negative effect on the concentration of I1 in the product, which seems to increase with time, reaching a higher level of 0.3 wt % after ∼150 h when a recycle ratio of 0.9 is applied. The origin of the higher I1 concentration is clear from the results of Figure 4. The longer transient period and the relative slow increase in I1 concentration is attributed to the longer response time of this recycle loop. The steady-state I2 concentration is almost not effected by the recycle, although it slightly decreases. This is because I2 originates from both C1 and C3 (through the intermediate C4) in the second reactor (see Reaction Set 2). We see in Figures 5(c) and 6(c) that for a high recycle ratio of 0.9 at some point during the campaign the impurities concentration in the final product (I1 + I2) exceeds the spec value, such that the productivity is no longer on-spec, which of H


Article

Industrial & Engineering Chemistry Research dynamic optimization problem, which is what we demonstrate next.

Note that on-spec product accumulation is considered only during the on-spec production epoch, whereas the consumption of C1 is considered over the full campaign duration. The quality constraints, q, defining an on-spec product are

■

DYNAMIC OPTIMIZATION FORMULATION In this section, we define an optimization problem to maximize the overall yield during the campaign, Y, given a target productivity level, P, for a predetermined horizon, [0, tf ]. Here, P is the overall accumulated on-spec product. The full campaign time horizon is split into three epochs: off-spec, on-spec, and off-spec again. These correspond to the intervals [0, ton), [ton, toff], and (toff, tf ], respectively. Production is enforced to be on-spec only during the second epoch by appropriate quality constraints, q. This results in an a priori determined sequence of modes, ensuring computable sensitivities of the resulting hybrid (discrete/continuous) and nonsmooth system.4,35,36 The start and end of the on-spec epoch, ton and toff, are included in the set of decision variables, enabling the algorithm to converge to the longest on-spec production period possible. The main difference between this approach and the classical start-up, steady-state (on-spec) operation and shutdown sequence is that the on-spec interval does not have to be steady. Hence, our formulation is superior to the traditional direct steady-state optimization method because it is a relaxation of the steady-state constraint. The optimization formulation takes the form

(56a)

q(p, t ) ≤ 0, ∀ t ∈ [ton , toff ]

(56b)

Pr P − Mon − spec(p , t f ) ≤ 0

(56c)

pl ≤ p ≤ p u

(56d)

(63)

The above quality constraints are accounted for by introducing new auxiliary differential variables, xα, which transform the path constraints to end-point equality constraints. This is accomplished by the following hybrid, nonsmooth, formulation l 0, t ∈ [0, ton) o o o o o dxα o max(0, q (p, t )), t ∈ [t , t ] , (p, t ) = m on off α ∈ {a , b , c , d} α o o dt o o o o 0, t ∈ (toff , t f ] o n xα(p, 0) = 0

where qα is the original path constraint. The new constraints are thus xα(p, t f ) = 0, α ∈ {a , b , c , d}

Pr 0.34Mon − spec(p , t )MWC1

M C1(p, t )MWAPI

, ∀ t ∈ (0, t f ]

ui(p, t ) = p(i − 1)n + k , ∀ t ∈ [tk − 1 , tk), k = 1 ,···, nt , i = 1 ,···, nu t

(66)

where nu is the number of control functions. We also optimize the durations of the subepochs, defined as the set ; : ; ≔ {τk : τk = tk − tk − 1 , k = 1 ,···, nt , }

(57)

C1

where τk ≥ 0 and ∑kn=t 1τk = tf. ; is added to the decision variables space: pk×(nu+1) = τk, for k = 1, ..., nt. The decision variables space (7 ) has a dimension of nt × (nu+1). We optimize the subepoch durations by introducing a timescale transformation of eq 3a, using the definition t ≔ τk · t*. Now, the integration is run up to predefined integer times, Formulation 674

where M is the total mass of the main reactant, C1, fed to the process, calculated by dM C1 (p, t ) = F C1(p, t ), dt M C1(p, 0) = 0

(58)

and is the total mass of on-spec product produced during the campaign, calculated by

dx (p, t *) = τk f(x(u , t *), y(u , t *), u(p, t *)), dt *

l 0, t ∈ [0, ton) o o o o o o Pr F (p, t ), t ∈ [ton , toff ] (p , t ) = m o o o o o o 0, t ∈ (toff , t f ] o n

Pr Mon ‐spec(p , 0) = 0

(65)

We discretize the entire time horizon into nt subepochs. Each of the main three epochs is divided into a predefined number of subepochs, n1, n2, and n3, respectively, such that ∑i3= 1ni = nt. The control functions, u, are defined as piecewise constant functions, obtaining a constant value at each subepoch

where p ∈  is the decision variables vector, which includes ton and toff. Here, pl and pu are the lower and upper bounds of these variables, respectively. The overall yield assumes the API mass fraction in the product is at its minimum value of 0.34 wt % during the on-spec production mode and is calculated according to eq 57

dt

(62)

qd(p, t ) = wSPr1 (p, t ) − 0.005 ≤ 0, ∀ t ∈ [ton , toff ]

np

Pr dMon ‐spec

(61)

(64)

model 3, ∀ t ∈ [0, t f ]

MPr on‑spec

qb(p, t ) = wIPr2 (p, t ) − 0.003 ≤ 0, ∀ t ∈ [ton , toff ]

∀ t ∈ [ton , toff ]

p∈7

Y (p , t ) =

(60)

qc(p, t ) = wIPr (p, t ) + wIPr2 (p, t ) − 0.005 ≤ 0, 1

max Y (p, t f ) s.t.

Pr qa(p, t ) = 0.34 − wAPI (p, t ) ≤ 0, ∀ t ∈ [ton , toff ]

(59)

∀ t * ∈ [k − 1, k), k = 1 ,···, nt

(67a)

0 = g(x(u , t *), y(u , t *), u(p, t *))

(67b)

x(u , 0) = x 0

(67c)

The original formulation, (56), thus translates to I


Article


Figure 8. Optimal solution of the case with no recycle, Rmem = 0: (a) optimal control profiles and (b) impurity levels in the final product. l u bounded between p5n = 10−6 and p5n = 0.99. The last t+k t+k decision variables are the subepoch durations, τk, with lower and upper bounds plnu × nt+k = 0 h and punu × nt+k = tf h, respectively. For all of the above, k = 1, ..., nt. The initial conditions (eq 3c) in all the cases we simulate represent empty vessels: Mi(0) = 0 or, in some cases, Mi(0) = ϵ to avoid numerical difficulties in finding the consistent initial condition. The reactors, R1 and R2, are assumed to be initially full with the main solvent (C2 in R1 and S1 in R2): wi(0) = 0, ∀i ≠ k, and wk(0) = 1, where wi is the mass fraction of species i, and k is the index of the solvent (most abundant component). In all the crystallizers, it is asserted in addition that no crystals of any size are present at the initial time. To solve the optimization problem, we use IPOPT with a C++ interface.33 We use DAEPACK32 to integrate the system of nonsmooth DAEs and calculate the required sensitivity information (the Jacobian). We set the tolerance in IPOPT to 10−3, and the integration tolerance is set to 10−6 in DAEPACK.

max Y (p, t f , n1 , n2 , n3) p∈7

s.t.

model 67, ∀ t * ∈ [0, nt ]

(68a)

nt

∑ τk − t f

=0

k=1

(68b)

xα(p, nt ) = 0, α ∈ {a , b , c , d}

(68c)

Pr P − Mon ‐spec(p , nt ) ≤ 0

(68d)

pl ≤ p ≤ p u

(68e)

which is the actual form implemented in the optimizer. To conclude, instead of treating ton and toff directly as decision variables in Formulation 56, we defined n1, n2, and n3, and included ; in the decision variables set in Formulation 68, which is equivalent. We control the dynamic performance of the manufacturing campaign by manipulating the flow in the key units: reactor R1 and crystallizers Cr1−Cr4. We achieve that through the respective valves, highlighted in Figure 1 and the recycle ratio, Rmem. We directly control the flow rate of C1 in the feed, FC1, eq 69, whereas the flow out of the crystallizers is determined by valve positions, eq 70 (see ref 5). The recycle ratio is also controlled directly, eq 71.

■

DYNAMIC OPTIMIZATION RESULTS We start by presenting the optimal results for a campaign operated without recycle. This case was extensively studied by Patrascu and Barton,5,17 but due to the minor changes in crystallizer Cr1 and crystallizer Cr2 models, and for the sake of completeness and proper comparison, we update the optimal solution for the current model. We present these results here for comparison with the recycling cases. In the no recycle case, we set n1 = 1, n2 = 2, and n3 = 1 so that nt = 4. Therefore, the first and last off-spec epochs are discretized into just one subepoch, while the on-spec epoch is split into two subepochs. Thus, we have a total of 24 decision variables (five control functions discretized into four subepochs, plus the four subepochs duration). In the simulations above, with no recycle, we obtained ∼80 kg total productivity with ∼65 kg being on-spec at a yield of ∼50%. Here, we set the minimal on-spec productivity to be 90 kg; i.e., in Formulation 68 P = 90 kg, and use the same campaign time, tf, of 200 h. Let us examine the optimal results. Figure 8 shows the control profiles and the final product impurity profiles. The objective function (the yield) and the productivity are presented in Figure 9. The campaign is characterized by

u1(p, t ) ≔ F C1(p, t ) = pk , ∀ t ∈ [tk − 1 , tk), k = 1 ,···, nt (69)

uj + 1(p, t ) ≔ uCrj(p, t ) = pj × n + k , ∀ t ∈ [tk − 1 , tk), k = 1 ,···, nt , j = 1 ,···, 4 t

(70)

u6(p, t ) ≔ Rmem(p, t ) = p5n + k , ∀ t ∈ [tk − 1 , tk), k = 1 ,···, nt t

(71)

The lower and upper bounds on the C1 feed flow rate are plk = 0.01 and puk = 0.5 kg/h, respectively. The outlet valve positions of the crystallizers are bounded between pjl × nt+k = 0.001 and puj × nt+k = 1, for j = 1, ..., 4. The recycle ratio is J


Article


Figure 9. Overall campaign performance for the optimal solution in the case with no recycle, Rmem = 0. The on-spec productivity equals the total productivity (90 kg), and the final yield is 67%.

Figure 10. Recycle back to Cr2, < R1 mem = 0 . Solution of the optimization problem 68 with P = 90 kg and campaign time of 200 h: (a) optimal control profiles and (b) impurity levels in the final product. The short and small deviations from the spec values observed at ∼125 h are within the numerical tolerance of the constraints.

before toff, i.e., still during the on-spec epoch. During this phase, the vessels are drained, and the C1 feed flow rate is turned down to its lower limit, but in a controlled way, such that the product is still on-spec. We want to emphasize that this result is an outcome of optimizing the yield; i.e., we did not enforce any shutdown procedure, e.g., draining the vessels, etc. Recycling to Cr2. In this case, we set n1 = 1, n2 = 3, and n3 = 1, such that, nt = 5. We also add one more dynamic control variable, the recycle ratio Rmem; thus, the optimization problem has a total of 35 decision variables. The optimal control profiles and the impurity levels in the final product are shown in Figure 10. Although we chose three subepochs for the onspec epoch, we see from the control profiles that their value during the first and second subepochs (τ2 and τ3) are very similar and only significantly change during the last on-spec subepoch (τ4), similar to the optimal solution to the previous “no recycle” case. This is explained by the impurity profiles obtained, Figure 10(b). One main difference is that the spike in the impurity concentration observed during the “start-up” phase has seemed to disappear. More importantly, the impurity levels are pushed all the way to both impurity specs for most of

three apparent phases: start-up, steady-state-like production, and shutdown. Interestingly, the optimal solution exhibits onspec productivity that equals the total productivity; that is, all the product is on-spec such that there is no wasted product. Also, it is apparent that the productivity obtained equals the minimum required, i.e., 90 kg, so constraint 68d is active. The overall yield is ∼67%, which is a significant improvement over ̈ the base case simulated in the Simulation Results (Naive Approach) section (Figures 5,6, and 7), where the on-spec productivity was 65 kg at a yield of 50%. The start-up phase is required to continue until a product with the minimal concentration of API is produced, resulting in ton = 13.7 h. The values of the controls are determined by the optimizer such that the peak in the impurity levels does not exceed the spec values, which is why the feed flow rate of C1 is lower during this phase. Examining Figure 8(b), the steadystate like phase is very obvious, where the impurity level of I2 is pushed all the way to the spec value and maintains almost steady-state dynamics for most of the campaign duration. The last epoch’s duration (τ4) is actually zero because the best solution is still on-spec during the “shutdown” phase, which in this case happens during the third subepoch, τ3, which is still K


Article


Figure 11. Recycle back to R1, < R1 mem = 1. Solution of the optimization problem 68 with P = 95 kg and campaign time of 200 h: (a) optimal control profiles and (b) impurity levels in the final product.

Figure 12. Overall transient campaign performance for the optimal solutions of the two recycle configuration: yield (a) and on-spec productivity (b). When recycling to Cr2, the on-spec productivity equals the total productivity (90 kg), and the final yield is 72%. When recycling to R1, the optimal solution results in on-spec productivity (equals the total productivity) of 102 kg, and the final yield is 82%.

to 0.2 kg/min for most of the campaign duration, as was set in the previous Simulation Results (Naive Approach) section. The only deviation from this value is during the start-up and shutdown phases. This is not surprising, because based on the reaction scheme and the kinetic parameters, this value is known to give the optimal yield in R1.12 The optimal recycle ratio for most of the campaign duration is around 0.75, as can be observed in Figure 10(a). This is a surprising result when one bears in mind Figure 7(b), where it appears that even a recycle ratio of 0.1 led to off-spec product. The high recycle ratio is made possible by simultaneously varying the residence times of R1 and the crystallizers, thereby keeping the impurity levels on spec. Recycling to R1. In this case, just as in the previous case, we set the number of subepochs to n1 = 1, n2 = 3, and n3 = 1, such that, nt = 5. Figure 11(a) shows the optimal control profiles, and Figure 11(b) shows the impurity concentrations in the final product for the optimal solution found with recycle back to the reactor R1. As opposed to the previous case, here the manipulated

the campaign duration (I2 is at 0.3 wt %, and I1 + I2 is at 0.5 wt %). The extra degree of freedom which was added to this configuration compared to the no recycle case (Rmem) has enabled the impurity levels to reach both impurity specifications for most of the campaign duration, as opposed to only I2 reaching its spec value in the no recycle case (remember that I2 is only formed downstream of Cr2). The recycle back to the crystallizer has thus increased the separation yield on the expense of increased impurity I1 concentration, which was possible because, as illustrated in Figure 8(b), it was below the spec limit for most of the campaign duration. The fact that the impurity levels are pushed all the way to the spec value for most of the campaign duration is a sign that better yields cannot be achieved through additional number of subepochs. In fact, a similar yield could have been obtained by using n2 = 2. This solution shows an improvement compared to the “no recycle” case in the overall yield, from 67% to 72%, maintaining the same on-spec productivity of 90 kg. Also note that the optimal reactant flow rate to the system (to R1), FC1, was determined by the optimizer to be very close L


Article

Industrial & Engineering Chemistry Research variables (controls) obtain very different values during the subepochs τ2 and τ3. The result is that there is no apparent steady-state during the campaign. That is, the impurity concentrations continue to change constantly during the entire campaign, staying below the spec levels throughout. Dynamically changing the residence times in the main reaction and separation units allows for a very high recycle rate, close to 1, without violating the impurity spec levels in the final product. First, impurity I2 reaches its spec value (at t ≈ 20 h), which has been achieved by a change in the control values at a prior time, around 10 h on stream. Once the combined impurities level reaches its spec level of 0.5 wt %, another change in the control values is triggered, which leads to a reduction in the concentration. The impurity I1 concentration continues to build in the final product, while impurity I2 concentration drops throughout the campaign duration until the final subepoch, in which, similarly to the other optimal solutions, the “shutdown” procedure is initiated, the feed drops to minimum, and the vessels are drained. The fact that this solution does not exhibit constant impurity levels at spec values (for some long period during the campaign, i.e., “the turnpike”) suggests that even a better solution may exist, if we add more subepochs as degrees of freedom. However, the benefit will probably be minor compared to the improvement that is already achieved here. This was demonstrated in ref 5. The campaign performance shows a significant improvement compared to the case without recycle and also a significantly better performance compared to the first recycle configuration, back to Cr2, Figure 12. The overall productivity for the campaign exceeds the target of at least 95 kg, and ∼102 kg of on-spec product is obtained (i.e., Constraint 68d is inactive). The overall yield is 82%. Interestingly, the overall reactant mass fed to the process is about the same in both optimal recycle configurations: 37.7 kg in the optimal case for recycle to Cr2 vs 37.4 kg in the optimal case for recycle to R1. The reason for the significant better performance when recycling to the reactor is that this configuration allows one to increase the I1 concentration while reducing the I2 concentration at very high recycle ratios. A hint for this was observed in the Simulation Results (Naive Approach) section, Figure 5, although in that case the decrease in I2 was relatively small since not all the degrees of freedom were used (the residence times in the crystallizers remained controlled at a constant set point). I2 originates from C1 through Reaction 2b in Reactor R2. By recycling C1 to R1, the overall conversion of C1 in the recycle loop increases (even though the conversion in reactor R1 is decreasing), which leads to a lower flow rate of C1 to R2 compared to the case of recycling to Cr2, Figure 13. This configuration enables high yields and productivities because it purges less C1 and C3, increases the concentration of I1 at the expense of I2 in the final product. This is not possible by recycling only to the crystallizer because by doing so one can only control the amount of impurities (C1 and I1) that are incorporated in the C3 crystals. Once these reach the maximum allowed by the specifications (on the final product) the yield cannot be improved further (Figure 10). The recycle to R1 also enables operating R1 at much shorter residence times (higher feed flow rates). This reduces the yield at short times (see, for example, Figure 12 at 50 h) but leads to lower I1 concentrations (at these short times) and eventually higher overall yield. This behavior is attributed to the dynamic design of the campaign; we can feed more reactant at the early stages

Figure 13. C1 flow rate going to the second reactor (R2) where I2 is formed, according to the optimal solutions of the two recycle configurations.

of the campaign at the expense of a lower temporary yield (Figure 11(a)) but lower the feed rate at a later stage and drain the vessels toward the end of the campaign.

■

CONCLUSIONS In this paper, we used nonsmooth DAEs, a new dynamic modeling paradigm, to simulate the dynamic performance of a continuous pharmaceutical manufacturing plant, which includes recycle. We embedded the dynamic model in an optimization problem and solved for the optimal dynamic performance of the plant in terms of yield and productivity. We have compared the theoretical results of two recycle configurations to a case with no recycle, inspired by an experimental pilot plant operated in the Novartis-MIT Center for Continuous Manufacturing. We illustrated that our formulation is superior to traditional steady-state optimization approaches. This is achieved by allowing an overall transient behavior (as opposed to forcing a true steady state), enforcing the quality constraints during an internal epoch, while optimizing its duration. The flow rate through key process units and the recycle ratio have been controlled assuming piecewise constant functions. Importantly, we optimized the timing of the control switchings by an appropriate time transformation, which allowed using a relatively small number of decision variables. We illustrated how the performance of the once-through operation could be improved by recycling the concentrated mother liquor back to the crystallizer (case 1) to increase the yield compared to a case with no recycle from 67% to 72% at the same productivity of 90 kg per campaign duration of 200 h. Yet recycling to the reactor (case 2) shows even a more significant improvement in terms of both the total on-spec productivity and the overall yield, obtaining a yield of 82% at a productivity of 102 kg per same campaign duration of 200 h. In the optimal solutions, the recycle ratio was kept at a high level, ∼0.75 and ∼0.95 for case 1 and case 2, respectively, for most of the campaign duration. These results could not have been predicted only by running simulations of the system, which showed that on-spec impurity levels are violated even for small recycle ratios, when the other controls are maintained at the value determined for the once-through case. This paper demonstrates the potential of implementing such nonsmooth dynamic formulations in online economic optimization and control systems, known as nonlinear model predictive control. In such systems, an updated optimal control profile will be calculated based on real-time measurements. These systems should also be robust and take parametric M


Article

Industrial & Engineering Chemistry Research uncertainty into account.21,22,26 In actual implementations the operating spec levels should be below the product spec levels to allow some room for the online control system to operate, without violating the product specifications. Same goes for safety specifications (not addressed in this paper).

■

(15) Kumar, A.; Vercruysse, J.; Vanhoorne, V.; Toiviainen, M.; Panouillot, P. E.; Juuti, M.; Vervaet, C.; Remon, J. P.; Gernaey, K. V.; Beer, T. D.; Nopens, I. Conceptual framework for model-based analysis of residence time distribution in twin-screw granulation. Eur. J. Pharm. Sci. 2015, 71, 25−34. (16) Su, Q.; Nagy, Z. K.; Rielly, C. D. Pharmaceutical crystallisation processes from batch to continuous operation using MSMPR stages: Modelling, design, and control. Chem. Eng. Process. 2015, 89, 41−53. (17) Patrascu, M.; Barton, P. I. Optimal campaigns in end-to-end continuous pharmaceuticals manufacturing. Part 1: Nonsmooth dynamic modeling. Chem. Eng. Process. 2018, 125, 298−310. (18) Franceschini, G.; Macchietto, S. Model-based design of experiments for parameter precision: State of the art. Chem. Eng. Sci. 2008, 63, 4846−4872. (19) Yu, L. X. Pharmaceutical quality by design: product and process development, understanding, and control. Pharm. Res. 2008, 25, 781− 91. (20) Kockmann, N.; Gottsponer, M.; Zimmermann, B.; Roberge, D. M. Enabling continuous-flow chemistry in microstructured devices for pharmaceutical and fine-chemical production. Chem. - Eur. J. 2008, 14, 7470−7477. (21) Braatz, R. D.; Nagy, Z. K. Robust nonlinear model predictive control of batch processes. AIChE J. 2003, 49, 1776−1786. (22) Mesbah, A.; Paulson, J. A.; Lakerveld, R.; Braatz, R. D. Model Predictive Control of an Integrated Continuous Pharmaceutical Manufacturing Pilot Plant. Org. Process Res. Dev. 2017, 21, 844−854. (23) Griffin, D. W.; Mellichamp, D. A.; Doherty, M. F. Reducing the mean size of API crystals by continuous manufacturing with product classification and recycle. Chem. Eng. Sci. 2010, 65, 5770−5780. (24) Alvarez, A. J.; Singh, A.; Myerson, A. S. Crystallization of cyclosporine in a multistage continuous MSMPR crystallizer. Cryst. Growth Des. 2011, 11, 4392−4400. (25) 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. (26) Boukouvala, F.; Niotis, V.; Ramachandran, R.; Muzzio, F. J.; Ierapetritou, M. G. An integrated approach for dynamic flowsheet modeling and sensitivity analysis of a continuous tablet manufacturing process. Comput. Chem. Eng. 2012, 42, 30−47. (27) Wang, J.; Lakerveld, R. Continuous Membrane-Assisted Crystallization to Increase the Attainable Product Quality of Pharmaceuticals and Design Space for Operation. Ind. Eng. Chem. Res. 2017, 56, 5705−5714. (28) Keshavarz, L.; Steendam, R. R. E.; Frawley, P. J. Impact of Mother Liquor Recycle on the Impurity Build-Up in Crystallization Processes. Org. Process Res. Dev. 2018, 22, 1541−1547. (29) Stechlinski, P.; Patrascu, M.; Barton, P. I. Nonsmooth differential-algebraic equations in chemical engineering. Comput. Chem. Eng. 2018, 114, 52−68. (30) Stechlinski, P. G.; Barton, P. I. Dependence of solutions of nonsmooth differential-algebraic equations on parameters. J. Differ. Equ. 2017, 262, 2254−2285. (31) RES Group Inc., 2019. http://www.resgroupinc.com/ (accessed June 2019). (32) Tolsma, J.; Barton, P. I. DAEPACK: An Open Modeling Environment for Legacy Models. Ind. Eng. Chem. Res. 2000, 39, 1826−1839. (33) 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. (34) Foley, M. A.; Jamison, T. F. Amide bond formation via reversible, carboxylic acid-promoted lactone aminolysis. Org. Process Res. Dev. 2010, 14, 1177−1181. (35) Galán, S.; Feehery, W. F.; Barton, P. I. Parametric sensitivity functions for hybrid discrete/continuous systems. Appl. Numer. Math. 1999, 31, 17−47. (36) Stechlinski, P. G.; Barton, P. I. Generalized Derivatives of Differential−Algebraic Equations. J. Optim. Theory Appl. 2016, 171, 1.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +1 617 2536526. ORCID

Paul I. Barton: 0000-0003-2895-9443 Present Address ‡

M. Patrascu: ExxonMobil Research and Engineering, Annandale, New Jersey 08801, United States

Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Plumb, K. Continuous processing in the pharmaceutical industry: Changing the mind set. Chem. Eng. Res. Des. 2005, 83, 730−738. (2) Cervera-Padrell, A. E.; Skovby, T.; Kiil, S.; Gani, R.; Gernaey, K. V. Active pharmaceutical ingredient (API) production involving continuous processes - A process system engineering (PSE)-assisted design framework. Eur. J. Pharm. Biopharm. 2012, 82, 437−456. (3) Gutmann, B.; Cantillo, D.; Kappe, C. O. Continuous-flow technology - A tool for the safe manufacturing of active pharmaceutical ingredients. Angew. Chem., Int. Ed. 2015, 54, 6688− 6728. (4) Sahlodin, A. M.; Barton, P. I. Optimal campaign continuous manufacturing. Ind. Eng. Chem. Res. 2015, 54, 11344−11359. (5) Patrascu, M.; Barton, P. I. Optimal campaigns in end-to-end continuous pharmaceuticals manufacturing. Part 2: Dynamic optimization. Chem. Eng. Process. 2018, 125, 124−132. (6) 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. (7) Lakerveld, R.; Benyahia, B.; Braatz, R. D.; Barton, P. I. ModelBased Design of a Plant-Wide 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) Kralj, J. G.; Sahoo, H. R.; Jensen, K. F. Integrated continuous microfluidic liquid-liquid extraction. Lab Chip 2007, 7, 256−263. (10) Quon, J. L.; Zhang, H.; Alvarez, A.; Evans, J.; Myerson, A. S.; Trout, B. L. Continuous crystallization of Aliskiren hemifumarate. Cryst. Growth Des. 2012, 12, 3036−3044. (11) Zhang, H.; Quon, J.; Alvarez, A. J.; Evans, J.; Myerson, A. S.; Trout, B. Development of continuous anti-solvent/cooling crystallization process using cascaded mixed suspension, mixed product removal crystallizers. Org. Process Res. Dev. 2012, 16, 915−924. (12) 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. Org. Process Res. Dev. 2014, 18, 402−409. (13) Wulkow, M.; Gerstlauer, A.; Nieken, U. Chem. Eng. Sci. 2001, 56, 2575−2588. (14) Li, Z.; Kind, M.; Gruenewald, G. Modeling the Growth Kinetics of Fluidized-Bed Spray Granulation. Chem. Eng. Technol. 2011, 34, 1067−1075. N


Optimal Dynamic Continuous Manufacturing of Pharmaceuticals with

Recommend Documents