A modeling strategy for integrated batch process development based on mixed-logic dynamic optimization

Computers and Chemical Engineering 94 (2016) 287–311

Contents lists available at ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

A modeling strategy for integrated batch process development based on mixed-logic dynamic optimization ˜ c , Wolfgang Marquardt b Marta Moreno-Benito a,∗ , Kathrin Frankl b , Antonio Espuna a b c

Centre for Process Systems Engineering, University College London, Torrington Place, London WC1E 7JE, UK Process Systems Engineering, RWTH Aachen, Turmstraße 46, D-52056 Aachen, Germany Department of Chemical Engineering, Universitat Politècnica de Catalunya, ETSEIB, Av. Diagonal 647, E-08028 Barcelona, Spain

a r t i c l e

i n f o

Article history: Received 3 January 2016 Received in revised form 10 June 2016 Accepted 31 July 2016 Available online 6 August 2016 Keywords: Batch process synthesis Plant allocation Dynamic optimization Generalized disjunctive programming Synchronization Multistage modeling

a b s t r a c t This paper introduces an optimization-based approach for the simultaneous solution of batch process synthesis and plant allocation, with decisions like the selection of chemicals, process stages, task-unit assignments, operating modes, and optimal control profiles, among others. The modeling strategy is based on the representation of structural alternatives in a state-equipment network (SEN) and its formulation as a mixed-logic dynamic optimization (MLDO) problem. Particularly, the disjunctive multistage modeling strategy by Oldenburg and Marquardt (2008) is extended to combine and organize single-stage and multistage models for representing the sequence of continuous and batch units in each structural alternative and for synchronizing dynamic profiles in input and output operations with material transference. Two numerical examples illustrate the application of the proposed methodology, showing the enhancement of the adaptability potential of batch plants and the improvement of global process performance thanks to the quantification of interactions between process synthesis and plant allocation decisions. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction Specialty chemical industries are meant to assimilate and execute many different processes along their lifecycle, as many specialty goods are subject to quickly changing markets or may be displaced by rapidly advancing technologies. The fast introduction of new products and processes into manufacturing facilities is crucial to be first in the marketplace. However, such processes should be also sustainable to ensure a long-term economic viability and maintain their competitiveness in front of stockholders. The investment of time and resources to design, construct, and validate a new plant each time that a different product is approved for commercialization is not a viable strategy in most cases. This way, batch systems are highly appropriate for specialty chemical sector because their inherent flexibility allows their adaptation and reconfiguration to incorporate new production lines with none or only partial plant modifications. After the discovery of a new product and acknowledgment of its market opportunity, the integrated batch process development problem addresses the synthesis of processing schemes and the allocation of manufacturing facilities (Rippin, 1983;

∗ Corresponding author. E-mail address: [email protected] (M. Moreno-Benito). http://dx.doi.org/10.1016/j.compchemeng.2016.07.030 0098-1354/© 2016 Elsevier Ltd. All rights reserved.

Stephanopoulos and Reklaitis, 2011). Both sub-problems have several degrees of freedom closely interrelated, which ultimately define the process and the required equipment investment in new or previously existing manufacturing facilities in grassroots or retrofit scenarios, respectively. For instance, the development of a batch process stage includes decisions like the operating mode for determining the process structure – i.e. parallel, series or single unit operation – that best suits the process requirements in the batch system. This decision is closely related to the task-unit assignment and equipment selection, especially in the case of existing plants with a given set of units and connecting pipelines, and is crucial to exploit the adaptability potential of batch systems. Simultaneously, each equipment unit has to accommodate the assigned tasks using a batch or semi-continuous procedure, which involves either a chain of batch operations or an intermittent use of continuously operated plant elements, respectively. Notably, batch operations are defined by their duration and dynamic feed-forward trajectories of control variables, and their feasibility depends on physical constraints like the capacity of the selected units. Additionally, it is necessary to determine the batch size and the number of batches, as they are responsible for the total demand satisfaction and interrelate with the equipment capacity and process conversion. Other degrees of freedom for full process development include the selection of process stages, the recirculation of intermediate flows, the installation of buffer tanks for working with different cycle times or

288

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

Notation Acronyms AIBN azobisisobutyronitrile AN acrylonitrile branch-and-bound B&B CHR convex-hull relaxation CNF conjunctive normal form dimethylformamide DMF DO dynamic optimization end-point constraints EPC GBD generalized Benders decomposition generalized disjunctive programming GDP IFS initial feasible solutions KPI key performance indicators mixed-integer dynamic optimization MIDO MILP mixed-integer linear programming MINLP mixed-integer non-linear programming mixed-logic dynamic optimization MLDO mathematical programming MP NCO necessary conditions of optimality NLP non-linear programming OA outer approximation PC path constraints piecewise constant PWC RTN resource-task network SEN state-equipment network STN state-task network VA vinyl acetate General variables ˙ z(t) derivatives of differential process variables z(t) differential process variables y(t) algebraic process variables udyn (t) dynamic control variables ustat time-invariant or static continuous control variables uint integer control variables logical or Booleans decisions uBool ubin binary decisions  algebraic time-invariant variables p process parameters General functions f DAE system path constraints g ge end-point constraints algebraic equations evaluated at the final time h l set of relations that define initial conditions m stage-to-stage continuity between consecutive mathematical stages DAE system in disjunctive equations fd gd path constraints in disjunctive equations gd,e end-point constraints in disjunctive equations hd algebraic equations evaluated at the final time in disjunctive equations ld set of relations that define initial conditions in disjunctive equations md stage-to-stage continuity between consecutive mathematical stages in disjunctive equations Bd equations system to define bypass stages in disjunctive equations  objective function

Problem sets and indices  set of technological specifications indexed by  j ⊆  technological specification of unit j ∈ U, |j | = 1, ∀ j ∈ U indexed by  i set of operating modes in task i ∈ PS indexed by C set of chemical compounds involved into the process indexed by c Cjs ⊆ C subset of potential reactants, solvents, or catalysts in unit j ∈ U indexed by c Di, ⊆ Ui subset of batch units Ui in task i ∈ PS whose input flow rate is a control variable in configuration ∈ i . It is defined such that |Di, | = DOFi, − |Ui |, where |Ui | represents the degrees of freedom removed by output flow rates indexed by j Ij ⊆ Kj subset of input stages for unit j ∈ U at Level 1 indexed by k set of all existing and potential equipment pieces, J J = U ∪ T ∪ Sp ∪ Mx indexed by j subset of equipment pieces within potential task i ∈ Ji ⊆ J PS indexed by j set of stages for unit j ∈ U at Level 1 indexed by k Kj L set of potential stages at Level 0, L = {1, . . ., Lmax } indexed by k and l La ⊆ L subset of active stages at Level 0 indexed by l subset of stages at Level 0 where unit j ∈ U can start Lj0 ⊆ L

its operation, Lj0 = {1, |Kj | − |Oj | + 1 | j = / j, j ∈ U} indexed by l set of pipelines for unit j ∈ U at Level 1 indexed by Mj m Mjin ⊆ Mj subset of input pipelines to unit j ∈ U at Level 1 indexed by m Mjout ⊆ Mj subset of output pipelines from unit j ∈ U at Level 1 indexed by m Mx set of existing and potential mixers indexed by j N set of pipelines at Level 0 indexed by n Ni ⊆ N subset of pipelines at Level 0 for task i ∈ PS indexed by n Niin ⊆ Ni subset of input pipelines to process stage i ∈ PS indexed by n Niout ⊆ Ni subset of output pipelines to process stage i ∈ PS indexed by n Ni,0 ⊆ N subset of pipelines at Level 0 whose flow rate is

f

Ni ⊆ N

restricted to zero in configuration ∈ i of task i ∈ PS indexed by n subset of pipelines at Level 0 for task i ∈ PS whose f

flow rate is fixed by the preceding task, |Ni | = 1 except for process stages preceded by a buffer j ∈ T indexed by n Nib ⊆ Ni bypass pipeline for process stage i ∈ PS, |Nib | = 1 indexed by n and b Nr ⊆ N subset of pipelines at Level 0 for recirculation indexed by n Njin ⊂ N subset of input pipelines to unit j ∈ J at Level 0 indexed by n Njout ⊂ N subset of output pipelines from unit j ∈ J at Level 0 indexed by n Oj ⊆ Kj subset of output stages for unit j ∈ U at Level 1 indexed by k P⊂C subset of desired products indexed by p

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

PS Q Sp T Tnr ⊆ T U Ui ⊆ U

set of potential process stages or tasks indexed by i set of ordered positions that can be assumed by unit procedures of j ∈ U indexed by q set of existing and potential splitters indexed by j set of existing and potential storage tanks indexed by j buffer tank for potential recirculation of flow n ∈ N r , |Tnr | = 1, ∀n ∈ N r indexed by j set of existing and potential batch units indexed by j subset of existing and potential batch units in process stage i ∈ PS indexed by j

Problem parameters Demp demand of product p ∈ P degrees of freedom with regard to the flow rates DOFi, at Level 0 at process stage i ∈ PS, according to each configuration ∈ i 0 lj,q starting stage of unit j ∈ U when the task-unit 0 = qth element of the ascendBoolean Wj,q is true, lj,q

Lmax pj,c

ing sort of Lj0 of unit j ∈ U maximum number of stages at Level 0 values for the set of process parameters pj in unit j ∈ U when potential chemical alternative c ∈ Cjs is selected

Problem variables (* Control variables) Batchp production size associated with each batch of product p ∈ P j Fm,k (t)(*) flow rate of input or output connection m ∈ Mjin ∪ Mjout in stage k ∈ Kj of unit j ∈ U at Level 1 flow rate of connection n ∈ N in stage l ∈ L at Level Fn,l (t) 0 j int k (t)(*) internal control variable in stage k ∈ Kj of unit j ∈ U at Level 1 and stage k ∈ L of unit j ∈ J\U at Level 0 NBp (*) number of batches of product p ∈ P Rn (*) recirculation Boolean for intermediate flow recycle n ∈ Nr j Sc (*) compound Boolean for reactant, solvent, or catalyst c ∈ Cjs in unit j ∈ U Shortfallp unaccomplished demand of product p ∈ P Sizej (*) discrete capacity of unit j ∈ U ∪ T final time at Level 0 tend Tf total time at Level 0 tj,end final time of unit j ∈ U at Level 1 Tj,f total time of unit j ∈ U model at Level 1 duration of stage l ∈ L at Level 0 tl (*) tj,s starting time of unit j ∈ U at Level 1 j duration of stage k ∈ Kj of unit j ∈ U at Level 1 tk starting time at Level 0 ts jk (t) material volume in stage k ∈ Kj of unit j ∈ U t Level 1 and stage k ∈ L of unit j ∈ J\U at Level 0 j V (*) technology Boolean for specification  ∈ j in processing unit j ∈ U Wj,q (*) assignment Boolean for processing order q ∈ Q in unit j ∈ U X i (*) configuration Boolean for alternative ∈ i in process stage i ∈ PS j xc,m,k (t) composition of compound c ∈ C in input or output connection m ∈ Mjin ∪ Mjout and stage k ∈ Kj of unit j ∈ U at Level 1

xc,n,l (t) Yj (*) Zi (*)

289

composition of compound c ∈ C in connection n ∈ N and stage l ∈ L at Level 0 equipment Boolean for processing or storage unit j ∈ U∪T process stage Boolean for task i ∈ PS

batch sizes in successive process stages, the selection from a set of alternative technologies, and the selection of chemicals. The mathematical formulation and the solution procedure for addressing this problem are challenging and computationally demanding activities for several reasons. The batch nature of the process involves not only a combinatorial assessment for linking the equipment and task networks (Gani and Papaeconomou, 2005) but also the adjustment of processing conditions along time in each process stage. This way, dynamic models are required for representing the process performance and dynamic profiles of control variables are necessary to fully explore the attainable region of the problem (Srinivasan et al., 2003). Additionally, the transition of batch operations in batch process stages entails the mathematical representation of discrete events (Barton et al., 1998). Besides, qualitative information has to be included regarding decisions like tasks selection, sequencing and splitting, equipment assignment, and chemicals selection. Finally, batch integrity has to be guaranteed for every structural processing alternative by synchronizing the material transference between batch and semi-continuous plant elements as a function of the selected processing scheme and the task performance achieved. In batch systems, the synchronization of the material transference between consecutive processing units requires establishing a link between the dynamic profiles of the input and output flows and the duration of the operations with material transference that occur at the same time. Given the complexity described above, the problem of batch process development is usually decomposed in sub-problems and solved independently. Most of the literature follows the natural decomposition into process synthesis and plant allocation sequence for simplifying the problem. This way, some studies address the synthesis of batch processes by disregarding the subsequent assignment of equipment items (Ahmad, 1997; Ali, 1999; Barton et al., 1999; Linninger et al., 1995; Papaeconomou, 2005; Sharif et al., 2000), while some contributions include allocation decisions through sequential decision-making procedures (Cavin et al., 2004; Iribarren, 1985; Iribarren et al., 1994; Mosat et al., 2008). Other works focus on the optimal operation of individual units through the application of tools like dynamic optimization (DO). DO or optimal control is a well-received research field in the last decades that addresses operating decisions like the dynamic profiles of control variables (Binder et al., 2001; Srinivasan et al., 2003), duration of batch operations (Vassiliadis et al., 1994), and structural decisions that determine the configuration of processing units (Jain et al., 2013; Oldenburg et al., 2003). The design of batch plants is a further area of extensive research (Reklaitis, 1989; Barbosa-Póvoa, 2007) where the allocation of manufacturing facilities and plant sizing are usually solved using predefined recipes or approximated models – e.g. the fixed time and size factor models – to simplify the process stage representation (Robinson and Loonkar, 1972). In this case, laboratory procedures used during the product development phase are scaled-up to pre-specify the information related to the batch process synthesis, like the chemicals, operating conditions or processing times involved in the process. However, divide and conquer strategies jeopardize the efficiency of the resulting process for several reasons. First, the sequential solution of batch process development involves losing a significant part of the interaction among the degrees of freedom, thus

290

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

restricting the evaluation of the system trade-offs according to industrial-scale objectives (Moreno-Benito et al., 2014). The potential of decisions like temperature profiles or task durations for improving global processing times, processing conditions, and tasks sequences is only relevant if production-scale restrictions and objectives are taken into account in the decision-making. Besides, the targets in bench-scale experiments – e.g. reaction yield – frequently differ from those of industrial-scale manufacture, which are usually global objectives that include all the different parts of the process – e.g. the total cost and environmental impact. Additionally, the direct implementation of fixed recipes in existing equipment may force its operation in extreme and even infeasible conditions (Grossmann et al., 1987). This way, it is crucial to take into account the physical restrictions of the given plant to ensure the feasible execution of the new processes. Overall, to fully optimize the processes and guarantee the best utilization of available equipment and new investments, process synthesis and plant allocation should be solved simultaneously. Significant efforts have been made for combining synthesis and allocation sub-problems with different degrees of integration over the last fifteen years. For instance, some studies address this problem using DO and mixed-integer dynamic optimization (MIDO) formulations to optimize the reference trajectories of control variables (Charalambides et al., 1995; Sharif et al., 1999), although they involve the pre-specification of the sequence of process stages and allocation decisions like the selection of the operating mode and the task-unit assignment. Other approaches simplify the problem through the approximation of the batch process behavior using algebraic models – e.g. the named screening models (Allgor et al., 1999) and posynomial functions (Iribarren et al., 2004). These approaches have been complemented with the optimization of the dynamic feed-forward trajectories of control variables (Allgor and Barton, 1999), albeit using an iterative evaluation of structural and performance decisions instead of a simultaneous solution procedure. In contrast, other works opt for constant set-points for defining the processing conditions while solving the synthesis and allocation problems simultaneously (Chakraborty and Linninger, 2002; Linninger and Chakraborty, 1999). Additionally, several authors seek the sustainable design of processing schemes in retrofit scenarios taking into account the previously existing equipment, while they use heuristic approaches based on path flow decomposition and indicator-based assessment (Banimostafa et al., 2012; Bumann et al., 2011; Carvalho et al., 2009; Halim et al., 2011; Simon et al., 2008). Besides, the representation of the process performance using dynamic models in allocation problems has received a strong impulse in the last two decades in the context of short-term scheduling (Bhatia and Biegler, 1996; Capón-García et al., 2013; Chu and You, 2014; Frankl et al., 2012; Nie et al., 2012). However, these contributions do not cover other synthesis and design decisions such as the selection of process technologies, chemicals, or equipment sizing, among others. Finally, the synchronization of process stages for linking operations with material transference in consecutive processing units is addressed in the context of scheduling and design of batch plants (Ferrer-Nadal et al., 2008; Furman et al., 2007), bearing in mind that the sequence of processing units is not known beforehand. The principal objective of the synchronization in these works is to elude physically inconsistent allocation solutions; however, material and energy balances are evaluated in a unique time point in each operation with material transference, instead of considering the dynamic profiles of the input and output flows over a time interval, which characterize batch processes. This paper presents a novel modeling strategy to optimize the integrated batch process development problem using a comprehensive mixed-logic dynamic optimization (MLDO) model that evaluates all the interactions between synthesis and allocation

decisions simultaneously. Our goal is the systematic optimization of batch processes according to global objectives that embrace the process as a whole, by taking the maximum advantage of new and previously existing equipment in grassroots and retrofit scenarios. This way, we provide a tool for the agile development, improvement, or adaptation of chemical processes in rapidly changing and competitive markets. The integrated batch process development problem studied here (Section 2) comprises process synthesis decisions like the determination of process stages or tasks and their splitting into subtasks, the selection of technologies and chemicals, the definition of feed-forward control trajectories, duration of batch operations, and intermediate flow recycles, and the synchronization of operations with material transference. Simultaneously, the problem includes plant allocation decisions, namely the equipment selection, the task-unit assignment, the definition of the operating mode in single, series, or parallel operation, and the batch sizing. Finally, equipment capacities may represent either model constraints in existing processing items or further decisions on plant investments. The combination of all these degrees of freedom in an integrated optimization problem is solved in this paper for the first time without simplifications. Essentially, the modeling strategy (Section 3) relies on the representation of processing alternatives in a state-equipment network (SEN) superstructure (Smith and Pantelides, 1995) and its mathematical formulation as a MLDO problem (Section 4). In particular, the combination of hybrid discrete/continuous modeling and generalized disjunctive programming proposed by Oldenburg and Marquardt (2008) is extended herein to incorporate: (i) The combination of single-stage and multistage models to represent continuous and batch processing elements; (ii) The definition of multiple representation levels to link the operations of the different batch units that happen at the same time, whose ordering is not known beforehand; and (iii) The synchronization of the flow rates and compositions over time and the duration of operations with material transference, as a function of the selected process structure and tasks performance. To the authors’ knowledge, the combination of these modeling elements to handle simultaneously all the decisions previously discussed has not hitherto been applied in the context of batch process design. In consequence, the major novelty of this work is the ordering and combination of the coexisting single-stage and multistage dynamic models that represent the multiple processing units that compose alternative processing schemes. Previous contributions on process optimization have not addressed the problem solved in this work. For instance, some authors only study one single process stage (e.g. Oldenburg et al., 2003) or do not use dynamic variables in the material transference between equipment units (e.g. Ferrer-Nadal et al., 2008). In other works, physical constraints of existing equipment units are not considered, what implies that particular equipment constraints do not restrict the operation of the selected tasks (e.g. Charalambides et al., 1995), or each batch process stage is composed of a single batch operation instead of a chain of operations (e.g. Corsano et al., 2004). With this last simplification, the whole process is represented by a unique multistage model by associating one modeling stage with each task. In our work, all these assumptions are eliminated by addressing the synchronization of coexisting single-stage and multistage models. The technical difficulty to perform such a synchronization is the need of linking dynamic profiles of the input and output flows – some of them being control variables, like the input/output flow rates – over the duration of the operations with material transference. This is opposite to linking a single time point and requires the definition of the corresponding

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

time interval in the time axis of the semi-continuous elements. Additionally, the sequence of processing units is not known beforehand, as it depends on structural decisions like the task selection or the operating mode. The resulting MLDO model of our overall approach can be solved through several solution methods; a direct-simultaneous methodology is used in this work (Section 5). The benefits associated with the proposed optimization-based strategy are illustrated with two numerical examples, namely the Denbigh competitive reaction system (Section 6) and an acrylic fiber production process (Section 7). 2. Problem statement This work solves the batch process synthesis and the plant allocation sub-problems simultaneously in single-product campaigns, in order to define the process and the required equipment investment in new or existing manufacturing facilities, while enhancing the process efficiency and plant flexibility compared to the use of fixed recipes. The entire problem is defined as follows. Given: • Planning data: set of products, intermediates, and raw materials, expected demand for final products, and maximum time horizon; • Plant diagram: the SEN superstructure of available and potentially installed equipment units for each process stage or task, pipelines and connection nodes like mixers and splitters; • Task network: potential and mandatory tasks, alternative chemicals involved in each process stage – i.e. reactants, solvents, or catalysts –, allowed technologies, and possible reuse of intermediates in following batches; • Batch process operation: allowed task-unit assignments, batch operations within each unit procedure, operation to operation switching conditions, and set of limiting processing conditions of each equipment unit; • Process dynamics: DAE systems to represent the process behavior in each unit procedure, initial conditions, and set of process variables and dynamic or time-invariant controls; • Data related to performance evaluation: decision criteria and parameters to evaluate the objective function – e.g. selling price of final product, direct cost of raw materials, capital costs, amortization rates, operating expenses in processing units, environmental impact indicators; the goal is to determine: • Synthesis of processing schemes: selection of process stages and splitting into subtasks, technological specification, chemicals involved – i.e. reactants, solvents, or catalysts –, reference trajectories of the feed-forward control variables, duration of batch operations composing each process stage, recirculation of intermediate flows to be used in the next batch, and material transfer synchronization between tasks – i.e. synchronization of flow rates, compositions, and starting and final times; • Allocation of manufacturing facilities: selection of processing and storage units, task-unit assignment, operating mode – i.e. single, series, or parallel operation –, and batch sizes; • Plant design: equipment sizing; such that the adopted performance metrics are optimized. This problem can be solved in grassroots and retrofit scenarios. In the second case, the already existing equipment units have a predefined size value unless their capacity expansion is considered. For the sake of completeness, it is worth noting that there are other process development decisions which are out of the scope of this

291

study, as it is the case of decisions associated with multiproduct and multipurpose campaigns, e.g. batch order. 3. Modeling strategy The solution of optimization problems with structural decisions is usually addressed through the representation of all processing alternatives in a superstructure, which is later formulated into a mathematical model, to be finally optimized (Grossmann and Guillén-Gosálbez, 2010). In this procedure, the modeling approach is a key issue for the successful solution of the problem. This work presents a novel modeling strategy for representing all the degrees of freedom and interactions that characterize the integrated batch process development in a holistic optimization model. The approach is based on a bi-level SEN representation and its formulation as a MLDO problem. The background of mathematical tools for understanding such a strategy is presented next. 3.1. Background 3.1.1. Superstructure representation Batch design and scheduling problems require superstructure representations that combine plant, process, and material states. The most popular proposals are the state-task network (STN) (Kondili et al., 1993), the resource-task network (RTN) (Pantelides, 1994), and the state-equipment network (SEN) (Smith and Pantelides, 1995). Among them, the SEN entails the representation of equipment items and states of transferred material. This representation enables the direct formulation of equipment equations and clearly illustrates the dependence of the profiles of material flows on selected processing paths and conversion of preceding process stages. Precisely, other contributions in the literature have applied SEN superstructures to solve scheduling problems with process dynamics (Nie et al., 2012; Nie, 2014). 3.1.2. Hybrid discrete/continuous models These are dynamic models that represent batch process performance and batch events by combining discrete and continuous variables (Barton and Pantelides, 1994). Depending on the type of discontinuities represented, they are classified as: (i) singlestage models, without discrete events, (ii) multistage models, with explicit discontinuities, and (iii) general hybrid discrete/continuous models, with implicit and explicit discontinuities. Hybrid discrete/continuous models have been widely used in DO problems for determining the optimal control profiles of individual processing units over time (Binder et al., 2001; Srinivasan et al., 2003). In predefined sequences of process stages, dynamic control trajectories have also been optimized in the field of batch process synthesis (Allgor and Barton, 1999; Charalambides et al., 1995; Sharif et al., 1999), batch plant design (Bhatia and Biegler, 1996; Corsano et al., 2007), and short-term scheduling (Capón-García et al., 2013; Chu and You, 2014; Frankl et al., 2012; Nie et al., 2012), among others. 3.1.3. Logic-based modeling The incorporation of logical variables and equations into mathematical models allows the combination of quantitative and qualitative information. The use of logics in optimization problems was initially proposed by Balas (1985) in the so-called disjunctive programming, and rapidly became an alternative to mixed-integer modeling for representing structural decisions (Grossmann, 1985; Duran and Grossmann, 1986). These kinds of problems were later formalized as generalized disjunctive programming (GDP) (Raman and Grossmann, 1991; Floudas and Grossmann, 1994), where qualitative decisions are mathematically represented using Booleans and available process knowledge and heuristic rules are formulated through logical propositions. Mixed-logic formulations with

292

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

dynamic equations have contributed to the solution of several problems such as model predictive control (Bemporad and Morari, 1999), simultaneous design and control (Bansal et al., 2002), and scheduling of continuous processes with grade transitions (Prata et al., 2008). In the context of batch process development, Linninger and Chakraborty (1999) and Chakraborty and Linninger (2002) addressed the problem of batch process synthesis and plant allocation including logical rules in the superstructure formulation, but simplifying the problem through approximated process performance models. In contrast, Oldenburg and Marquardt (2008) and Oldenburg et al. (2003) presented a MLDO formulation that combined hybrid discrete/continuous models and GDP for solving the optimal control profiles, configuration, sequencing, and operation of batch equipment items with structural decisions. However, the formulation only considered individual units rather than task sequences, so the ordering and synchronization of processing units were not approached. 3.1.4. Coexisting single-stage and multistage models for synchronization purposes Single-stage and multistage models can be used to represent batch and semi-continuous procedures in the superstructure, respectively. Since it is necessary to guarantee the lot integrity, these kinds of models have to be synchronized in operations with material transference as a function of the selected sequence of processing units. This synchronization involves establishing a link between the dynamic profiles of the input and output flows and the duration of the operations with material transference. However, currently, there are neither modeling frameworks nor software tools which can handle coexisting multistage and single-stage models. Previous works on optimization of sequences of processing elements with dynamic profiles have not addressed the ordering and combination of several multistage models associated with batch units. This is either because the physical constraints of existing equipment units were not considered (Charalambides et al., 1995), or because batch process stages were simplified as a single batch operation and the whole process could be represented by a unique multistage model – i.e. by associating one modeling stage with each task (Corsano et al., 2004). 3.2. Bi-level SEN representation and MLDO modeling strategy The proposed modeling strategy is based on the combination of processing alternatives in a SEN superstructure, which is divided into two representation levels and formulated mathematically as a MLDO problem. A basic SEN superstructure illustrates the interconnection between the existing and potential equipment units and the states of material flows according to the plant diagram, e.g. Fig. 1a. This superstructure representation enables the formulation of physical plant restrictions and procedures of processing elements as well as the definition of the profiles of material flows as a function of the selected units, their order, and their process conversion. The mathematical formulation of the superstructure and all the associated decisions into a MLDO model comprises disjunctions and logical propositions, which define the structural decisions, as well as hybrid discrete/continuous models – referred to as multistage models in this paper –, which determine the process performance in each unit. This formulation is inspired by the approach by Oldenburg and Marquardt (2008) for solving the optimization of individual units. Here, the MLDO formulation is extended to combine, organize, and synchronize the multiple equipment units in the SEN superstructure according to the processing alternative that is selected. With this purpose, each unit is represented by either a single-stage or a multistage model associated with a semi-continuous or batch procedure, respectively.

Fig. 1. Basic SEN superstructure (a) and bi-level superstructure (b) of the motivating example 1.

Therein, each mathematical stage corresponds to a batch operation. Moreover, in order to perform the coordination and timing of consecutive units, these multistage and single-stage models are nested in their entirety inside specific terms of the disjunctive equations, which determine their usage and their order. For visualizing this concept, the plant elements and their models are distributed into two representation levels by dividing the SEN superstructure into Level 1, containing batch units, and Level 0, with connection nodes, storage tanks, and semi-continuous processing units, as illustrated in Fig. 1b for the previous example. Independently of their representation level, all the models are optimized at the same time. Overall, this modular bi-level representation allows the ordering of processing elements as a function of decisions like the task-unit assignment because, this way, the time coordinates of batch units within the integrated optimization problem can be moved in relation to the time horizon of other batch units and the complete process. The novel coexistence and synchronization of single-stage and multistage models are illustrated in Fig. 2a. The time links across the different models are established through disjunctive equations and logical propositions in the formulation. Additionally, for enabling the relations across the mathematical stages of the different models using standard optimization tools, the equations that compose the optimization problem are reformulated as following detailed: 1. Transformation of single-stage models at Level 0 into multistage models. The connection between the models of the different processing units requires the conversion of single-stage models at Level 0 into multistage models that include all the possible batch operations at Level 1, as represented in Fig. 2b. With this purpose, the sets of equations of the original singlestage models are replicated for every possible new stage, and the continuity of process variables is guaranteed by incorporating stage-to-stage relations. 2. Normalization of mathematical stages. Since the several processing elements at Levels 0 and 1 may be composed of a

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

293

Fig. 2. Multistage and single-stage models of the motivating example 1 to represent batch and semi-continuous plant elements respectively.

diverse number of mathematical stages with different duration (see Fig. 2b) which should be solved simultaneously, it is necessary to normalize their time axis. This way, all the stages of every model can be solved at the same time, as shown in Fig. 2c. The normalization of time intervals is carried out by multiplying the differential equations by the stage duration and redefining their derivation interval from 0 to 1. Except for the actual first stage, the initial conditions of all other stages should be treated as degrees of freedom whose values are subject to the fulfillment of stage-to-stage continuity equations. 3. Explicit definition of stage durations. The normalization of the time intervals defined in the previous step has a collateral advantage of expressing the stage durations as explicit variables. This way, it is possible to match the timespans of the mathematical stages of synchronized operations. Beyond these reformulations of the model equations, the identity of each mathematical stage at Levels 0 and 1 is preserved because it is a key issue for the synchronization of the selected processing units.

4. Mathematical formulation The variables and equations that compose the optimization model for integrated batch process development are detailed in this section according to the modeling strategy mentioned above (Section 3). The key sets, parameters and variables required to understand the formulation are first presented.

Fig. 3. Synchronization of batch stages and flow rates at Levels 0 and 1 for configurations ˛, ˇ, , and of the motivating example 1, with |Kj | = 3, ∀ j ∈ {U1 , U2 }.

4.1. Notation 4.1.1. Mathematical stages Mathematical stages k ∈ Kj at Level 1 should be linked with stages l ∈ L at Level 0 to synchronize the duration of batch operations in units j ∈ U. However, the number of total mathematical divisions that are necessary at Level 0 for this purpose depends on the selected tasks i ∈ PS and operating modes ∈ i . For instance, the utilization of a single batch unit U1 in the superstructure of Fig. 1 could involve the three operations load, hold, and unload of unit U1 , whereas two additional stages would be required in series operation to complete the hold and unload operations of a second batch unit U2 , as illustrated in Fig. 3a (configuration ˛) and d (configuration ). The maximum number of stages Lmax concerning all processing alternatives is determined by the structural system composed of all the batch units in series. The total time Tf at Level 0 is then divided into Lmax intervals with different durations. Out of

294

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

Fig. 5. Petri net representing process and bypass stages k ∈ Kj in the model of batch unit j ∈ U at Level 1. Circles symbolize stages and bars represent stage-to-stage transitions. Fig. 4. Input/output variables and stages of batch unit j ∈ U models, considering min ∈ Mjin input and mout ∈ Mjout output flows, and c ∈ C chemicals.

the resulting set of stages L = {1, . . ., Lmax }, only the subset La composed by the so-called active stages is in force in each structural alternative. The connection between Levels 0 and 1 is especially relevant to stages that represent operations with material transference. Particularly, these are termed input and output stages and are denoted by Ij ⊆ Kj and Oj ⊆ Kj respectively. 4.1.2. Input and output flows Additionally, input and output flows of batch units are represented at both Level 0 (denoted by n ∈ Njin ∪ Njout ) and Level 1 (denoted by m ∈ Mjin ∪ Mjout ) and are linked between the two levels in the formulation. Essentially, these flows are characterized by j j flow rates Fm,k (t) and compositions xc,m,k (t), which are referred to as input and output variables. As illustrated in Fig. 4, the inflows and outflows of batch units should be restricted to zero in those stages that are not input Kj \Ij or output Kj \Oj stages respectively, to guarantee the absence of material transmission. 4.1.3. Control variables There are four types of decision or control variables: 1. Dynamic control variables udyn (t): flow rates in input and outj put connections Fm,k (t), m ∈ Mjin ∪ Mjout and internal control j

variables int k (t) over operations k ∈ Kj in batch units j ∈ U; 2. Static or time-invariant decision variables ustat : duration of batch operations tl represented by mathematical stages l at Level 0; 3. Integer decision variables uint : number of batches NBp of product p ∈ P and discrete sizing Sizej of equipment units j ∈ U; 4. Boolean decision variables uBool : selection of tasks Zi , i ∈ PS, processing and storage units Yj , j ∈ U ∪ T, operating modes X i , ∈  in process stages i ∈ PS, processing orders Wj,q , q ∈ j

Q, technological specifications V ,  ∈ , and chemical comj

pounds Sc , c ∈ C of processing units j ∈ U, and recirculation of intermediate flows Rn , n ∈ Nr . The consideration of all these decisions may lead to overspecified problems in some cases. Therefore, it is possible to improve the efficiency of the search procedure and reduce convergence issues, by removing the extra degrees of freedom, which are mostly connected with logical decisions and dynamic control variables in operations with material transference. The reduction of the number of control variables is analyzed at the end of this section based on the optimization model. 4.1.4. Other variables and parameters The problem comprises other types of variables, namely differ˙ ential variables z(t) and their derivatives z(t), algebraic variables y(t), time-invariant or static variables , and parameters p. When

these variables refer specifically to equipment units or mathematical stages, they may be written with subindices j or k respectively. 4.2. Batch procedures at Level 1 4.2.1. Batch units Let us consider discontinuous units j ∈ U located at Level 1 of the superstructure. Their operation and process performance are represented by a |Kj |-stage model. For each unit, a two-term disjunction controlled by Boolean Yj is introduced. If Yj is true, this processing item is selected to allocate one task, and the |Kj |-stage model in the corresponding disjunctive term is activated. On the contrary, if Yj is false, the unit is not selected and a bypass strategy (Oldenburg and Marquardt, 2008) is applied in order to define an equivalent number of |Kj | mathematical stages, termed bypass stages, where no process takes place. Bypass stages are represented in the Petri net of Fig. 5. They are characterized by stage durations j tk equal to zero and a set of equations that override the dynamic model of the batch process. Overall, the disjunctive equation for batch unit j ∈ U is defined by the following form, taking into account that the time interval has been normalized (see Fig. 2c):

⎡

⎤

Yj

⎢ fj,kd (z˙ j,k (t), zj,k (t), yj,k (t), udyn (t), ustat , uint , pj ), t ∈ [0, 1], ∀ k ∈ Kj , ⎥ j,k ⎢ ⎥ ⎢ ⎥ ljd (z˙ j,1 (0), zj,1 (0)), ⎢ ⎥ ⎢ d ⎥ stat int ⎢ gj,k (zj,k (t), yj,k (t), udyn ⎥ (t), u , u , p ) ≤ 0, t ∈ [0, 1], ∀ k ∈ K , j j j,k ⎢ ⎥ ⎢ ⎥ dyn d,e stat int gj,k (zj,k (1), yj,k (1), uj,k (1), u , u , pj ) ≤ 0, ∀ k ∈ Kj , ⎢ ⎥ ⎢ ⎥ d ⎢ ⎥ zj,k+1 (0) − mj,k (zj,k (1)) = 0, ∀ k ∈ {1, . . ., |Kj | − 1}, ⎢ ⎥ ⎢ ⎥ dyn j = hdj (zj,|Kj | (1), yj,|Kj | (1), uj,|K | (1), ustat , uint , pj ), ⎢ ⎥ j ⎢ ⎥  ⎣ ⎦ k ∈ Kj

⎡

(1)

j

tk , t j,end = t j,s + T j,f

T j,f =

⎤

¬Yj

⎢ ⎢ Bd (z˙ j,k (t), zj,k (t), yj,k (t), udyn (t), ustat , uint , j , pj ) = 0, j,k ⎢ j ⎢ j ⎢ tk = 0, ∀ k ∈ Kj , ⎣ j,s

T j,f = 0, t j,end = 0, tj = 0

⎥

t ∈ [0, 1], ⎥ ⎥

⎥, ⎥ ⎦

∀ j ∈ U,

d , g d , g d,e , and md are the DAE system, path constraints where fj,k j,k j,k j,k (PC), end-point constraints (EPC), and stage-to-stage continuity in stage k ∈ Kj of unit j ∈ U. ljd are the relations to determine the ini-

tial conditions, hdj is the set of equations to calculate time-invariant variables  j evaluated at the final time (t = 1) of the last stage |Kj |, and Bjd are the equations that define the system in bypass stages. The dynamic model of unit j ∈ U comprises time-dependent derivatives z˙ j,k (t) and differential zj,k (t), algebraic yj,k (t), and condyn

j

trol uj,k (t) variables – including input/output flow rates Fm,k (t) and j xc,m,k (t)

j

and internal control variables int k (t) – in compositions stages k ∈ Kj . It is also characterized by process parameters pj and

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

295

Fig. 6. Basic SEN superstructure of the motivating example 2.

time-invariant variables  j , necessary for the objective function and other key performance indicators (KPIs), like product selectivity or process conversion in j ∈ U. To complete batch unit models at Level 1, input/output variables should be dismissed in those stages which are not input/output stages respectively, as illustrated in Fig. 4, by restricting their value to zero. With this purpose, the following equation is defined:

⎡

⎤

Yj

(2)

j

Fmout ,k (t) = 0, ∀ mout ∈ Mjout , ∀ k ∈ Kj \Oj j

Additionally, the volume k (t) at unit j ∈ U cannot surpass the equipment capacity Sizej in any stage k ∈ Kj , either Sizej is a free decision variable to design a newly installed unit or an expansion or it is a fixed variable associated with an existing item. This constraint is formulated by: t ∈ [0, 1], ∀ k ∈ Kj , ∀ j ∈ U.

(3)

4.2.2. Technological specifications In some cases, different processing technologies  can be used in the same unitary operation. The selection between alternative technologies requires the representation of different equipment arrangements that determine the particular physicochemical equations, parameters, and batch operations sequences for each case. Hence, every alternative has to be associated with a specific processing item j ⊆ , |j | = 1, ∀ j ∈ U, in order to characterize its DAE system through Eq. (1). For instance, Fig. 6 shows a SEN with two technological alternatives A and B in equipment units U1,A and U1,B respectively. The technology  ∈ j associated with j

equipment unit j ∈ U is represented by Boolean V , which indicates

j whether such specification is selected (V

j

∀  ∈ j , j ∈ U.





c ∈ Cs



j

Sc pj = pj,c

j

,

∀ j ∈ U.

(5)

4.3.1. Processing order of batch units The time axis for the models of batch units j ∈ U that are selected (Eqs. (1) and (2) with Yj = true) should be moved over the time axis of the entire process at Level 0 according to their processing order, as illustrated in Fig. 3. Then, each batch unit model is synchronized with the models of other plant elements by relating the set of stages Kj at Level 1 to specific stages l ∈ L at Level 0. Assignment Booleans Wj,q are defined to lead the synchronization by indicating the order q ∈ Q of unit j ∈ U with reference to other batch units. Once a unit is selected, this should be assigned one and only one processing order q ∈ Q: Yj ⇔  Wj,q , q∈Q

(4)

4.2.3. Chemical compound alternatives The selection of chemicals like reactants, solvents, or catalysts c ∈ Cjs that are involved in the process of unit j ∈ U is represented by

(6)

• The starting time tj,s of unit j ∈ U (Level 1) is calculated from the 0 general starting time ts and duration tl of stages that precede lj,q (Level 0) by:

⎡



q∈Q

⎣

⎤

Wj,q t j,s = t s +

0 −1 ⎦ , ∀ j ∈ U; lj,q

l=1

(7)

tl

• The duration t j of stage k ∈ Kj of unit j ∈ U (Level 1) has to be k equal to the duration tl of the corresponding stage l ∈ L (Level 0) by:

⎡

j

Boolean Sc , which indicates whether a chemical compound is used j j (Sc = true) or not (Sc = false). The SEN of Fig. 6 illustrates an example with two possible compounds c1s and c2s in unit U2 . Unlike the case of alternative technologies  ∈ j , the choice of a chemical compound c ∈ Cjs does not necessarily affect the complete set of equations in the model of unit j ∈ U (Eq. (1)) but only the set of parameters pj . In this case, such parameters are defined by the values pj,c associated

∀ j ∈ U,

0 ∈ L0 . This way, the relaand its operation should start in stage lj,q j tion between stage k ∈ Kj at Level 1 and its corresponding stage l ∈ L at Level 0 can be expressed as a function of the processing order 0 − 1 ∈ L. On this q ∈ Q of unit j ∈ U through the function l = k + lj,q basis, the variables of models at Levels 0 and 1 are connected as follows:

j

= true) or not (V = false) and is related to the equipment Boolean Yj as follows:

Yj ⇔ V ,

j

Sc . This constraint is defined as follows:

4.3. Synchronization

⎢ j ⎥ ⎢ F in (t) = 0, ∀ min ∈ Mjin , ∀ k ∈ Kj \Ij , ⎥ , ∀ j ∈ U. ⎣ m ,k ⎦

jk (t) ≤ Sizej ,

with alternative c ∈ Cjs according to the chemical selection Boolean



q∈Q

⎣

⎤

Wj,q j , l−l0 +1

tl = t

j,q

⎦ , ∀ j ∈ U; (8)

0 , . . ., |K | + l0 − 1} ∀ l ∈ {lj,q j j,q

• The profiles of input and output variables of unit j ∈ U (Level j j 1) – i.e. Fm,k (t) and xc,m,k (t), ∀c ∈ C, ∀m ∈ Mjin ∪ Mjout – have to be equal to the profiles of the analogous variables at Level

296

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

0 – i.e. Fn,l (t) and xc,n,l (t), ∀c ∈ C, ∀n ∈ Njin ∪ Njout . Additionally, the latter ones have to be fixed to zero in the remaining nonsynchronized stages, as well as in the case that unit j ∈ U is not selected. The resulting equation for every pair of equivalent flows at Levels 0 and 1 (n, m) ∈ (Njin , Mjin ) ∪ (Njout , Mjout ) reads as:

⎡

⎤

Wj,q

j ⎢ Fn,l (t) = F j (t), xc,n,l (t) = x (t), t ∈ [0, 1], ⎢ m,l−l0 +1 c,m,l−l0 +1 j,q j,q ⎢ 0 0  ⎢ ∀ l ∈ {lj,q , . . ., |Kj | + lj,q − 1}, q∈Q ⎢ ⎣ Fn,l (t) = 0, xc,n,l (t) = 0, t ∈ [0, 1], 0 0 ∀ l ∈ L\{lj,q , . . ., |Kj | + lj,q − 1}

⎡ ¬Y

⎤

⎢

⎥, ⎦

j

 ⎣ Fn,l (t) = 0,

xc,n,l (t) = 0, t ∈ [0, 1], ∀ l ∈ L

Fig. 7. Petri net representing active stages l ∈ La = {1, . . ., |La |} and bypass stages l ∈ L\La = {|La | + 1, . . ., Lmax } for semi-continuous elements at Level 0. Gray elements represent feasible paths for other locations of stage |La |.

⎥ ⎥ ⎥ ⎥ ⎥ ⎦

⎡ ⎣



∈ i

∀ j ∈ U.

(9)

⎤

Xi Fn,l (t) = 0, t ∈ [0, 1], ∀ n ∈ Ni,0 , ∀ l ∈ L

⎦,

∀ i ∈ PS.

(16)

4.5. Plant elements at Level 0 4.4. Process stages 4.4.1. Process stages Logical variable Zi represents the selection of task i ∈ PS and, therefore, controls whether the input flow n ∈ Niin is processed or

whether it is bypassed through pipeline b ∈ Nib to the following task. In the problem superstructure of Fig. 6, this decision is related to splitter Sp3 . The two exclusive alternatives are defined through the following disjunction:





Zi Fb,l (t) = 0



¬Zi



Fb,l (t) = Fn,l (t)

, t ∈ [0, 1],

(10)

∀ l ∈ L, ∀ n ∈ Niin , ∀ b ∈ Nib , ∀ i ∈ PS.

following proposition:  Xi ,

∈ i

(11)

The principal purpose of configuration Booleans is to govern the selection of equipment items j ∈ U, represented by variable Yj , and their processing order q ∈ Q, represented by variable Wj,q . A logical proposition is used to define each configuration. For example, configurations i = {˛, ˇ, , } which correspond to single operation in U1 , single operation in U2 , parallel operation in U1 and U2 , and series operation in U1 followed by U2 in the example of Fig. 1, are defined by: X˛i

⇔ WU1 ,1 ∧ ¬YU2 ,

(12)

Xˇi ⇔ WU2 ,1 ∧ ¬YU1 ,

(13)

Xi ⇔ WU1 ,1 ∧ WU2 ,1 ,

(14)

X i

(15)

⇔ WU1 ,1 ∧ WU2 ,2 .

4.5.1. Semi-continuous elements The time variables at Level 0 are related as follows:



tl , (17)

l∈L t end = t s + T f ,

where Tf is the total batch processing time which starts at initial time ts , ends at final time tend , and is divided in Lmax intervals with duration tl . Following the bypass strategy, the stage duration is zero in extra stages l ∈ L\La , whereas it is confined between lower tL and upper tU bounds in active ones l ∈ La : t L ≤ tl ≤ t U ,

∀ l ∈ {1, . . ., |La |},

bypass stages:tl = 0,

∀ l ∈ {|La | + 1, . . ., Lmax }.

(18)

Additionally, material balances in connecting nodes at Level 0 relate flow rates Fn,l (t) and compositions xc,n,l (t), c ∈ C of connections n ∈ N in every stage l ∈ L. In particular, mixers j ∈ Mx are characterized by several input connections nin ∈ Njin , |Njin | > 1, and a single output nout ∈ Njout , |Njout | = 1, and are described by:



Fnin ,l (t) = Fnout ,l (t),

t ∈ [0, 1],

nin ∈ N in j

Additionally, configuration Booleans X i enforce specific flow distributions with each operating mode. To do so, the flow rate of a specific set of connections n ∈ Ni,0 in each configuration ∈ i is restricted to zero, according to:

are not required at Level 0 are overridden by forcing stage durations tl to zero and withdrawing dynamic equations and constraints of Level 0 elements.

Tf =

4.4.2. Operating mode Out of the set of allowed operating modes or configurations ∈ i in a process stage i ∈ PS, with the corresponding Boolean variables X i , only one can be selected. This is defined through the

Zi ⇔

The bypass strategy, previously applied to batch units j ∈ U (see Fig. 5), is also used at Level 0 to define the set of extra stages L\La , out of the total number of stages Lmax , that should be disregarded in each structural alternative. The reason is that the set of active stages La representing the sequence of batch operations in the complete process depend on the selected tasks Zi and their configuration X i . As shown in Fig. 7, the set of stages l ∈ L\La that

∀ nout ∈ Njout , ∀ l ∈ L, ∀ j ∈ Mx,



Fnin ,l (t)xc,nin ,l (t) = Fnout ,l (t)xc,nout ,l (t),

(19)

t ∈ [0, 1],

nin ∈ N in j

∀ c ∈ C, ∀ nout ∈ Njout , ∀ l ∈ L, ∀ j ∈ Mx,

(20)

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

whereas splitters j ∈ Sp are characterized by a single input connection nin ∈ Njin , |Njin | = 1, and several outputs nout ∈ Njout , |Njout | > 1, and are described by:



Fnin ,l (t) =

Fnout ,l (t),

t ∈ [0, 1],

nout ∈ N out j

∀n

in

∈

∀ l ∈ L, ∀ j ∈ Sp,

Njin ,

(21)

297

4.5.2. Recirculation of intermediate flows The recirculation of an intermediate material to be used in a following batch is represented by connection n ∈ Nr . In particular, this decision is controlled by Boolean Rn , which indicates if the intermediate flow is recirculated (Rn = true) or not (Rn = false). As a result, the flow rate can either acquire a value between the lower FnL and upper FnU bounds or be canceled:





Rn FnL ≤ Fn,l (t) ≤ FnU , t ∈ [0, 1], ∀ l ∈ L

xc,nin ,l (t) = xc,nout ,l (t),



t ∈ [0, 1],

∀ c ∈ C, ∀ nin ∈ Njin , ∀ nout ∈ Njout , ∀ l ∈ L, ∀ j ∈ Sp.

(22)





¬Rn

,

Fn,l (t) = 0, t ∈ [0, 1], ∀ l ∈ L

∀ n ∈ Nr .

(25)

It is not necessary to remove these equations in bypass stages because all flow rates would be forced to zero through previous Eqs. (9), (19) and (21). Storage units j ∈ T are also defined at Level 0 by |L|-stage modd , ld , g d , g d,e , hd , els. Such models comprise the general functions fj,l j j,l j j,l

The recirculation of intermediate flows additionally requires the installation of buffer tanks j ∈ Tnr in the superstructure – e.g. tank r in the example of Fig. 6 – in order to store the intermediate T11 material that is later supplied in a subsequent batch. The logical proposition reads as:

and mdj,l for storage tanks j ∈ T and mathematical stages l ∈ La . Since storage tanks operate continuously, all the equations are d = fd d = gd repeated in active stages as defined by fj,l , gj,l , and j,l+1 j,l+1

Rn ⇔ Yj ,

d,e gj,l

d,e gj,l+1 , ∀l

{1, . . ., |La | − 1}.

= ∈ Following, the bypass strategy is applied to two situations: First, if storage unit j ∈ T is not selected, its complete model is deactivated like in the case of batch units j ∈ U. Second, the bypass method is required to differentiate active La and bypass L\La stages according to Fig. 7, provided that this unit is selected. The resulting formulation of storage tanks j ∈ T reads as:

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

⎤

Yj dyn d fj,l (z˙ j,l (t), zj,l (t), yj,l (t), uj,l (t), ustat , uint

, pj ), t ∈ [0, 1], ∀ l ∈ {1, . . ., |L |}, a

ljd (z˙ j,1 (0), zj,1 (0)), dyn

d gj,l (zj,l (t), yj,l (t), uj,l (t), ustat , uint , pj ) ≤ 0, t ∈ [0, 1], ∀ l ∈ {1, . . ., |La |}, dyn

d,e gj,l (zj,l (1), yj,l (1), uj,l (1), ustat , uint , pj ) ≤ 0, ∀ l ∈ {1, . . ., |La |},

zj,l+1 (0) − mdj,l (zj,l (1)) = 0, ∀ l ∈ {1, . . ., |La | − 1}, dyn

j = hdj (zj,|La | (1), yj,|La | (1), uj,|La | (1), ustat , uint , pj ), bypass stages: dyn Bjd (z˙ j,l (t), zj,l (t), yj,l (t), uj,l (t), ustat , uint , j , pj ) = 0, t ∈ [0, 1],

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦



(23)

¬Yj

In order to define mathematically the sequence of storage and usage of intermediate material, initial conditions zj,l (t = 0) in the first stage l = 1 of unit Tnr are set to be equal to the final conditions zj,l (t = 1) in the last stage l = Lmax , what reads as:



Rn zj,1 (0) = zj,|Lmax | (1), ∀ j ∈ Tnr



,

∀ n ∈ Nr .

(27)

4.6.1. Number of batches The batching problem consists in the division of the total product demand into a number of batches of a specific production size. Common approaches for solving scheduling problems first address the batching activity and following solve the allocation, timing, and task sequencing sub-problems. In this work, the simultaneous solution of batch process development incorporates the optimization of the number of batches by means of the following equation: NBp Batchp + Shortfallp ≥ Demandp ,

∀ p ∈ P,

(28)

where Demandp is the total demand of product p ∈ P, NBp represents the number of batches, Batchp stands for the batch production size, and Shortfallp is the unaccomplished demand.

⎤

⎣ Bjd (z˙ j,l (t), zj,l (t), yj,l (t), udyn (t), ustat , uint , j , pj ) = 0, ⎦ , ∀ j,l

(26)

4.6. Batching

∀ l ∈ {|La | + 1, . . ., Lmax }

⎡

∀ j ∈ Tnr , ∀ n ∈ N r .

j ∈ T,

t ∈ [0, 1], ∀ l ∈ L

4.7. Objective function and degrees of freedom

dyn

where z˙ j,l (t), zj,l (t), yj,l (t), and uj,l (t) are the time dependent derivatives and differential, algebraic, and control variables in mathematical stage l ∈ L, ustat are time-invariant control variables, namely stage durations, uint are integer decisions such as the storage tank size, pj are process parameters, and  j are time-invariant variables which may contribute to the evaluation of the objective function or KPIs in unit j ∈ T. Finally, a constraint to limit the maximum volume should be also formulated, which reads as:

The optimization problem presented above is solved to minimize the objective function :



˙ minimize  z(t), z(t), y(t), udyn (t), ustat , uint , uBool , , p , (29) udyn (t), ustat uint , uBool

(24)

˙ defined by derivatives z(t), differential variables z(t), algebraic variables y(t), time-invariant variables  and parameters p, j j according to the control variables udyn (t) = {Fm,k (t), int k (t)},

is the inventory volume of unit j ∈ T in stage l ∈ L and where Sizej is the maximum capacity of the tank. The decision on the selection of intermediate tanks is crucial for allowing the recirculation of intermediate flows as well as for defining more than one cycle time or batch size for different parts of the process.

Sc , Rn }. The degrees of freedom of the optimization model are the number of decision variables remaining after subtracting the number of equations and predefined values. This parameter determines the number of practical control variables in the reality. In the

jl (t) ≤ Sizej , jl (t)

t ∈ [0, 1], ∀ l ∈ L, ∀ j ∈ T,

j

ustat = {tl }, uint = {NBp , Sizej } and uBool = {Zi , Yj , X i , Wj,q , V , j

298

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311 j

j

case of Booleans uBool = {Zi , Yj , X i , Wj,q , V , Sc , Rn }, the number of degrees of freedom is limited by logical propositions (Eqs. (4), (6), (9), (11), (12)–(15), and (26)). As for the dynamic control profiles of j flow rates in batch units Fm,k (t) ⊆ udyn (t), these profiles are interrelated through global balances in mixers and splitters (Eqs. (19) and (21)). Therefore, the evaluation of degrees of freedom of dynamic profiles between consecutive units is especially relevant because the input and output flow rates of different units j ∈ U that are synchronized cannot behave as decision variables simultaneously. In fact, the system over-specification deteriorates the performance of the solution search procedure. The degrees of freedom related to flow rates are denoted by DOFi, and differ for each configuration ∈ i and process stage i ∈ PS because they depend on the total number of flow rates at Level 0 |Ni |, the number of equations in splitters and mixers |Spi ∪ Mxi |, and the number of predefined flow rates, namely the f flow restrictions |Ni,0 | in configuration ∈ i and input flows |Ni | defined in preceding tasks. Overall, the parameter DOFi, puted as follows: DOF i,

is com-

= No. variables − No. equations − No. fixed variables f

= |Ni | − |Spi ∪ Mxi | − |Ni,0 | − |Ni |.

(30)

To adjust the number of control variables according to the degrees of freedom of each structural alternative, outflows of batch units are always considered as free-decision variables: j

Fmout ,k (t) ∈ udyn (t),

5.1. MLDO reformulation

t ∈ [0, 1],

∀ mout ∈ Mjout , ∀ k ∈ Oj , ∀ j ∈ Ui , ∀ i ∈ PS,

(31)

while only the inflows of the subset of batch units Di, ⊆ Ui are defined as control variables in configuration ∈ i of process stage i ∈ PS. The set Di, is defined such as |Di, | = DOFi, − |Ui |, where |Ui | corresponds to the degrees of freedom assigned to the output flows in Eq. (31). The resulting equation reads as:

⎡

Xi

equations (Luus, 1990), (ii) indirect methods through necessary conditions of optimality (NCO) derived from Pontryagin’s formulation (Bryson and Ho, 1975), and (iii) direct methods which convert the time-continuous optimization problem into a finitedimensional nonlinear programming problem by discretization. Notably, the advantage of the direct approaches is that they do not require an explicit derivation of the necessary conditions of optimality and the adjoined equations and, after the discretization of time-dependent variables of the original problem, provide a mathematical programming (MP) problem. There exist different direct methods, depending on what parts of the problem are discretized. In the specific case of the direct-simultaneous strategy (Neuman and Sen, 1973; Tsang et al., 1975), there is a full discretization of the dynamic model, as both process and control variable profiles are approximated by a series of finite points along the time horizon. The resulting problem is a mixed-integer non-linear programming (MINLP) model that can be solved without requiring sensitivities, second order derivatives computation, and obligatory variable continuity, as opposed to direct-sequential approaches (Sargent and Sullivan, 1978). However, their disadvantage is that the model size can increase exponentially in the discretization step. This work employs the classical direct-simultaneous strategy detailed in the next subsections. The reader interested in alternative stochastic and hybrid solution strategies is addressed to Moreno-Benito et al. (2015).

⎤

⎢ j ⎥  ⎢ F in (t) ∈ udyn (t), t ∈ [0, 1], ⎥ , ∀ i ∈ PS. ⎦ ∈ i ⎣ m ,k

(32)

∀min ∈ Mjin , ∀k ∈ Ij , ∀j ∈ Di,

5. MLDO solution Different kinds of solution strategies can be used to tackle the integrated batch process development formulated as a MLDO problem (Oldenburg and Marquardt, 2008; Stein et al., 2004). The majority of deterministic methods guarantee to find a local optimum and most advanced algorithms are also able to provide global solutions under fairly general assumptions, e.g. (Beaumont, 1990; Sahinidis, 1996). In particular, deterministic approaches are based on the reformulation of the mixed-logic problem into a mixedinteger one, the reformulation of Booleans as continuous variables, or the problem solution with logic-based solution techniques. The methods of the first group are also known as classical solution approaches and are the most extended ones. This is because the problem reformulation allows exploiting the advantages of mixedlogic modeling – e.g. the natural inference of previous knowledge and rules about the process and their direct conversion into mathematical equations – as well as using well established MIDO solution strategies. Among the classical methods for solving MIDO problems, we find three different groups, namely: (i) Dynamic Programming (DP) based on Hamilton–Jacobi–Bellman formulation to transform the original problem into a system of partial differential

To address the relaxation of the original MLDO model into MIDO one in classical solution methods, several transformations are  required in  the formulation. First, Boolean variables uBool ∈ true, false are replaced by binaries ubin ∈





0, 1 . Thus, the vector of logical decision variables uBool = j

j

{Zi , Yj , X i , Wj,q , V , Sc , Rn } is transformed into a vector of binary j

j

variables ubin = {zi , yj , xi , wj,q , v , sc , rn }, whose 0 and 1 values correspond to prior false and true values respectively. Second, disjunctive equations that enclose other constraints are transformed into mixed-integer equations. In the MLDO formulation presented in the previous section, these equations are Eqs. (1), (2), (5), (7)–(10), (16), (23), (25), (27) and (32). There are several alternatives for performing this transformation, like the big-M (Raman and Grossmann, 1991) and the convex-hull reformulation (CHR) (Türkay and Grossmann, 1998). In this work, we use binary multiplication. This technique is based on the decomposition of all the variables that are contained in the disjunctive terms into variable contributions, except for control variables udyn (t), ustat , and uint . In particular, the variables that are decomposed are of the fol˙ lowing types: derivatives z(t), differential z(t), algebraic y(t), and time-invariant  variables. After their decomposition, each variable contribution is associated with the constraints inside a disjunctive term and its corresponding Boolean variable. For instance, let us consider the following disjunctive equation:



d∈D

Ad Fd (x, u) ≤ 0

,

Ad ∈ {true, false}, x ∈ R, u ∈ R, (33)

where D is the set of disjunctive terms in the equation, Ad are the Boolean variables that control each of disjunctive terms d ∈ D, Fd are the functions nested in each disjunctive term d ∈ D, x are the dependent variables, and u are the decision variables. Through substituting the Booleans Ad ∈ {true, false} by the binaries ad ∈ {0, 1}

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

and applying binary multiplication, the variables x are decomposed into variable contributions xd and Eq. (33) is reformulated as: Fd (xd , u) ≤ 0, x=



xd ad ,

d∈D

(34)

ad = 1,

d∈D

ad ∈ {0, 1}, x ∈ R, u ∈ R. Finally, logical propositions, namely Eqs. (4), (6), (11)–(15) and (26) of the formulation, are expressed as linear constraints. Such a transformation can be done systematically by formulating the conjunctive normal form (CNF) of the original logical equations. This way, it is possible to obtain an expression like C1 ∧ C2 ∧ . . . ∧ CN , where Cn , n ∈ {1, . . . N} are the clauses that must be true in the problem, related by “and” operators (∧). This procedure involves the application of a series of pure logical operations, which were first formalized by Clocksin and Mellish (1981), as reported by Raman and Grossmann (1991). The operations are detailed in Appendix A. Next, the resulting true-clauses Cn can be transformed into algebraic equality or inequality equations using a number of equivalences between logical and algebraic expressions (Raman and Grossmann, 1991), which are summarized in Table 5 of Appendix A. 5.2. MIDO solution The solution of the resulting MIDO problem is addressed through the full discretization of the model using the orthogonal collocation method. This technique was originally introduced by Cuthrell and Biegler (1989) to solve optimal control problems with discontinuous control profiles, and renders attractive stability, symmetry, and accuracy properties (Cuthrell and Biegler, 1989, Section 3.1). The strategy consists of dividing the time axis into several intervals – referred to as finite elements e ∈ {1, . . ., Ne } – and specific time points – termed collocation points m ∈ {1, . . ., M} – and approximating the state and control variable profiles in these locations, as represented in Fig. 8. The position of the collocation points is determined by computing the roots of an orthogonal polynomial, like the shifted Legendre polynomial of order M. Then, the state and control variable profiles are approximated using monomial basis representations, such as Lagrange polynomials. Further details about the calculation of the collocation points and the construction of the orthogonal polynomials used in the numerical examples are provided in Appendix B. As illustrated in Fig. 8, boundary conditions are defined for differential variables z(t) in each finite element e to enforce their continuity from the previous element e − 1 in time te with the symmetry properties of the method. For this, it is necessary to establish a larger order Mth

299

in the polynomials that define the differential variables, compared to the polynomials of algebraic y(t) and control variables udyn (t), with order (M − 1)th (see Eqs. (B.2) and (B.3) in Appendix B). In contrast, continuity across neighboring elements is not necessary for algebraic and control variables. In fact, control variables are approximated to a piecewise constant (PWC) function. Overall, the quality of the approximation of the dynamic model is determined by the number of collocation points M and finite elements Ne , and so is the complexity of the resulting MINLP problem. The number of collocation points is a critical parameter because it defines the order of the resulting polynomials. It is usually selected considering the trade-off between being large enough for securing a good approximation of the discretized variables and small enough for limiting the increase in the computational demand associated with high-order polynomials. As for the number of finite elements, this factor can be adjusted according to the length of the time interval to be discretized. In this work, the full discretization approach is applied to each of the mathematical stages k ∈ Kj , j ∈ U and k ∈ L that compose the multistage models of the formulation using a predefined number of finite elements and collocation points. Since the duration of the batch operations is a decision variable, the length of the finite elements can vary after scaling the time intervals of the optimal solution according to their final values. In order to guarantee an accurate discretization, a large number of finite elements is selected in this work, taking into account the upper bounds of the batch times. Alternatively, an adaptive strategy could be implemented, by adjusting the number of finite elements for the discretization of each mathematical stage of the formulation iteratively according to their specific duration. This way, it could be possible to reduce the model size. The originated MINLP problem can be solved through several search algorithms available in literature (Grossmann, 2002). Their strategies are either based on enumeration, such as branchand-bound (B&B) methods with non-linear programming (NLP) sub-problems (Nemhauser and Wolsey, 1988), or on decomposition, such as outer approximation (OA) (Duran and Grossmann, 1986) or generalized Benders decomposition (GBD) (Geoffrion, 1972). On the one hand, enumeration strategies like the B&B nonlinear global optimization solver BARON (Sahinidis, 1996) can guarantee global solutions under general assumptions like the availability of finite lower and upper bounds on the variables and expressions in the MINLP to be solved. However, these approaches also have the disadvantage of rendering costlier computational loads because they operate through the systematic enumeration of integer combinations and have to address a NLP problem for each integer solution. On the other hand, the advantage of decomposition approaches is precisely their higher efficiency, since the problem is iteratively solved through mixed-integer linear programming (MILP) master problems with linearized equations and NLP primal

Fig. 8. Finite elements discretization for differential z(t), algebraic y(t), and control udyn (t) variable profiles, where continuity is only enforced for the first ones. Diamonds ˙ in the finite elements e ∈ {1, . . ., Ne } and boundary and represent the approximations of z(t), y(t), and udyn (t) and triangles represent the estimates of the derivatives z(t) collocation points m ∈ {0, 1, . . ., M, M + 1}.

300

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

problems with fixed integer variables. The limitation, in this case, is that decomposition approaches only guarantee local optima unless the nonlinear objective function and constraints are convex. Many of the MINLP solution algorithms have been implemented in commercial software and operate through modeling interfaces like GAMS (Rosenthal, 2012). After comparing the performance of different solvers for optimizing the numerical examples of Sections 6 and 7, the OA solver DICOPT has been selected as the most suitable alternative to be used in this study due to its efficiency. The principal disadvantage met using global enumeration solvers is the difficulty to find feasible solutions after more than 50 h of computation, in front of the attainment of optimal solutions, which are possibly local but already improve the results of non-integrated problems, in not more than 3 h with DICOPT. Global optimality cannot be guaranteed with the OA algorithm due to the existence of non-convex terms in the model, e.g. bilinear functions in mixers. Hence, it is desirable to repeat the optimization procedure using different initial feasible solutions (IFSs), with the intention of identifying accidental local optima. In this work, in order to generate IFSs effectively, the reduced problem is solved in a preliminary step by fixing feasible combinations of structural variables, namely the selection of process stages Zi , equipment units Yj , j

operating modes X i , task-unit assignments Wj,q , technologies V , j

compounds Sc , and recirculation flows Rn . Additionally, constant profiles of the control variables can be used. The combination of feasible structural decisions in the IFSs is chosen randomly, given that the starting point for the solver does not regulate the optimal configuration obtained with the full-space MINLP search. Ultimately, the objective is to start the search from different search space locations. 6. Numerical example 1: Denbigh reaction system The numerical example presented in this section addresses the introduction of a new chemical into an already existing plant by addressing the integrated batch process development problem. The optimal solution is compared to a standard predefined recipe to quantify the potential of the proposed approach. Specifically, the process consists of a competitive reaction mechanism, the Denbigh reaction system (Denbigh, 1958):

144 h. Additionally, a penalty is applied to the product shortfall. This way, the optimization problem is defined by: minimize  = Cost A + Penalty, udyn (t), ustat ,

(36)

uint , uBool j

j

where udyn (t) comprise input F1,k (t) and output F2,k (t) flow rates j

and reaction temperature k (t) in stages k ∈ {1, 2, 3} of batch units ustat

j ∈ {U1 , U2 }, is the duration tl of mathematical stages l, uint is the number of batches NBS , and uBool represent the Booleans for equipment units Yj , unit ordering Wj,q , and configurations X 1 . The first term of the objective function CostA stands for raw material expenses of the complete production campaign and is computed in base of a raw material price pˆ A of 4.8 cD /kg (for further details, see Eqs. (S16) and (S21) in Section S1 of the Supplementary Material), whereas the Penalty term represents the economic charge of unaccomplished demand, calculated considering twice the selling price pˆ S of product S, with a value of 43.1 cD /kg (see Eqs. (S17) and (S23)). The complete MLDO model for solving this example is provided in the Supplementary Material, including the specifications of the already existing equipment items, which are incorporated into the optimization problem as operating constraints. 6.2. Problem solution The MLDO problem is reformulated into a MINLP one using 32 finite elements and three collocation points in normalized Legendre roots, according to the direct-simultaneous strategy explained in Section 5. This is implemented in GAMS v.24.7.1 (Rosenthal, 2012) and solved through the decomposition-based OA solver DICOPT (Duran and Grossmann, 1986), using CONOPT 3.17A and CPLEX 12.6 to handle the NLP and MILP sub-problems respectively. In this example, four IFSs are first generated using constant profiles for the j j j control variables F1,k (t), F2,k (t), and k (t) within each mathematical stage k ∈ {1, 2, 3} of units j ∈ {U1 , U2 } and fixing the configuration Booleans X 1 = true for each of the four possible configurations ∈ {˛, ˇ, , }. Next, the complete MINLP problem is solved with the prior IFSs and with no IFS. As seen in Table 1, which summarizes the features of the resulting MINLP model implemented in GAMS, the best solution is obtained by solving the full MINLP problem initialized with IFS 4. 6.3. Results and discussion

(35) whose kinetic and thermodynamic data are available in Table S1 (see Supplementary Material). The SEN superstructure of the example corresponds to the reactor network of Fig. 1, comprising two batch reactors U1 and U2 with a capacity of 1 m3 . Since the investment in new equipment units and pipelines is initially not contemplated, the operating modes that can be selected are: operation in one single unit U1 (configuration ˛) or U2 (configuration ˇ), parallel in-phase operation in U1 and U2 (configuration ), or series operation in U1 followed by U2 (configuration ). Batch procedures in the reactors are composed of load, hold, and unload operations and the control variables are the input and output flow rates and the reaction temperature. 6.1. Optimization problem The production target is to fulfill a demand of 21 t of product S minimizing the raw material expenses within a time horizon of

The optimal solution of the integrated problem is compared with a fixed recipe that has a predefined production size of 300 kg/batch of product S and requires 70 batches and a batch cycle time of 2.06 h/batch to fulfill the demand of 21 t in its entirety within the time horizon of 144 h. The operating mode is fixed to configuration ˇ, setting Boolean decisions Xˇ1 , YU2 , and WU2 ,1 with a true value. Unit U2 is defined to operate with maximum input and output flow rates (7.7 m3 /h) during load and unload operations and maximum reaction temperature (110 ◦ C) to guarantee that the total demand is satisfied. The optimization problem is solved for three cases with different degrees of freedom to track the effect of the structural and the performance decisions: • Dynamically optimal recipe: optimization of dynamic control profiles with fixed configuration ˇ; • Structurally optimal recipe: optimization of structural decisions j j with constant profiles of control variables F1,1 (t), F2,3 (t), and j

k (t), k ∈ {1, 2, 3}, j ∈ {U1 , U2 };

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

301

Table 1 Features of the MINLP models implemented in GAMS in the Denbigh example, obtained using 32 finite elements and three collocation points in the discretization step: IFS solutions with fixed configurations X 1 = true, ∈ {˛, ˇ, , } and complete MINLP problem initialized with IFS 4.

IFS 1 (fixed ˛) IFS 2 (fixed ˇ) IFS 3 (fixed ) IFS 4 (fixed ) MINLP with IFS 4

No. equations

No. continuous variables

No. discrete variables

Non-zero elements

Non-linear terms

103,209

98,339

6

346,791

159,578

102,633

98,321

10

345,639

159,578

Solution time (s) 259 226 581 276 4294

optimizing the dynamic profiles with a predefined configuration ˇ, whereas the expenditure decreases as far as a 24% when qualitative decisions are considered as degrees of freedom with constant profiles. Going a step further, this improvement becomes a 25% when all decisions are optimized simultaneously in the optimal recipe i. The reduction of raw material expenses is related to the increase in the product selectivity, from 0.448 in the fixed recipe to 0.507, 0.593, and 0.597 in the three optimal solutions (see Table 2). Such an increase is reflected in a molar fraction of product S at the final time which is higher in the optimal solution (Fig. 10c4) than in the fixed one (Fig. 10a4). At the same time, the consumption of intermediate R rises significantly in the optimal recipe, whereas the generation of byproducts like T is considerably lower.

Fig. 9. Optimal solution in the Denbigh example: operating mode and synchronization of unit procedures.

• Optimal recipe i: simultaneous optimization of dynamic profiles and structural decisions according to the problem statement. In all cases, the number of batches and the duration of batch operations are free decision variables. 6.3.1. Results The optimal operating mode for both optimizations that include structural decisions is series configuration , as shown in Fig. 9a. This configuration activates and synchronizes the 3-stage models of units U1 and U2 , such that the unload operation of the first unit corresponds to the load operation of the second one. For the case of the optimal recipe i, the synchronization of batch operations is detailed in Fig. 9b, which illustrates the mathematical stages at Levels 1 and 0 for the production of one batch. The trajectories j j j of dynamic control variables F1,1 (t), F2,3 (t), and k (t), k ∈ {1, 2, 3}, j ∈ {U1 , U2 } are presented in Fig. 10(b1–c1, b2–c2), showing the evolution of the control variable profiles range between their lower and upper bounds. The synchronization of operations with material transference is illustrated through the equivalence between the U1 from unit U and the input flow rate F U2 to output flow rate F3,2 1 1,1 unit U2 . Additionally, Fig. 10(b3–c3, b4–c4) shows process variables like the reaction volume and molar compositions. All of them are compared to the trajectories of the fixed recipe in batch unit U2 (Fig. 10(a1–a4)). 6.3.2. Improvement extent The incentives to improve the process performance through the proposed modeling strategy are assessed by comparing the raw material savings obtained with the fixed recipe. The results, summarized in Table 2, show a 12% reduction in the consumption of A by

6.3.3. Key decision variables Essentially, the improvement of product selectivity is due to the optimization of the dynamic temperature profile and the arrangement of the two reaction units in series (configuration ). On the one hand, the temperature profile in the optimal recipe starts at the lower bound 50 ◦ C in the unit procedure of U1 (Fig. 10b1), in order to enhance the selectivity of R with respect to T by favoring reaction 1 (Ea,1 = 1000 kcal/kmol, see Table S1 in the Supplementary Material) with respect to reaction 2 (Ea,2 = 2580 kcal/kmol). The profile gradually increases to the upper bound of 80 ◦ C in U1 , thus favoring reaction 3 (Ea,3 = 1800 kcal/kmol) with respect to reaction 4 (Ea,4 = 1210 kcal/kmol) to improve the selectivity of the desired product S instead of U. In the unit procedure of U2 , the temperature profile rises even more (Fig. 10c1), reaching the upper bound of 110 ◦ C in this unit. On the other hand, production is restricted in all recipes by the maximum time horizon of 144 h (see Table 2) and by the available equipment capacities of 1 m3 in U1 and U2 (Fig. 10(a3–c3)). Therefore, series configuration favors the process performance by increasing the reaction volume. 6.3.4. Dynamic controls in operations with material transference In contrast, the use of dynamic control variables in operations with material transference barely affects the solution. The feed-forward trajectories of flow rates in the optimal recipe are lead to their upper bound most of the load/unload time intervals (Fig. 10(b2–c2)). Thus, in this example, load and unload operations tend to be as rapid as the pumping capacity permits. The dynamic feed-forward profile of reaction temperature is also kept constant in operations with material transference (Fig. 10(b1–c1)). In contrast, the dynamic control profiles are exploited in other batch operations. 6.3.5. Trade-offs among KPIs Since heating costs are not included in the objective function, the interaction between the selectivity of product S and the processing expenses does not affect the optimal temperature profiles in this example. Regardless, processing costs associated with heat consumption also diminish (see Table 2) thanks to the reduction in energy consumption in side reactions 2 and 4, which are endothermic. The rise of S selectivity also goes hand in hand with a reduction

302

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

Fig. 10. Control and process variable profiles in the fixed and optimal recipe i in the Denbigh example.

of reactant consumption and to a decrease in the total amount of reaction mixture required in the system. Consequently, the heat needed for achieving the reaction temperature is smaller and the optimal recipe benefits from a collateral reduction of heating costs of a 15%. Besides, since new investments are not considered, null amortization values are computed in all recipes. Moreover, due to a significant increase in occupation costs by using two reactors, the total profit is higher in the dynamically optimal solution with fixed configuration ˇ than in the optimal recipe i with operating mode . This way, improvements of the 11.9% and 6.5% compared

to the profit of the fixed recipe are obtained, respectively. These results reflect the duplication of labor costs, such as unit cleaning and rental. 6.4. Adaptability potential In order to evaluate the potential of the proposed strategy to enhance the batch plant flexibility and adaptability, the optimal recipe is solved according to other decision criteria, namely the profit (case ii) and profitability (case iii) maximization. These

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

303

Table 2 KPIs and individual economic weights in the Denbigh example: fixed recipe, dynamically optimal recipe in U2 , structurally optimal recipe with constant control profiles, and simultaneous optimal recipes for all degrees of freedom in the problem statement in cases: (i) raw material cost minimization, (ii) profit maximization, (iii) profitability maximization, and (iv) profit maximization with capacity expansion. Items in bold correspond to the objective function. KPI

Fixed recipe

Dynamically optimal

Structurally optimal

Optimal recipe i

Optimal recipe ii

Optimal recipe iii

Optimal recipe iv

Equipment configuration No. batches Batch size [kg/batch] Total processing time [h] Batch processing time [h/batch] Batch cycle time [h/batch] Total profit [D ] Total revenue Raw material cost Processing cost in U Occupation cost in U Amortization in U Penalty Profit per batch [D /batch] Profitability [D /h] Selectivity of S [kmol S/kmol total] Total energy consumption [kWh]

ˇ 70 300 144.03 2.06 2.06 4,764 9,046 2,250 967 1,065 0 0 68 33 0.448 38,692

ˇ 54 389 144.03 2.67 2.67 5,333 9,046 1,989 899 825 0 0 99 37 0.507 35,956

47 447 144.03 6.00 3.06 5,080 9,046 1,699 827 1,440 0 0 108 35 0.593 33,088

47 447 144.03 5.98 3.06 5,096 9,046 1,687 822 1,440 0 0 108 35 0.597 32,897

 23 913 144.03 6.26 6.26 5,841 9,046 1,696 789 720 0 0 254 41 0.595 31,558

 31 677 66.49 2.14 2.14 4,861 9,046 2,286 955 944 0 0 157 73 0.441 38,212

˛ 13 1615 144.03 11.08 11.08 5,965 9,046 1,692 729 535 125 0 459 41 0.596 29,159

two objectives additionally include product revenue (RevenueS ), processing cost (Costj,p ), occupational expenses (Costj,o ), and amortization cost (Costj,a ) of unit j, and are defined by the minimization problems: minimize  = −Profit, udyn (t), ustat , uint , uBool = −RevenueS + Penalty + Cost A +



(37)

(Cost j,p + Cost j,o + Cost j,a ),

j ∈ {U1 ,U2 }

and minimize  = −Profitability = −Profit/T total , udyn (t), ustat ,

(38)

uint , uBool considering the economic data detailed in Table S2 (see Supplementary Material). Moreover, the profit maximization problem is solved allowing the re-sizing of batch reactors to evaluate the additional process improvement that can be obtained with plant modifications (case iv). In this last case, equipment capacities larger than the original 1 m3 volume become degrees of freedom and have associated amortization cost. The resulting KPIs in these cases are presented in Table 2. The most relevant variations between optimal recipes i–iv are the selected operating modes, batch processing times, the number of lots, and temperature profiles, the last ones illustrated in Fig. 11. Additionally, the optimal recipe with equipment capacity expansion involves the re-sizing of unit U1 with a new volume of 3.5 m3 . Overall, by including occupation costs into the objective function, the series configuration is replaced by parallel one  because the simultaneous operation of U1 and U2 in configuration  involves a that fewer batches are needed – e.g. 23 lots in case ii instead of 47 in case i. This way, the sum of start-ups in each unit decreases and the occupation costs, which have a term that depends on the number of start-ups, are reduced significantly in optimal recipes ii and iii (see Table 2 and Eqs. (S2) and (S19) in the Supplementary Material). As for the processing mode with capacity expansion of U1 , the optimal configuration is ˛ with the single operation of U1 , making the smaller unit U2 available to other processes.

Fig. 11. Temperature profiles in batch units in the Denbigh example for the fixed and optimal recipes in cases: (i) raw material cost minimization, (ii) profit maximization, (iii) profitability maximization, and (iv) profit maximization with capacity expansion.

Regarding the temperature profiles, the profit maximization (case i) leads to the maximum temperature of 80 ◦ C in unit U1 and a compromise in U2 , operating at 94 ◦ C rather than at the upper bound of 110 ◦ C. In contrast, profitability maximization reaches the upper bound in unit U2 , since the processing time reduction is contemplated indirectly in the objective function (case iii); hence, the optimal solution takes into account the compromise between the economic gain and the reduction of the processing time. As a result, the processing time in optimal recipe iii achieves a lower value – almost three times less – compared to optimal recipe ii with profit maximization alone. Finally, the higher processing capacity in optimal recipe iv permits the functioning at the minimum temperature in recipes with a much longer processing time and batch size (see Table 2). Essentially, an improvement of the 23% in the total profit is noted in optimal recipe ii compared to the use of the fixed recipe. This

304

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

Fig. 12. Process stages in acrylic fiber production system: general processing scheme (solid lines) and potential processing alternatives (dashed lines).

increase is related to a reduction of raw material, occupation, and processing costs, due to a higher S selectivity, lower energy consumption, and reduction in the number of unit start-ups. Moreover, optimal recipe iii involves a rise of the 121% in the profitability objective with respect to the fixed recipe, even though this solution exhibits the highest total cost and the lowest selectivity of all the optimal solutions. Besides, the consideration of expanding the capacity of processing units in case iv hardly improves the total profit compared to the optimal solution without equipment modifications. The reason is that the original two-reactor plant was already equipped with the necessary capacity to fulfill the requested product demand, provided that the processing units are properly arranged. Nevertheless, the plant modification affects the expense distribution, dropping all the costs except for the amortization.

7. Numerical example 2: acrylic fiber production system This case study addresses the development of an industrial-size batch copolymerization process to produce an acrylic fiber with a specific composition and quality in a grassroots scenario. With this purpose, the presented modeling strategy is applied for constructing the integrated optimization model that simultaneously quantifies the interaction between process synthesis and plant allocation decisions like the selection of process stages, equipment technology, chemicals, and solvent recovery and reuse, among others. The study of the integrated problem solution is motivated by the complex reaction mechanisms and physicochemical phenomena – for example, the auto-acceleration and the Trommsdorf effect – which determine the compromise between productivity and polymer quality. Then, the choice of processing conditions, such as reactor temperature and monomer feed rate trajectories, and processing times is crucial to determine the resulting polymer properties and production yield (Giudici, 2000). Additionally, there exist trade-offs between polymerization and downstream tasks, which also contribute to the achievements of the overall economic and environmental production targets. Particularly, polymer-solvent separation can be omitted when high conversions are achieved (Bajaj et al., 1996), the magnitude of the environmental impact of organic and aqueous solvents depends on their molecular composition (Gol’dfein and Zyubin, 1990), and cleaning technologies influence the polymerization processes (Capón-García et al., 2011). The process for acrylic fiber production comprises a primary stage for bulk copolymer production and a secondary stage for transforming the copolymer into spun format, following the general process stages of the block diagram of Fig. 12 (see Section S2 of the Supplementary Material for further details). Multiple operating modes can be considered for each particular task, which determine

the required processing units and their synchronization. Moreover, the following qualitative decisions can be made: 1. Solution or suspension copolymerization in the reaction (stage 1); 2. Organic or aqueous solvent in solution copolymerization (stage 1); 3. Selection of separation (stage 2); 4. Recirculation of the solvent recovered in the separation (stage 2) to the copolymerization reaction (stage 1) in a subsequent batch; 5. Selection of washing and filtration (stage 3); 6. Recirculation of the washing water of the filtration (stage 3) to the copolymerization reaction (stage 1); 7. Operating mode in the polymer wash and filtration (stage 7), i.e. single unit F71 or series F71 followed by F72 ; 8. Operating mode in the solvent recovery (stage 8), i.e. single unit C81 or series C81 followed by C82 ; 9. Recirculation of the solvent recovered in the separation (stage 8) to the copolymerization reaction (stage 1) or the repulping (stage 4). The resulting SEN superstructure representing all the possible process stages and structural alternatives for acrylic fiber production is presented in Fig. 13. 7.1. Optimization problem The goal is to produce 5 t of acrylic fiber composed of 85% of acrylonitrile (AN) and 15% of vinyl acetate (VA) in bulk format minimizing the total operational and capital cost. Also, a maximum deviation of the 0.025 in the composition should be fulfilled to guarantee polymer quality constraints. A single-product campaign is assumed to produce the total amount of final product in a time horizon of 144 h. Particularly, the problem tackles the optimization of the batch production process and the batch plant according to the following objective function: minimize  = Cost M1 + Cost M2 + Cost I + Cost S udyn (t), ustat uint , uBool

(39)

+Cost a + Cost p + Cost waste ,

where M1 refers to AN monomer, M2 to VA monomer, I to initiator azobisisobutyronitrile (AIBN), and S to organic solvent dimethylformamide (DMF), aqueous solvent sodium thiocyanate (NaSCN(aq)), or suspension medium water. udyn (t) represents j the control profiles of input and output flow rates (Fin1,k (t), j

j

j

Fin2,k (t), Fout,k (t)) and cooling temperature ( cool,k (t)) in the solution and suspension copolymerization reactors j ∈ {R11 , R12 }, which follow three operations k ∈ {1, 2, 3} load, hold, and unload, respectively. It also represents the control profiles of the heat j (Qheat,k (t)) supplied to the evaporator j ∈ {E2 }, which comprises

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

305

Fig. 13. Superstructure of the acrylic fiber example, composed of: solution and suspension polymerization reactors R11 and R12 , evaporator E2 , washing and filtering units F3 , F5 , F71 and F72 , repulping unit R4 , spinneret S6 , distillation columns C81 and C82 , and buffer tanks T2 , T3 , T4 , T81 and T84 . Equipment elements in black represent the addressed primary polymerization process for bulk polymer production. Equipment elements in gray represent the secondary part of the process for spun production.

three operations k ∈ {1, 2, 3} load, unload of distillate, and unload of both fractions, respectively. ustat denotes the composition of raw j materials in the reactors’ input stages (cc,in1,k , c ∈ {M1 , M2 , I}, k ∈

{1, 2}) and the duration of batch operations (tl , ∀ l) and uint stands for the capacity (Sizej ) of processing units j ∈ {R11 , R12 , E2 } and the number of batches NB. Regarding the logical decisions uBool , these include the selection of separation stage (Zi , i ∈ {2}), solution or suspension technologies j in copolymerization reactors (V ,  ∈ {solu, susp}, j ∈ {R11 , R12 }), organic or aqueous solvent in solution copolymerization reacj tor (Sc , c ∈ {DMF, NaSCN(aq)}, j ∈ {R11 }), and solvent recovery and reuse after the separation (R2,1 ). These decisions are related to the selection of procesing units (Yj ), such as polymerization reactors j ∈ {R11 , R12 }, separation unit j ∈ {E2 }, and recirculation buffer tank j ∈ {T2 }. Finally, the multiple terms of the objective function  refer to: (i) the expenses of raw material, i.e. monomers M1 and M2 , initiator I, and solvent or suspension medium S (Cost M1 , Cost M2 , Cost I , Cost S ), (ii) equipment amortization (Costa ), (iii) total processing cost (Costp ), and (iv) waste disposal costs (Costwaste ). The definition of these terms and the corresponding economic parameters can be found in Section S2 of the Supplementary Material, (Eq. (S39) and Table S11). The MLDO model additionally includes disjunctive equations and logical propositions to relate the mentioned qualitative alternatives, as well as multistage models representing potential units. In particular, the process stages that have a critical impact on the cost function are the copolymerization reaction – determining the selectivity and production yield – and the separation stage – governing the potential recovery of unreacted monomer,

solvent, and suspension medium. Their process performance is represented by dynamic models with discrete events and is controlled by dynamic feed-forward control trajectories and duration of batch operations. The model and process parameters are available in Section S2 of the Supplementary Material. 7.2. Problem solution The MLDO problem is solved through a direct-simultaneous strategy. The dynamic model is fully discretized using eight finite elements and three collocation points in normalized Legendre roots. Piece-wise constant functions are used to define the profiles of control variables. Moreover, given the low influence of dynamic profiles in operations with material transference observed in example 1, the dynamic models are approximated in load and unload operations, namely stages k ∈ {1, 3} in reactors j ∈ {R11 , R12 } and evaporator j ∈ {E2 }. The resulting MINLP model is computational intense for global optimizers due to the presence of non-convex terms, even if it is moderate in size, as it can be observed in Table 3. Therefore, the problem is solved with the OA solver DICOPT in GAMS v.23.8.2 (Rosenthal, 2012), using CONOPT 3.15D and CPLEX 12.4 optimizers in the NLP and MILP sub-problems respectively. Several IFSs are used to initiate the search procedure from different initial points. Specifically, the starting points correspond to several subsystems of the integrated problem and are generated by fixing the Boolean variables that define specific structural alternatives of the process. The subsystems are determined by the following combinations of technologies and solvent:

306

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

Table 3 Features of the MINLP models implemented in GAMS in the acrylic fiber example, obtained using eight finite elements and three collocation points in the discretization step: IFS solutions with problem subsystems 1a, 1c, 2b, and 3b, and whole MINLP problem initialized with IFS 3.

IFS 1 (subsystem 1a) IFS 3 (subsystem 1c) IFS 5 (subsystem 2b) IFS 7 (subsystem 3b) MINLP with IFS 3

No. equations

No. continuous variables

No. discrete variables

Non-zero elements

Non-linear terms

Solution time (s)

16,860

8285

4

66,608

15,996

85

16,860

8285

14

66,608

15,996

338

R11 , S R11 = (1) Solution polymerization and organic solvent, i.e. Vsolu DMF true; polymerization and aqueous solvent, i.e. (2) Solution R11 , S R11 Vsolu = true; NaSCN(aq) R12 = true. (3) Suspension polymerization, i.e. Vsusp

Additionally, these alternatives are evaluated without the recirculation of distillate from the evaporator E2 (R2,1 = false in subsystems 1a, 2a, 3a) and with recirculation (R2,1 = true in 1b, 2b, 3b). In the specific case of solution polymerization with organic solvent, there exists the possibility of not having separation stage too (Z2 = false in subsystem 3c). The remaining Booleans are computed through the logical propositions of the problem. Overall, seven sets of fixed logical variables are used to generate IFSs. The computational performance to calculate some of these IFSs is summarized in Table 3.

7.3. Results and discussion 7.3.1. Trade-off among process stages The diagram of Fig. 14 presents the optimal objective function of each subsystem in front of the optimal solution of the integrated MINLP problem for comparison. The optimal solution among the three polymerization and solvent alternatives 1, 2, or 3 in the reaction stage depends on decisions associated with other process stages like the recirculation choice. Suspension polymerization is chosen when there is no recirculation of the distillate (subsystem 3a), whereas solution polymerization with an organic solvent is the best alternative otherwise (subsystem 1b). At the same time, the selection of separation stage depends on decisions made in the reaction, like the polymerization temperature and monomer dosage, which determine the achieved monomer conversion, and

Fig. 14. Total cost for each subsystem and the integrated problem with and without recirculation and avoiding separation stage (i.e. R2,1 = true, R2,1 = false, Z2 = false) in the acrylic fiber example. In black, integrated solution.

the solvent selected. The only possibility for disregarding the separation is to utilize solution polymerization with DMF achieving a minimum conversion (subsystem 1c). The proposed modeling approach integrates all these degrees of freedom and permits the simultaneous evaluation of the trade-offs between process stages to obtain the optimal solution. 7.3.2. Cost balance The KPIs and cost contributions of the subsystem solutions and the integrated problem are summarized in Table 4. Essentially, the most significant weights are equipment amortization and raw material expenses, which are considerably higher than the cooling water costs in the reaction and the heating costs in the separation. The combined expenditure of monomers M1 and M2 is the most relevant one, followed by organic and aqueous solvent costs in R11 . Therefore, the possibility of recirculating solvent in subsystems 1b and 2b encompasses important raw material savings without loss of the reaction effectiveness. Nevertheless, despite the raw material savings obtained with recirculation, the optimal solution gives priority to the total cost minimization by eliminating the equipment investment for the evaporator E2 and the recirculation storage T2 . 7.3.3. System adaptability Comparing these results with preliminary calculations with fixed batch size, it was observed how the results presented here have closer values of the objective function for the different subsystems, thanks to the consideration of the number of batches as a degree of freedom. However, this similarity does not represent homogeneous optimal solutions but the adaptability of the model to provide favorable solutions by evaluating the decisions trade-offs. Notably, previous research on flexible plant design (Moreno-Benito et al., 2014) enhanced the batch plant flexibility by combining process synthesis decisions, like dynamic control profiles, with plant allocation decisions, like batch sizing. 7.3.4. System response The optimal processing scheme corresponds to solution polymerization without separation, illustrated in Fig. 15, and is characterized by the control and processing variable profiles represented in Fig. 16. The optimal temperature trajectory in the polymerization reactor follows a monotonically increasing function, whereas the M1 dosage features a local maximum in the second half of the process that guarantees the desired polymer composition of 85% of M1 . Thanks to the dosage profile and the long task duration (i.e. 18.4 h) the process generates slight excess of monomers and fulfills the minimum monomer conversion of the 75% required for this processing structure. This solution is compared to an equivalent optimization with halved unitary amortization costs. In this case, the optimal soluR12 = tion is suspension polymerization with recirculation (i.e. Vsusp true, R2,1 = true), since the raw material savings are prioritized in front of amortization costs, which are larger than in the original economic scenario. The control and processing variable profiles, illustrated in Fig. 17, are characterized by a reduction of the

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

307

Table 4 KPIs and individual economic weights in the acrylic fiber example considering: subsystem alternatives 1a–c, 2a–b, 3a–b, and the whole MINLP problem. Items in bold correspond to the total cost minimization objective. Subsystems 1a Polymerization type Solvent type Process structure (Z2 ,R2,1 ) Equipment size [m3 ] No. batches Batch size [kg/batch] Total processing time [h] Batch proc. time [h/batch] Batch cycle time [h/batch] Total cost [D ] Amortization R11 or R12 Water consumption cost Monomer M1 cost Monomer M2 cost Initator I cost Solvent/suspension S cost Amortization E2 Energy cost Waste disposal cost Amortization T2 T4 Amortization Cost per batch [D /batch] Monomer conversion

(1,0) 1.5 8 625 119.9 18.0 16.7 46,919 12,234 113 12,048 5,082 3,355 1,192 8,368 202 1,016 – 3,309 5,865 0.44

Integrated 1b Solution DMF (1,1) 1.5 9 556 134.9 18.2 16.9 41,057 12,234 114 13,750 6,448 2,876 692 8,368 257 1,130 3,309 3,309 4,562 0.37

1c

2a

2b

(0,0) 2 6 833 89.9 18.3 18.3 38,614 14,083 113 8,568 2,267 7,449 2,824 – – – – 3,309 6,436 0.75

Solution NaSCN(aq) (1,0) (1,1) 2.5 2 6 9 833 556 89.9 86.2 19.9 11.9 18.5 10.6 44,432 42,904 15,716 14,083 106 104 8,412 16,233 1,619 8,526 2,258 1,114 732 236 10,750 9,633 268 452 1,262 1,530 – 3,309 3,309 3,309 7,405 4,767 0.83 0.29

3a Suspension – (1,0) (1,1) 1.5 1.5 17 22 294 227 144.0 144.0 11.6 8.4 10.4 7.2 41,567 42,187 12,234 12,234 372 369 10,418 13,151 1,637 2,013 1,653 2,145 99 129 8,368 8,368 655 977 1,980 2,571 – 3,309 3,309 3,309 2,445 1,918 0.71 0.56

3b Solution DMF (0,0) 2 6 833 89.9 18.4 18.4 38,620 14,083 113 8,683 2,154 7,296 2,981 – – – – 3,309 6,437 0.75

Fig. 15. Optimal structure and logical variables in the acrylic fiber example.

monomer dosage along the second stage of R12 , due to a higher input of monomer in the first load operation, coming from raw material and recirculation flow rates. There is also a higher excess of monomer M1 to guarantee the desired composition. This excess is recovered in the distillate fraction in E2 . Even though the polymerization time is shorter in this case, with a cycle time of 9.1 h, the solution exhibits a polymerization reaction significantly longer than the separation. Consequently, the consideration of more polymerization reactors with a parallel out-of-phase configuration would be challenging to reduce the cycle time in case of desiring larger production rates. To conclude, the major strength of the proposed methodology is the holistic evaluation of the decision criteria. The

optimization model considers the trade-offs between the various process stages within the complete process, as well as the interactions between synthesis and allocation degrees of freedom. Additionally, it could expedite the analysis of processing alternatives according to diverse processing, economic, or sustainable production drivers. The primary copolymerization stage for acrylic fiber production has been solved successfully, so it could be challenging to consider the secondary polymerization process to exploit the methodology. Additionally, further degrees of freedom could be incorporated into the model, like the use of several reactors and their arrangement in parallel out-of-phase configuration to allow the reduction of the cycle time and increase the production in a limited time horizon.

308

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

8. Conclusions and future directions

Fig. 16. Control and process variable profiles of the polymerization reactor R11 for the optimal solution in the acrylic fiber example: (1) flow rates/dosage and temperature and (2) molar compositions and reaction volume.

This work introduces a modeling strategy to tackle batch process development and improvement at production scale, pursuing optimal solutions that enhance process efficiency and guarantee the best utilization of available equipment and new investments. With this purpose, multiple degrees of freedom associated with process synthesis and plant allocation sub-problems are solved simultaneously, while taking into account the plant design and retrofit. The principal novelty of this strategy is the synchronization of unit procedures with dynamic profiles as a function of alternative processing schemes in an integrated optimization problem. From a modeling perspective, significant progress has been made to coordinate coexisting multistage and single-stage models through the combination of hybrid discrete/continuous models, dynamic optimization, and mixed-logic modeling, in a formulation that can be handled by current optimization tools. Overall, the resulting MLDO model manages the problem complexity using mixed-logic equations, while its reformulation into a MINLP problem enables its optimization through well-established search algorithms. The notable improvement of plant adaptability gained through the presented approach justifies and stimulates the research in this direction. The complex mathematical implications and the risk of obtaining mathematically intractable problems can be addressed with appropriated solution strategies like the use of initial feasible solutions and the definition of proper variable bounds. The promising results of the numerical example show improvements between the 23% and the 121% in the objective function compared to the use of fixed recipes in all the optimization scenarios i–iv of the Denbigh example. Also, the rapid evaluation of diverse objective functions proves the potential of the proposed strategy to adapt master recipes according to in-time needs, as they could be changes in economic scenarios or decision criteria, and to analyze the system response in front of uncertain parameters.

Fig. 17. Control and process variable profiles of the polymerization reactor R12 and the separator E2 for the optimal solution with reduced amortization costs in the acrylic fiber example: (a1) flow rates/dosage and temperature, (b1) vapor flow, and (a2–b2) molar compositions and volumes.

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

This strategy also enables the quantification of interactions between synthesis and allocation sub-problems. A larger influence on structural decisions with constant set-points optimization is identified in the Denbigh example in comparison to the optimization of dynamic profiles with fixed equipment configuration. Both cases render improvements over the fixed recipe of the 24% and the 12% respectively. In contrast, the consideration of dynamic flow rate profiles in the operations with material transference has a minor impact, since most of the corresponding stages are characterized by constant profiles in the optimal solutions of the numerical cases tested. Finally, the acrylic fiber example shows the trade-offs among decisions in neighboring process stages. Therein, the influence of the distillate recirculation on reaction decisions like the selection of the polymerization technology and solvent is illustrated. Likewise, reaction decisions like the temperature and dosage profiles, raw monomer concentration, and processing times determine the possibility of omitting separation stage if high monomer conversions are achieved. The large advantage of the presented integrated strategy is its potential to assess all processing alternatives simultaneously and avoid methodologies based on the identification of challenging solutions and their independent solution, with the probability of missing optimal solutions. Future work can be focused on the incorporation of additional degrees of freedom in the problem, like out-of-phase parallel operation and multi-product/multi-purpose production campaigns. The modular representation of unit procedures in the presented strategy should allow the direct extension of the formulation to achieve this purpose. Equally, the refinement of solution strategies is highly relevant to propitiate the use of global optimization solvers and to face problems with a larger number of process stages and products. In this regard, the study of adaptive approaches that adjust the number of finite elements in the discretization of the dynamic models according to the selected batch times is suggested. The solution of this problem using hybrid approaches that combine deterministic solvers with evolutionary algorithms is another promising area of research.

Acknowledgements The authors would like to thank the financial support received from the Spanish Ministry of Economy and Competitiveness and the European Regional Development Fund (research project SIGERA DPI2012-37154-C02-01 and FPU Research Fellowship). Marta Moreno would also like to thank Prof. Luis Puigjaner for the fruitful discussions held during the development of this work.

309

Table 5 Basic logical operators expressed using logical and algebraic equations (Raman and Grossmann, 1991), where Ad , d ∈ {1, . . ., N} are Boolean variables and ad , d ∈ {1, . . ., N} are the equivalent binaries. Logical operator

Logical expression

Algebraic equation

Conjunction “and” (∧) Disjunction “or” (∨) “exclusive or” () Implication (→)

A1 ∧ A2 ∧ · · · ∧ AN A1 ∨ A2 ∨ · · · ∨ AN A1  A2  · · ·  AN A1 →A2 or ¬A1 ∨ A2 (A1 → A2 ) ∧ (A2 → A1 ) or(¬A1 ∨ A2 ) ∧ (¬A2 ∨ A1 )

a1 ≥ 1, a2 ≥ 1, . . ., aN ≥ 1 a1 + a2 + · · · + aN ≥ 1 a1 + a2 + · · · + aN = 1 1 − a1 + a2 ≥ 1 or a1 − a2 ≤ 0 a1 − a2 ≤ 0, a2 − a1 ≤ 0 or a1 = a2

Equivalence (⇔)

3. Distributing the “or” operator (∨) over the “and” operator (∧) recursively by using the following equivalence: (A1 ∧ A2 ) ∨ A3 ⇔ (A1 ∨ A3 ) ∧ (A2 ∨ A3 ).

(A.4)

Table 5 summarizes the relation between logical and algebraic expressions, according to Raman and Grossmann (1991). Appendix B. Application of orthogonal collocation on finite elements Orthogonal collocation on finite elements (Cuthrell and Biegler, 1989) is used to perform the full discretization of the timecontinuous model. The collocation points used in the numerical examples are calculated with a shifted Legendre polynomial P() of order M = 3, which is defined as: P() = 20  3 − 30  2 + 12  − 1.

(B.1)

The roots of the above polynomial provide the normalized locations  m = {0.1127, 0.5000, 0.8873} of the collocation points m ∈ {1, . . ., M}, in addition to the boundary points located at  0 = 0 and  M+1 = 1. Then, Lagrange polynomials are used for the approximation of the state and control variables. In particular, the polynomials of the differential z(t), algebraic y(t), and control udyn (t) variables, ˙ have the following form: and derivatives z(t) z(t) =

M 

ze,m m (t),

m=0 M

y(t) =



˙ z(t) =

M 

ze,m ˙ m (t),

m=0

ye,m m (t),

m=1

udyn (t)

=

M 

(B.2) dyn ue,m

m (t),

m=1 dyn

Appendix A. Reformulation of logical propositions into linear constraints The operations for formulating the conjunctive normal form (CNF) of a logical proposition were formalized by Clocksin and Mellish (1981). Considering Boolean variables Ad , d ∈ {1, 2, 3}, these operations are: 1. Replacing the implication by its equivalent disjunction: A1 ⇒ A2 ⇔ ¬A1 ∨ A2 ;

(A.1)

2. Moving the negation inward by applying DeMorgan’s Theorem: ¬(A1 ∧ A2 ) ⇔ ¬A1 ∨ ¬A2 ,

(A.2)

¬(A1 ∨ A2 ) ⇔ ¬A1 ∧ ¬A2 ,

(A.3)

where ze,m , ye,m , and ue,m represent the approximations of the differential, algebraic, and control variables in finite elements e ∈ {1, . . ., Ne } and collocation or boundary points m ∈ {0, 1, . . ., M, M + 1}, M = 3 with normalized locations  m = {0, 0.1127, 0.5000, 0.8873, 1}. m (t) and m (t) represent the basis functions defined by:

 (t − te,m ) (t − te,m ) , m (t)= . (te,m − te,m ) (te,m − te,m ) m = 0, m = 1, / m / m m = m =

m (t)=

M 

M

(B.3)

The polynomials of differential variables z(t) are one order higher Mth than the polynomials of algebraic y(t) and control variables udyn (t), with order (M − 1)th. The purpose is to define boundary conditions of the approximation of differential variables ze,0 in each finite element e, thanks to the symmetry properties of the method. This is necessary to enforce their continuity from previous element e − 1 in time te . Thus, accurate initial conditions ze,0 = z(te−1 ) = ze−1,M+1 can be calculated through the evaluation of

310

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311

the first polynomial in Eq. (B.2) at the final point  M+1 = 1 of each finite element: ze,0 =

M 

ze−1,m m (M+1 = 1),

∀ e ∈ {2, . . ., Ne }.

(B.4)

m=0

Finally, control variables are approximated to a piecewise condyn dyn stant (PWC) function by establishing the relation ue,m = ue . Appendix C. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.compchemeng. 2016.07.030. References Ahmad, B.S., 1997. Synthesis of Batch Processes with Integrated Solvent Recovery (Ph.D. thesis). Massachusetts Institute of Technology. Ali, S.A., 1999. Synthesis of Batch Processing Schemes for the Production of Pharmaceuticals and Specialty Chemicals (Ph.D. thesis). Massachusetts Institute of Technology. Allgor, R.J., Barton, P.I., 1999. Mixed-integer dynamic optimization I: problem formulation. Comput. Chem. Eng. 23 (4–5), 567–584. Allgor, R., Evans, L., Barton, P., 1999. Screening models for batch process development: part I. Design targets for reaction/distillation networks. Chem. Eng. Sci. 54 (19), 4145–4164. Bajaj, P., Sen, K., Bahrami, S.H., 1996. Solution polymerization of acrylonitrile with vinyl acids in dimethylformamide. J. Appl. Polym. Sci. 59 (10), 1539–1550. Balas, E., 1985. Disjunctive programming and a hierarchy of relaxations for discrete optimization problems. SIAM J. Algebraic Discrete Methods 6 (3), 466–486. Banimostafa, A., Papadokonstantakis, S., Hungerbühler, K., 2012. Retrofit design of a pharmaceutical batch process improving green process chemistry & engineering principles. In: Karimi, I.A., Srinivasan, R. (Eds.), 11th International Symposium on Process Systems Engineering. Computer Aided Chemical Engineering, vol. 31. Elsevier, pp. 1120–1124. Bansal, V., Perkins, J.D., Pistikopoulos, E.N., 2002. A case study in simultaneous design and control using rigorous, mixed-integer dynamic optimization models. Ind. Eng. Chem. Res. 41 (4), 760–778. Barbosa-Póvoa, A.P., 2007. A critical review on the design and retrofit of batch plants. Comput. Chem. Eng. 31 (7), 833–855. Barton, P.I., Pantelides, C.C., 1994. Modeling of combined discrete/continuous processes. AIChE J. 40 (6), 966–979. Barton, P.I., Allgor, R.J., Feehery, W.F., Galan, S., 1998. Dynamic optimization in a discontinuous world. Ind. Eng. Chem. Res. 37 (3), 966–981. Barton, P.I., Ahmad, B.S., Cheong, W., Tolsma, J., 1999. Synthesis of batch processes with integrated solvent recovery. In: Sikdar, S.K., Diwekar, U. (Eds.), In: Tools and Methods for Pollution Prevention, vol. 62. Kluwer Academic Publishers, Dordrecht, pp. 205–231. Beaumont, N., 1990. An algorithm for disjunctive programs. Eur. J. Oper. Res. 48 (3), 362–371. Bemporad, A., Morari, M., 1999. Control of systems integrating logic, dynamics, and constraints. Automatica 35 (3), 407–427. Bhatia, T., Biegler, L.T., 1996. Dynamic optimization in the design and scheduling of multiproduct batch plants. Ind. Eng. Chem. Res. 35 (7), 2234–2246. Binder, T., Blank, L., Bock, H.G., Bulirsch, R., Dahmen, W., Diehl, M., Kronseder, T., Marquardt, W., Schlöder, J.P., von Stryk, O., 2001. Introduction to Model Based Optimization of Chemical Processes on Moving Horizons. Springer-Verlag, Berlin, pp. 295–339. Bryson, A.E., Ho, Y.-C., 1975. Applied Optimal Control: Optimization, Estimation and Control. John Wiley & Sons, Inc. Bumann, A.A., Papadokonstantakis, S., Fischer, U., Hungerbühler, K., 2011. Investigating the use of path flow indicators as optimization drivers in batch process retrofitting. Comput. Chem. Eng. 35 (12), 2767–2785. ˜ A., Puigjaner, L., 2011. Multiobjective Capón-García, E., Bojarski, A.D., Espuna, optimization of multiproduct batch plants scheduling under environmental and economic concerns. AIChE J. 57 (10), 2766–2782. ˜ A., 2013. Integrating process Capón-García, E., Guillén-Gosálbez, G., Espuna, dynamics within batch process scheduling via mixed-integer dynamic optimization. Chem. Eng. Sci. 102, 139–150. Carvalho, A., Matos, H.A., Gani, R., 2009. Design of batch operations: systematic methodology for generation and analysis of sustainable alternatives. Comput. Chem. Eng. 33 (12), 2075–2090. Cavin, L., Fischer, U., Glover, F., Hungerbühler, K., 2004. Multi-objective process design in multi-purpose batch plants using a Tabu Search optimization algorithm. Comput. Chem. Eng. 28 (4), 459–478. Chakraborty, A., Linninger, A.A., 2002. Plant-wide waste management. 1. Synthesis and multiobjective design. Ind. Eng. Chem. Res. 41 (18), 4591–4604. Charalambides, M.S., Shah, N., Pantelides, C.C., 1995. Synthesis of batch reaction/distillation processes using detailed dynamic models. Comput. Chem. Eng. 19 (S1), 167–174.

Chu, Y., You, F., 2014. Integrated scheduling and dynamic optimization of network batch processes. In: American Control Conference (ACC), pp. 5024–5029, http://dx.doi.org/10.1109/ACC.2014.6858593, ISSN 0743-1619. Clocksin, W.F., Mellish, C.S., 1981. Programming in Prolog. Springer-Verlag, New York. Corsano, G., Aguirre, P.A., Iribarren, O.A., Montagna, J.M., 2004. Batch fermentation networks model for optimal synthesis, design, and operation. Ind. Eng. Chem. Res. 43 (15), 4211–4219. Corsano, G., Montagna, J.M., Iribarren, O.A., Aguirre, P.A., 2007. Heuristic method for the optimal synthesis and design of batch plants considering mixed product campaigns. Ind. Eng. Chem. Res. 46 (9), 2769–2780. Cuthrell, J., Biegler, L., 1989. Simultaneous optimization and solution methods for batch reactor control profiles. Comput. Chem. Eng. 13 (1–2), 49–62. Denbigh, K., 1958. Optimum temperature sequences in reactors. Chem. Eng. Sci. 8 (1-2), 125–132. Duran, M., Grossmann, I.E., 1986. An outer-approximation algorithm for a class of mixed-integer nonlinear programs. Math. Program. 36, 307–339. Ferrer-Nadal, S., Capón-García, E., Mendez, C.A., Puigjaner, L., 2008. Material transfer operations in batch scheduling. A critical modeling issue. Ind. Eng. Chem. Res. 47 (20), 7721–7732. Floudas, C.A., Grossmann, I.E., 1994. Algorithmic approaches to process synthesis: logic and global optimization. In: Biegler, L.T., Doherty, M.F. (Eds.), Foundations of Computer-Aided Process Design. , pp. 198–221. Frankl, K., Beenken, J., Marquardt, W., 2012. Integrated scheduling and control of continuous-time blending processes. In: Karimi, I.A., Srinivasan, R. (Eds.), 11th International Symposium on Process Systems Engineering. Computer Aided Chemical Engineering, vol. 31. Elsevier, pp. 1090–1094. Furman, K.C., Jia, Z., Ierapetritou, M.G., 2007. A robust event-based continuous time formulation for tank transfer scheduling. Ind. Eng. Chem. Res. 46 (26), 9126–9136. Gani, R., Papaeconomou, I., 2005. Conceptual, design and synthesis of batch processes. In: Chemical Industries. CRC Press, pp. 43–82. Geoffrion, A.M., 1972. Generalized Benders decomposition. J. Optim. Theory Appl. 10, 237–260. Giudici, R., 2000. Polymerization reaction engineering: a personal overview of the state of art. Latin Am. Appl. Res. 30, 351–356. Gol’dfein, M.D., Zyubin, B.A., 1990. Kinetics and mechanism of the processes of preparing fibre-forming polymers based on acrylonitrile – review. Polym. Sci. U.S.S.R. 32 (11), 2145–2166. Grossmann, I.E., Guillén-Gosálbez, G., 2010. Scope for the application of mathematical programming techniques in the synthesis and planning of sustainable processes. Comput. Chem. Eng. 34 (9), 1365–1376. Grossmann, I., Westerberg, A., Biegler, L., 1987. Retrofit design of processes. In: Reklaitis, G., Spriggs, H. (Eds.), Proceedings of the First International Conference on Foundations of Computer Aided Process Operations. , ISBN 0 444 98925 0, pp. 403–442. Grossmann, I.E., 1985. Mixed-integer programming approach for the synthesis of integrated process flowsheets. Comput. Chem. Eng. 9 (5), 463–482. Grossmann, I.E., 2002. Review of nonlinear mixed-integer and disjunctive programming techniques. Optim. Eng. 3, 227–252. Halim, I., Carvalho, A., Srinivasan, R., Matos, H.A., Gani, R., 2011. A combined heuristic and indicator-based methodology for design of sustainable chemical process plants. Comput. Chem. Eng. 35 (8), 1343–1358. Iribarren, O., Malone, M., Salomone, H., 1994. A heuristic approach for the design of hybrid batch-continuous processes. Chem. Eng. Res. Des. 72 (A3), 295–306. Iribarren, O.A., Montagna, J.M., Vecchietti, A.R., Andrews, B., Asenjo, J.A., Pinto, J.M., 2004. Optimal process synthesis for the production of multiple recombinant proteins. Biotechnol. Prog. 20 (4), 1032–1043. Iribarren, O.A., 1985. Batch Chemical Processes Design (Ph.D. thesis). University of Massachusetts. Jain, S., Kim, J.-K., Smith, R., 2013. Process synthesis of batch distillation systems. Ind. Eng. Chem. Res. 52 (24), 8272–8288. Kondili, E., Pantelides, C., Sargent, R., 1993. A general algorithm for short-term scheduling of batch operations – 1. MILP formulation. Comput. Chem. Eng. 17 (2), 211–227, ISSN 0098-1354. Linninger, A.A., Chakraborty, A., 1999. Synthesis and optimization of waste treatment flowsheets. Comput. Chem. Eng. 23 (10), 1415–1425. Linninger, A.A., Stephanopoulos, E., Ali, S.A., Han, C., Stephanopoulos, G., 1995. Generation and assessment of batch processes with ecological considerations. Comput. Chem. Eng. 19 (Suppl. 1), 7–13. Luus, R., 1990. Application of dynamic programming to high-dimensional nonlinear optimal control problems. Int. J. Control 52 (1), 239–250. ˜ A., Puigjaner, L., 2014. Flexible batch process and plant Moreno-Benito, M., Espuna, design using mixed-logic dynamic optimization: single-product plants. Ind. Eng. Chem. Res. ˜ A., Puigjaner, L., 2015. Integrated Moreno-Benito, M., Dombayci, C., Espuna, process and plant design optimisation of industrial scale batch systems: addressing the inherent dynamics through stochastic and hybrid approaches. Chem. Eng. Trans. 45, 1789–1794. Mosat, A., Cavin, L., Fischer, U., Hungerbühler, K., 2008. Multiobjective optimization of multipurpose batch plants using superequipment class concept. Comput. Chem. Eng. 32 (3), 512–529. Nemhauser, G.L., Wolsey, L.A., 1988. Integer and Combinatorial Optimization. John Wiley & Sons, Inc. Neuman, C., Sen, A., 1973. A suboptimal control algorithm for constrained problems using cubic splines. Automatica 9 (5), 601–613.

M. Moreno-Benito et al. / Computers and Chemical Engineering 94 (2016) 287–311 Nie, Y., Biegler, L.T., Wassick, J.M., 2012. Integrated scheduling and dynamic optimization of batch processes using state equipment networks. AIChE J. 58 (11), 3416–3432. Nie, Y., 2014. Integration of Scheduling and Dynamic Optimization: Computational Strategies and Industrial Applications (Ph.D. thesis). Carnegie Mellon University. Oldenburg, J., Marquardt, W., 2008. Disjunctive modeling for optimal control of hybrid systems. Comput. Chem. Eng. 32 (10), 2346–2364, ISSN 0098-1354. Oldenburg, J., Marquardt, W., Heinz, D., Leineweber, D.B., 2003. Mixed-logic dynamic optimization applied to batch distillation process design. AIChE J. 49 (11), 2900–2917. Pantelides, C.C., 1994. Unified frameworks for the optimal process planning and scheduling. In: Proceedings of the 2nd Conference on the Foundations of Computer Aided Operations, p. 253. Papaeconomou, I., 2005. Integration of Synthesis and Operational Design of Batch Processes (Ph.D. thesis). Department of Chemical Engineering, Technical University of Denmark. Prata, A., Oldenburg, J., Kroll, A., Marquardt, W., 2008. Integrated scheduling and dynamic optimization of grade transitions for a continuous polymerization reactor. Comput. Chem. Eng. 32 (3), 463–476. Raman, R., Grossmann, I.E., 1991. Relation between MILP modelling and logical inference for chemical process synthesis. Comput. Chem. Eng. 15 (2), 73–84. Reklaitis, G., 1989. Progress and issues in computer-aided batch process design. In: 3rd International Conference on Foundations of Computer-Aided Process Design, Snowmass, Colorado, July 09–14, pp. 241–275. Rippin, D.W.T., 1983. Simulation of single- and multiproduct batch chemical plants for optimal design and operation. Comput. Chem. Eng. 7 (3), 137–156. Robinson, J.D., Loonkar, Y.R., 1972. Minimizing capital investment for multi-product batch-plants. Process Technol. 17 (11), 861. Rosenthal, R.E., 2012. GAMS: A User’s Guide. Scientific Press. Sahinidis, N.V., 1996. BARON: a general purpose global optimization software package. J. Glob. Optim. 8 (2), 201–205.

311

Sargent, R.W.H., Sullivan, G.R., 1978. The development of an efficient optimal control package. In: Stoer, J. (Ed.), Optimization Techniques. Lecture Notes in Control and Information Sciences, vol. 7. Springer-Verlag, Berlin Heidelberg, pp. 158–168. Sharif, M., Shah, N., Pantelides, C.C., 1999. Design of integrated batch processes with discrete and continuous equipment sizes. Comput. Chem. Eng. 23 (Suppl.), S117–S120. Sharif, M., Samsatli, N.J., Shah, N., 2000. Abstract design in the development of pharmaceutical processes. In: Pierucci, S. (Ed.), European Symposium on Computer Aided Process Engineering – 10. Computer Aided Chemical Engineering, vol. 8. Elsevier, pp. 685–690. Simon, L.L., Osterwalder, N., Fischer, U., Hungerbühler, K., 2008. Systematic retrofit method for chemical batch processes using indicators, heuristics, and process models. Ind. Eng. Chem. Res. 47 (1), 66–80. Smith, E., Pantelides, C., 1995. Design of reaction/separation networks using detailed models. Comput. Chem. Eng. 19 (Suppl. 1), 83–88. Srinivasan, B., Palanki, S., Bonvin, D., 2003. Dynamic optimization of batch processes – I. Characterization of the nominal solution. Comput. Chem. Eng. 27 (1), 1–26. Stein, O., Oldenburg, J., Marquardt, W., 2004. Continuous reformulations of discrete-continuous optimization problems. Comput. Chem. Eng. 28 (10), 1951–1966. Stephanopoulos, G., Reklaitis, G.V., 2011. Process systems engineering: from Solvay to modern bio- and nanotechnology. A history of development, successes and prospects for the future. Chem. Eng. Sci. 66 (19), 4272–4306. Türkay, M., Grossmann, I.E., 1998. Tight mixed-integer optimization models for the solution of linear and nonlinear systems of disjunctive equations. Comput. Chem. Eng. 22 (9), 1229–1239. Tsang, T., Himmelblau, D., Edgar, T., 1975. Optimal control via collocation and nonlinear-programming. Int. J. Control 21 (5), 763–768. Vassiliadis, V.S., Sargent, R.W.H., Pantelides, C.C., 1994. Solution of a class of multistage dynamic optimization problems. 1. Problems without path constraints. Ind. Eng. Chem. Res. 33 (9), 2111–2122.