Optimal Control of Methanol Synthesis Fixed-Bed Reactor - American

Aug 9, 2013 - ABSTRACT: The possibility to apply the nonlinear model predictive control (NMPC) to the fixed-bed tubular reactor for methanol synthesis...
0 downloads 0 Views 2MB Size
Article pubs.acs.org/IECR

Optimal Control of Methanol Synthesis Fixed-Bed Reactor Flavio Manenti* and Giulia Bozzano Dipartimento di Chimica, Materiali e Ingegneria Chimica “Giulio Natta” Politecnico di Milano, Piazza Leonardo da Vinci 32, 20133 Milano, Italy ABSTRACT: The possibility to apply the nonlinear model predictive control (NMPC) to the fixed-bed tubular reactor for methanol synthesis is broached. The distributed nature of the methanol synthesis reactor implies a simultaneous management of temporal and spatial dynamics to control the conventional physical time-dependent parameters (temperature, flows, etc.) but also the axial migration of the temperature hot-spot typical of methanol synthesis. The approach can be easily extended to other systems characterized by exothermic reactions. Special attention is given to the prompt response of NMPC to ensure the online effectiveness. Theoretical and practical information is provided.

1. INTRODUCTION

2. EXTENDED NMPC TECHNIQUES Whenever new techniques or methodologies are proven to be effective in their original field, the interest starts increasing in other areas that are relatively close. This situation is what happened to NMPC in the past decade, with different modifications and advances so as to extend its uses to other research topics. For instance, the extension of NMPC from continuous to batch processes has been one of the recent advances in the field.11 Pistikopoulos’s idea of a multiparametric NMPC12,13 that leads to an explicit representation of NMPC has extended the application domain to microchips. Another variant is proposed by Lima and co-workers,14 where the usual differential-algebraic constraints are replaced by a Takagi− Sugeno15 fuzzy model to characterize in an easy way strong nonlinear systems such as polymer/copolymer reactors. The key to the success of the method is the relevant reduction of the computational complexity. An interesting approach to extend the NMPC application field is to combine a series of models and use each of them according to some predefined situations and operating conditions. For example, Dougherty and Cooper16 introduced a multimodel adaptive control strategy for multivariable NMPC, and Guiamba and Mulholland17 developed and implemented an adaptive linear dynamic matrix control algorithm to control an integrated pump-tank system consisting of two input and two output variables. Cameron and co-workers proposed a multimodel approach18,19 to achieve robust control and global stability. Dones et al.9 proposed and validated in an industrial plant a self-adaptive NMPC approach, which automatically selects the best level of detail of the mathematical model. A NMPC using a simplified model is much more flexible (and feasible) than a NMPC using a detailed model. Enabling the algorithm to adjust the model to the changes in the dynamics makes it possible to boost and optimize the performances of the NMPC. A conscious approach has been proposed,20 where the NMPC is able to self-select the best set of algorithms according to the problem, the operating conditions, and the control

Nonlinear model predictive control (NMPC) is still mainly perceived as an academic technique, although it is one of the most appealing approaches to manage process units and plant sections. Benefits of NMPC have been demonstrated in several contributions, especially applied to distillation columns and polymer reactors, where the strong nonlinear behavior allows one to highlight its potential.1−5 Nevertheless, many issues are still open for its wide application such as (i) the computation time that exceeds the real-time requirements, (ii) the insufficient accuracy or loss of stability, and (iii) the limitations in achieving some performance specifications. Such issues are dramatically amplified for industrial applications, which are often characterized by a distributed nature. Actually, distributed systems usually lead to the solution of partial differential and differentialalgebraic equations, and this is the main reason for which they have received much less attention than lumped systems. Specifically, distributed systems require additional and deeper studies on the stability of the NMPC, the implementation of very high performing solvers to ensure the online feasibility, and especially the need to extend to the concept of NMPC to more than a single coordinate. Whereas the former points have been already investigated in the recent literature,6−10 the extension of NMPC to the multidimensional problem of controlling the process along more coordinates is not yet well-defined. Within this context, the present activity is aimed at proposing the theoretical and practical guidelines to extend the traditional concept of NMPC and apply it for the first time to the tubular fixed-bed boiler-reactor for methanol synthesis, demonstrating its industrial effectiveness. The paper is organized as follows. Section 2 gives an overview on the recent advances to extend the NMPC applicability. Nonlinear techniques for optimal control are described in section 3, including receding horizon methodology and practical rules for tuning. Essential modifications to the classical NMPC are explained in section 4. The application of the extended NMPC to the methanol synthesis reactor is discussed in section 5. Numerical results and NMPC tuning are given in section 6. © 2013 American Chemical Society

Received: Revised: Accepted: Published: 13079

May 12, 2013 August 4, 2013 August 9, 2013 August 9, 2013 dx.doi.org/10.1021/ie401511e | Ind. Eng. Chem. Res. 2013, 52, 13079−13091

Industrial & Engineering Chemistry Research

Article k + hp

scenario. In practice, the NMPC is able to select the best differential solver to speed up calculations and preserve the accuracy of predictions, but also to switch from efficient to robust optimizers and vice versa whenever long computational times or hard minimization problems are automatically detected. At last, Christofides et al.21 extended the application of NMPC to several distributed systems, where the convolution model is generally constituted by partial differential equations. Analogously, the present paper proposes to extend Christofides et al.’s approach to a Lurgi-type fixed-bed boiler-reactor for methanol synthesis22,23 focusing the attention on the distributed nature of the selected system as well as on the promptness of the response due to the industrial feasibility of such an approach.

min

+

F(x̅ (τ ), u̅ (τ )) dτ

l=k



ω2, l(u̅ (l i) − u̅ (l i − 1))2 }

(4)

where xj̅ − xset j is the deviation between the jth controlled variable and its set point, u̅l − utar l is the deviation between the lth (i−1) manipulated variable and its steady-state target, u̅(i) is the l − uj̅ incremental variation between ith and (i − 1)th time intervals of the lth manipulated variable, coefficients of the diagonal semipositive definite matrices ω are the weighting factors, HP is the prediction horizon, HC is the control horizon, and min and max superscripts indicate lower and upper bounds, respectively, for both manipulated and controlled variables. Set points and targets are assigned by the user, and they usually stem steadystate economic optimizations. The ultimate goal of NMPC is to find the vector u̅ that minimizes the gap between controlled variables and their set points and between manipulated variables and their targets, by satisfying the constraints to which the objective function is subject. MPC based on linear/linearized models had great success as testified by the spreading of dedicated commercial tools.24,25 Nevertheless, since the 1990s, interest has started moving toward nonlinear models (leading to the nonlinear model predictive control, NMPC) mainly pushed by the fact that processes underwent more and more tight environmental regulations and market requirements and process specifications became narrower and narrower. As explained elsewhere,26,27 the NMPC does not refer to a specific control strategy, rather to an explicit use of convolution models to minimize a proper objective function, which usually consists of a set of quadratic terms. Some authors are used to inserting even economic terms within the NMPC objective function, obtaining a kind of integration between the optimal control (quadratic nature) and the dynamic optimization (economic nature). Actually, NMPC is a single level of a more structured process control and optimization hierarchy,28−33 usually named the integration pyramid, where every hierarchical level is characterized by specific targets and constraints.34,35 A controller that handles nonlinearities is more effective in achieving such constraints for chemical processes.26,36,37 The development of first-principles mathematical models is crucial as it is practically infeasible (or anyhow very complex) to extend process identification from a linear to a nonlinear case. The apparent advantage is that the model is closer to the real plant behavior, and process parameters can be evaluated without claiming massive plant data. The main disadvantage of NMPC is its numerical complexity, leading to the need of very high performing solvers for differential systems and robust (but even efficient) optimizers to ensure online feasibility (for more detail, see also ref 38) and the need to focus attention on the routine stability so as to prevent any kind of ringing and divergence problems (for more details, see also refs 6 and 8).

(2)

subject to ẋ = Ax̅ + Bu̅ y = Cx

ω1, l(u̅ l − ultar)2

subject to: ⎧ x min ≤ x ≤ x max j ̅j ⎪ j ⎪ ⎪ ulmin ≤ u̅ l ≤ ulmax ⎨ ⎪ Δu (i),min ≤ Δu(i) = u(i) − u(i − 1) ≤ Δu (i),max l ̅l ̅l ̅l ⎪ l ⎪ ⎩Convolution Model

(1)

t + HP



l=k

where

∫t

k + hp − 1

ωj(x̅ j − xjset)2 +

k + hp − 1

Due to more and more frequent market dynamics and competitiveness, conventional control systems consisting of proportional−integral−derivative loops are progressively less effective at handling grade/load changes of chemical processes. In addition, if we consider the pressing need to reduce losses during process transients for environmental laws, safety, and material exploitation, it is necessary to adopt more flexible and predictive control techniques. Generally, the optimization problem of NMPC is formulated as follows:

J(x(t ), u̅ ( ·); HC , HP): =



u̅ (k)....u̅ (k + hc − 1) j=k+1

3. NONLINEAR MODEL PREDICTIVE CONTROL

min J(x(t ), u̅ ( ·); HC , HP) u(·)

{

(3)

HC and HP are control and prediction horizons of the NMPC, respectively. Their role is explained later. The notations u̅ and x ̅ are necessary to distinguish virtual variables that are inside the NMPC routine from real variables (x,u). Such a separation is necessary as each kind of model is unable to reproduce the plant behavior exactly, and a deviation occurs during the calculations of NMPC. It is however worth underlining that very detailed convolution models usually provide a more accurate prediction, but their solution requires higher computational efforts (not always23); a good compromise between the model detail and the computational effort is required for the implementation of NMPC. Actually, a reasonably detailed model should require relatively small computational efforts, and its deviation with respect to the real plant should be negligible according to the control discretization adopted on the plant. The generic objective function is usually stated as follows: 13080

dx.doi.org/10.1021/ie401511e | Ind. Eng. Chem. Res. 2013, 52, 13079−13091

Industrial & Engineering Chemistry Research

Article

Figure 1. Receding horizon and online application: qualitative scheme.

Thus, while the NMPC must be invoked at each time interval (Figure 1), a series of associated devices has to be promptly solved too for its effectiveness, such as data reconciliation to purify raw data39−41 from error in measures as well as gross errors usually affecting plant measures.42,43 Specifically, on the basis of reconciled measures acquired at time t0, NMPC predicts the future behavior of the system on the selected prediction horizon HP and determines the optimal series of vectors of manipulated variables U̅ (t0) = [u1̅ (t0), ..., uN̅ (t0)] by minimizing an opportune objective function (eq 4). Only the first vector u̅1 is implemented in the real plant so as to handle disturbances, model−plant deviations, and uncertainty in general, and the procedure is iterated at the next time interval t1 according to a methodology that undergoes what is called a receding horizon for NMPC and dynamic optimization problems and a rolling horizon for scheduling and planning problems.44 3.1. Tuning. Control and prediction horizons HC and HP are two important parameters of NMPC, and their adequate selection is essential to guarantee the NMPC effectiveness and stability.45,46 Dubljevic and co-workers have pointed out that the issue of stability is usually addressed by means of imposing suitable penalties and constraints on the state at the end of the optimization horizon.8 In the present work, we do not impose stability constraints in the optimization problem, but we have increased the HP length and appropriately selected the sampling time of NMPC, which cannot be smaller than the plant sampling time to avoid the partial optimization of the same sampling time. Moreover, a specific antiringing term (see also eq 4) has been adopted to further prevent the algotithm instability. It is worth remarking that such a nonoptimal condition would be iteratively implemented at each calculation of NMPC. Next, a reasonable HP should be on the same order of magnitude of the plant characteristic time. According to the selection of HP and HC, the system is either stable with smooth control actions or close to instability or unstable.47 More and more frequently, industrial best practice proposes the relation: HP = 5HC

as good compromise between computational efforts and prediction capacity, whereas it should be better to tune NMPC according to the specific process unit. In any case, relation 5 can be adopted as a first guess to tune the system. 3.2. Solution Strategies. Two solutions are available for NMPC, and they differ in the discretization of constraints within the optimization problem.48,49 In the sequential approach, the dynamic model is separate from the optimizer. Simulation and optimization are performed sequentially, and the convolution model is numerically integrated along HP at each iteration of the optimizer. This approach is also called the partially parametrized method since manipulated variables only are discretized and the feasible path method since a feasible solution of the dynamic system is needed at each iteration. In the simultaneous approach, the dynamic model is fully discretized as per manipulated variables, and it is directly integrated in the optimization problem. This approach is sometimes called the f ully parametrized method. It is also named the infeasible path method since it is not necessary to have feasible solutions of the dynamic model (except for the last iteration) as the numerical integration and the optimization simultaneously converge. The simultaneous approach leads to a sparse structured system that can be integrated by specific solvers,50 but the dimensions of the resulting nonlinear programming problem are unavoidably enlarged. On the other hand, the sequential approach keeps the size of the optimal problem small, and it is easier to implement, but it generally requires higher computational effort due to the larger number of iterations.

4. EXTENSION OF NMPC ALGORITHM TO CONTROL DISTRIBUTED SYSTEMS More than other scientific areas, process engineers are used to linking process control to time. This is because process control is almost always adopted to regulate unsteady operating conditions of static structures such as process units, reactors, distillation columns, and so on but also since most of the dynamic models show an evolution along time or, if not, because the relevant coordinate along which the system is to be controlled is time. Thousands of time-evolving applications are proposed in the

(5) 13081

dx.doi.org/10.1021/ie401511e | Ind. Eng. Chem. Res. 2013, 52, 13079−13091

Industrial & Engineering Chemistry Research

Article

equations (PDEs and PDAEs) to implement NMPC on plug flow reactors and subsequently on other units.53−55 At the same time, other research groups56,57 worked on computational effort reduction by proposing methods to significantly reduce problem size and therefore the amount of degrees of freedom of the resulting nonlinear programming problem to speed up NMPC solution by reasonably preserving accuracy and optimality of the procedure. Merging these advances, an extended NMPC acting on more coordinates is no longer infeasible, and the methanol synthesis reactor is the ideal validation field with the tasks of (i) optimizing the methanol production, (ii) controlling the reactor temperature, and (iii) monitoring the hot-spot axial migration. 5.1. Methanol Synthesis. The modeling issues for the methanol synthesis reactor are already reported elsewhere with deep detail22,23 on the dynamic characterization and special focus on the numerical methods. We report only a brief introduction of the system for self-consistency of the paper. Methanol is still largely produced by syngas (CO and H2 mixture) obtained via steam reforming operations.58,59 It consists of four main phases: (i) feed treatment (purification by sulfur and water impurities), (ii) reforming to convert methane and steam in syngas, CH4 + H2O ⇄ CO + 3H2 and CO + H2O ⇄ CO2 + 3H2, (iii) methanol synthesis to convert syngas to methanol, and (iv) product purification and storage. It is well-known that the crucial point of the overall process modeling and simulation is the synthesis reactor. According to the current industrial trend, the shell and tube reactor-boiler60 is selected. The tube side is filled by the catalyst, whereas the shell side is fed by water to regulate the reactor temperature and to generate high-pressure steam. A shell and tube reactor favors process controllability61 and improves the process efficiency with respect to traditional intercooled reactors. The overall reaction from carbon monoxide to methanol CO + 2H2 ⇄ CH3OH is an exothermal process with ΔH298K = −90.55 [kJ/mol] and equilibrium ΔG0|408K = 0 and favored at low temperatures to the detriment of the reaction rate. To get a reasonable industrial conversion, it is necessary to operate with a specific catalyst at high pressure because of the decrease in the total amount of moles throughout the synthesis. Methanol synthesis from carbon dioxide and the reverse watergas shift must be both considered:

literature for grade changes, market uncertainties, regulation and servomechanism problems, plant startup, and emergency shutdown, but it is also worth underlining the possibility to extend process control to any kind of coordinate where the system could evolve. For instance, our previous works22,23 have shown that the modeling of the industrial methanol synthesis fixed-bed reactor may involve several coordinates where the system could evolve. In other words, the multidimensional optimization problem of NMPC is furthermore extended from the single coordinate of time-dependent phenomena to m coordinates along which the system could evolve. The real control issue in this case is the need to find as many manipulated variables as the number of controlled variables for each coordinate along which they could evolve. As a consequence, vectors of controlled and manipulated variables that are usually adopted in formulating NMPC problems should be replaced by matrices: x̅ ∈ RnCV → X̅ ∈ RnCV × m (6) u̅ ∈ RnMV → U̅ ∈ RnMV × m where nCV and nMV are the number of controlled and manipulated variables, respectively. Hence, the optimization is not performed anymore to find out the optimal vector of vectors {u̅1,u̅2,....,uN̅ } of degrees of freedom, where N is the number of sampling times included in the HC, but the optimal vector of matrices {U̅ 1, U̅ 2, ...., U̅ N}:

min

{U1, U2,...., UN }

J(X̅ , {U̅ 1, U̅ 2, ...., U̅ N }; h C, hP)

(7)

Note that it is not necessary to adopt the same control and prediction horizons for all coordinates of process dynamics, but the vectors hC and hp, which are both of dimension m (number of coordinates), could be separately assigned according to the specific coordinate. In addition, it is worth underlining that if any controlled variable is of interest only in correspondence with the end of a specific coordinate (in other words, its profile along such a coordinate is not of interest), no manipulated variables are needed to regulate it along the same coordinate on the condition that an ad hoc Mayer term appears in the objective function of the extended NMPC. The resulting nonlinear programming is furthermore complex: whereas the traditional NMPC has nMV·N degrees of freedom, the extended problem could have up to nMV· N·m degrees of freedom. To better explain and implement the extended NMPC, the Lurgi methanol synthesis reactor is adopted as a case study.

CO2 + 3H 2 ⇄ CH3OH + H 2O ΔH298K = −49.43 [kJ/mol]

CO2 + H 2 ⇄ CO + H 2O

(8)

ΔH298K = + 41.12 [kJ/mol] (9)

5. APPLICATION TO THE METHANOL SYNTHESIS REACTOR Computational effort made the NMPC useless for some applications of chemical engineering for a long time. It is not a coincidence that optimal control is not applied to systems usually well-characterized by partial differential equation systems, such as crystallizers, fixed-bed reactors, and solid-state polymerizers to quote a few, if not in a simplified way. Recent advances in modern computers and programming (i.e., multiprocessor machines, object-oriented programming, and parallel computing38) have pushed these applications to look forward the NMPC as an appealing and interesting methodology also in these cases. Christofides and co-workers51,52 started replacing ordinary differential and differential-algebraic equation (ODE and DAE) systems that are usually behind NMPC by distributed models consisting of partial differential and partial differential-algebraic

In addition, some side reactions could take place: CO + 3H2 ⇄ CH4 + H2O (methanation), nCO + 2nH2 → CnH2n+1OH + (n− 1)H2O, and 2CH3OH → CH3OCH3 + H2Oeven though they are inhibited by the high catalyst selectivity under the nominal operating conditions.62 It is worth underlining that the methanation reaction, which occurs at temperatures higher than 300 °C and is catalyzed by Fe-based complexes, might bring to the reactor instability and runaway for its relevant exothermicity ΔH298K = −209.45 [kJ/mol]: this justifies the internal copper covering for the tube side. Even though bases for an optimal schedule of programmed maintenances are provided in this paper, the life-cycle of the commercial catalyst is still not considered here as its deactivation phenomena are on the order of years, and it is significantly far from NMPC and dynamic optimization time-scales of interest. Nominal operating con13082

dx.doi.org/10.1021/ie401511e | Ind. Eng. Chem. Res. 2013, 52, 13079−13091

Industrial & Engineering Chemistry Research

Article

ditions are 490−515 K and 77 bar; the molar fraction entering the reactor is reported in Table 1.

Water: εb ρgas

Table 1. Feed Molar Composition components

molar fraction

CO CO2 H2 H2O CH3OH N2 CH4

0.046 0.094 0.659 0.0004 0.005 0.093 0.1026

dωH2O, n

M ωH2O, n − ωH2O, n − 1 Δz A int ωH2O, n + 1 − 2ωH2O, n + ωH2O, n − 1 + +ρg (Δz)2

=−

dt

+ PM H2OrH2O, n

Carbon monoxide: εb ρgas

dωCO, n dt

M ωCO, n − ωCO, n − 1 + +ρg A int Δz ωCO, n + 1 − 2ωCO, n + ωCO, n − 1

=−

5.2. First-Principles Dynamic Modeling. Looking at the methanol system described above, the corresponding compound−element matrix of Table 2 for methanol synthesis has

(Δz)2 + PMCO( −rCH3OH, n + rH2O, n)

Table 2. Compound−Element Matrix C O H

CO

CO2

H2

CH3OH

H2O

1 1 0

1 2 0

0 0 2

1 1 4

0 1 2

εb ρgas

dωCO, n dt

M ωCO, n − ωCO2, n − 1 + +ρg Δz A int ωCO2, n + 1 − 2ωCO2, n + ωCO2, n − 1

=−

(Δz)2 + PMCO2( −rH2O, n)

r1

p1

0 1 1 1

0 0 0 5

0 0 10 10

(13)

Hydrogen:

Table 3. Parameters of the Objective Function F1 for the OneCoordinate Analysis q1

(12)

Carbon dioxide:

rank equal to Table 3.23 It means that two reactions are enough to characterize the overall system, and only water and methanol mass balances can be proposed in an independent form.

open-loop without penalty optimization of methanol production anti-ringing

(11)

εb ρgas

dω H 2 , n dt

M ω H 2, n − ω H 2, n − 1 + +ρg Δz A int ωH2, n + 1 − 2ω H2 , n + ωH2, n − 1

=−

(Δz)2 + PM H2( −2rCH3OH, n − rH2O, n)

(14)

Energy balance: The diffusion term leads to a second-order, parabolic PDE system, the features of which allow one to prevent possible nonphysical oscillations generated by a first-order system. The central form leads to a sparse, diagonal-block band matrix that can be very efficiently solved by specific classes of the BzzMath library,63,64 able to automatically exploit the system sparsity and structure, leading to reasonable computational times in spite of the system complexity. Furthermore, comparing a first order PDE system (without diffusion terms) to the second-order parabolic one, it is possible to underline even a significant reduction in the computational time, as several numerical instabilities are properly removed. CPU times are halved in practice. Thus, the dynamic model is provided in mass balances and the finite difference form: Methanol: εb ρgas

dωCH3OH, n dt

mix

=k

cat

Tn + 1 − 2Tn + Tn − 1



(Δz)2

U (Tshell − Tn) + [rCH3OH( −ΔHr1) A int

+ rH2O( −ΔHr2)]

(15)

where

M ωCH3OH, n − ωCH3OH, n − 1 =− Δz A int

rCH3OH = ρcat (1 − εb)η1(r1 + r3)

(16)

rH2O = ρcat (1 − εb)η2(r2 + r3)

(17)

and65−67 ⎡ fCH OH ⎤ k1K CO⎢fCO fH1.5 − 0.53 ⎥ f H KP1 ⎦ 2 ⎣ 2 r1 = ⎡ 0.5 ⎛ K H2O ⎞ ⎤ (1 + K COfCO + K CO2fCO )⎢f H + ⎜ 0.5 ⎟fH O ⎥ 2 ⎣ 2 ⎝ K H2 ⎠ 2 ⎦

+ +ρg ωCH3OH, n + 1 − 2ωCH3OH, n + ωCH3OH, n − 1 (Δz)2 + PMCH3OHrCH3OH, n

dTn dt T − Tn − 1 M − cp n mix Δz A int

[εb ρgas c p + (1 − εb)ρcat c p ]

(10)

(18) 13083

dx.doi.org/10.1021/ie401511e | Ind. Eng. Chem. Res. 2013, 52, 13079−13091

Industrial & Engineering Chemistry Research

Article

⎡ fCH OH fH O ⎤ k 3K CO2⎢fCO fH1.5 − 31.5 2 ⎥ fH KP3 ⎦ ⎣ 2 2 2 r2 = ⎡ 0.5 ⎛ K H2O ⎞ ⎤ (1 + K COfCO + K CO2fCO )⎢f H + ⎜ 0.5 ⎟fH O ⎥ 2 ⎣ 2 ⎝ K H2 ⎠ 2 ⎦

This motivates the wide margin TFIXED‑BED ≤ 540 K assumed to operate methanol synthesis reactors, even if both catalyst sintering and the methanation reaction take place at temperatures higher than 570 K. A smaller margin could be used if a more accurate way to estimate hot-spot temperature and, at the same time, to control its position were available. This is possible by exploiting the intrinsic features of the convolution model: with the system being governed by PDEs that characterize it along time and spatial coordinates, it is possible to dimensionally extend the optimization problem from the common time coordinate only to a multicoordinate problem involving even spatial objectives. Thus, rather than controlling only the hot-spot temperature, while the system undergoes disturbances and load changes, it is also possible to monitor and control the hot-spot axial position. At the same time, beyond these control objectives, there is also the possibility to insert economical terms: methanol production can be maximized while the problem of monitoring and control of the temperature hotspot and its position is solved. The matrix X̅ ∈ R(nCV)×(m=2) of the controlled variables in the extended NMPC form is

(19) fCO fH O ⎤ ⎡ k 2K CO2⎢fCO fH − K 2 ⎥ ⎣ 2 2 P2 ⎦ r3 = ⎡ 0.5 ⎛ K H2O ⎞ ⎤ (1 + K COfCO + K CO2fCO )⎢f H + ⎜ 0.5 ⎟fH O ⎥ 2 ⎣ 2 ⎝ K H2 ⎠ 2 ⎦

(20)

5.3. NMPC Implementation. Another important point is to have a model able to exploit at best the instrumentation installed by the field and, if needed, to infer some measures that are neither measurable for any reason nor so accurate.68,69 The temperature profile of methanol synthesis reactors is measured by multipoint thermocouples placed in some tubes of the tube bundle. These tools provide a discrete temperature profile, usually every 0.5 m or more. This inevitably means that the hot-spot position can easily fall far from temperature measurement points, and values coming from the field may provide a picture of the plant that is completely different from the current operating condition. In other words, because of their discrete nature, measures acquired by multipoint thermocouples may be unable to provide the real hot-spot value as the same hot-spot position probably easily falls within two points of the thermocouple (Figure 2), by losing the information about the maximum temperature of the catalytic bed and very probably bringing an earlier catalyst deactivation due to catalyst sintering.

⎡y 0 ⎤ CH3OH ⎥ X̅ = ⎢ ⎢⎣ Thotspot zhotspot ⎥⎦

(21)

where rows are controlled variables and columns are dimensions along which the system evolves. It is worth noting that the methanol fraction yCH3OH is controlled only along the time coordinate, and, since there is no interest in the methanol fraction profile along the reactor but only in its final value exiting the reactor, the methanol production can be optimized throughout the operations without the need of any specific manipulated variable. To do so, it is enough to insert an opportune term in the objective function. Thus, a constrained, nonlinear, multidimensional optimization problem is obtained. The shell temperature Tshell and the recycle flow rate Mrec are the two degrees of freedom belonging to the matrix U̅ ∈ R(nMV)×(m = 2) adopted to control both hot-spot temperature and position and, at the same time, to optimize the methanol production. ⎡ Tshell 0 ⎤ ⎥ U̅ = ⎢ ⎢⎣ M rec 0 ⎥⎦

(22)

Note that neither along which coordinates the manipulated variables evolve nor that the number of manipulated variables along a certain coordinate corresponds to the amount of controlled variables for the same coordinate is important. The

Figure 2. Multipoint thermocouple measures (light dots) and temperature axial profile. The hot-spot (dark dot) is often inaccurately measured.

Figure 3. Adaptive grid to improve the accuracy to infer hot-spot temperature and position. 13084

dx.doi.org/10.1021/ie401511e | Ind. Eng. Chem. Res. 2013, 52, 13079−13091

Industrial & Engineering Chemistry Research

Article

only important thing is that manipulated variables are properly selected by a parametric sensitivity in order to promptly act on the system. In the selected case in the study, the second column of the matrix U̅ is identically zero as no manipulated variables operate along the spatial coordinates, but only along the time coordinate, even though the control objective is to match the hot-spot temperature along time and its position along the reactor axis. For the sake of completeness, we mention that the second column of the matrix U̅ can be filled in systems adopting intermediate fresh feed flow rates and/or intermediate cooling, but this is not the case for the methanol synthesis reactor.

6. OPTIMAL CONTROL The resulting dynamic model of the methanol synthesis reactor involves PDEs, and it is solved by using the finite difference method as above. It is then discretized along the spatial coordinate obtaining a series of ODE systems. Since the model can be used to infer the hot-spot measurement, it is suitable to make thicker and thicker points of the finite difference grid in order to have a more accurate estimation of the hot-spot temperature and its position (Figure 3). Unfortunately, computational times are unavoidably increased by introducing more grid points. An effective solution is to adopt adaptive grids:70,71 smaller elements are used in the first part of the reactor (0−2 m) and larger elements in the remaining part so as to balance the total number of elements. Thus, a grid with 100 points is adopted, and the resulting constrained nonlinear programming involves 700 equations as constraints. The prediction horizon is HP = 450 s. The control horizon is HC = 90 s, and the sampling time is 30 s leading to N = 3 intervals within the HC. The one-coordinate case (time only), where the optimal control manages the reactor temperature by acting on the shell side temperature, involves nMV·N = 1 × 3 = 3 degrees of freedom. In the extended case, where the shell temperature and the recycle flow rate manage the hot-spot temperature and its position, the degrees of freedom are nMV·N = 1 × 3 = 3. Note that m = 1 since the matrix U̅ has no manipulated variables along the spatial coordinate in the case study. The optimization of methanol production is assumed as test-case in the following paragraphs where both the problems of one coordinate and extended NMPC are solved by progressively increasing their complexity in the objective function: • open-loop analysis • closed-loop without penalty terms • optimization of the methanol production • optimization of the methanol production with antiringing terms 6.1. One-Coordinate (Time) Analysis. The first problem deals with the maximization of the methanol molar fraction acting on the shell-side temperature. Figure 4 shows the improvement in the methanol fraction exiting the reactor of the closed loop simulation with respect to the open loop simulation. The objective function F1 for this analysis is as follows, and its corresponding parameters are reported in Table 3:

Figure 4. Methanol molar fraction. Optimized variable.

As mentioned, the selected horizons for the NMPC are HP = 450 s (process dynamics are on the order of some minutes) and HC = 90 s, selected as a good compromise between control and computational efforts. The sampling time is 30 s, which allows one to ensure the online feasibility of the NMPC (actually, the worst numerical conditions indentified in the tests lead to almost 26 s of computational effort). The methanol fraction outflowing the reactor is not a real controlled variable, even though it has all the features of controlled variables, being within the objective function and significantly sensitive to manipulated variables. Figure 4 shows that all objective functions considered are able to move the system from the open-loop condition to an optimized condition. By considering the necessity to define an upper bound for hot-spot temperature (540 K), the objective function must include a penalty term to account for the possibility that the system is forced to overcome such a bound. Figures 5 and 6 both

Figure 5. Hot-spot temperature.

show the importance of the penalty term since the objective function without it takes the system to overcome the bound: Thotspot = 542 K, leading the system toward methanation and catalyst sintering conditions. At last, the antiringing term is aimed at preventing enduring oscillations in the shell temperature, which could easily lead to relevant thermo-mechanical stress in the reactor (or in the steam drum) and to sudden pressure drops. It is worth underlining that this manipulated variable is intentionally expressed in degrees Celsius rather than Kelvin to prevent any possibility that the normalized term is negligible with respect to other terms of the objective function (eq 23). Manipulating only the shell temperature, it is not possible to

⎡ ⎛ T i − T i − 1 ⎞2 2 ⎢ F1(x̅ , u̅ ) = ∑ q1(yC H OH, i − 0.1) + r1⎜ shell i − 1 shell ⎟ 3 ⎢ Tshell ⎝ ⎠ i=1 ⎣ N

⎛ Thotspot, i − 540 ⎞⎤ + p1 ⎜ ⎟⎥ 540 ⎝ ⎠⎥⎦

(23) 13085

dx.doi.org/10.1021/ie401511e | Ind. Eng. Chem. Res. 2013, 52, 13079−13091

Industrial & Engineering Chemistry Research

Article

the objective function F2 is formulated as follows, and its parameters are reported in Table 4: N

F2(X̅ , U̅ ) =

∑ [q1(yCH OH,i − 0.1)2 i=1

3

⎛ Tshell, i − Tshell, i − 1 ⎞2 ⎟⎟ + q2(zhotspot, i − 1)2 + r1⎜⎜ Tshell, i − 1 ⎝ ⎠ ⎛ M rec, i − M rec, i − 1 ⎞2 ⎟⎟ + r2⎜⎜ M rec, i − 1 ⎝ ⎠ ⎛ Thotspot, i − 540 ⎞ + p1 ⎜ ⎟] 540 ⎝ ⎠

(24)

Figure 6. Shell temperature. Manipulated variable.

Table 4. Parameters of the Objective Function F2 for the Multicoordinate Analysis

simultaneously control the temperature hot-spot and its position (Figure 7). Thus, having the need to operate close to a desired

open-loop without penalty optimization of methanol production anti-ringing

q1

q2

r1

r2

p1

0 1 1 1

0 5 5 5

0 0 0 5

0 0 0 2

0 0 10 10

In this case, by having introduced the recycle mass flow rate Mrec, even the hot-spot position is controlled. To implement the present control system, the following two simplifications were adopted: only the fresh feed coming from the upstream reformer is subject to composition disturbances, and an ideal flash separating the liquid compounds is assumed downstream of the methanol synthesis reactor. The scheme of the plant section is reported in Figure 8. A step disturbance at the inlet syngas composition was implemented to check the closed-loop effectiveness of the system, and one of them is provided in Table 5. The selected perturbation is richer in carbon monoxide in spite of hydrogen leading to a temperature increase within the reactor. Figures 9−13 underline the importance of all terms involved in the objective function F2; actually, the penalty term ensures that the upper bound of 540 K of the hot-spot temperature is satisfied (Figure 11), while the antiringing terms regarding manipulated variables U̅ ensure control stability. The tuning of parameters r1 and r2 assigned in these simulations are discussed later. In Figure 9, the function without the penalty term seems the best performing one in terms of methanol molar fraction. Unfortunately, no penalty terms in the objective function bring the reactor temperature to significantly overcome the upper bound of 540 K as shown in Figure 11. Moreover, the methanol yield is only apparently increased since Figure 12 clearly shows that the corresponding flow rate entering the reactor is practically halved with respect to the other cases. At last, having supposed that thermocouples measure the reactor temperature every 0.5 m, Figure 13 shows that the objective function including penalty and antiringing terms is able to properly control even the hotspot position by keeping it in correspondence with 1 m. 6.3. Tuning and Simplification of the Objective Function F2. An important aspect is to study the possibility of simplifying the objective function to reduce the computational effort. In other words, it is important to check if all the terms included in eq 24 are relevant for the optimization problem. Specifically, it is important to see if the antiringing term for the shell side temperature r1((Tshell,i − Tshell,i−1)/(Tshell,i−1))2 is still

Figure 7. Hot-spot position.

temperature, the hot-spot position can freely move across the reactor axis. As a result, even though short- and medium-term operations are well-accomplished by controlling the reactor temperature, the long-term operations are surely nonoptimal as uncertain information about the hot-spot position are available and the catalyst could easily deactivate randomly in accordance with the operating conditions. 6.2. Extended NMPC (Time and Spatial) Analysis. As mentioned above, hot-spot position is a crucial parameter to be monitored and controlled in order to prolong the catalyst lifecycle but also to have higher accuracy and reliability in the reactor temperature acquired by the field as one can force the hot-spot to be very close to a specific point of multipoint thermocouples. In practice, there is not a desired position for the temperature hotspot, but it is important to continuously change the hot-spot position to preserve the catalyst integrity. Consequently, the wide margin of reactor temperature upper bound could be raised to Tupperbound > 540 K to increase the methanol yield. In the hotspot following, the upper bound of hot-spot temperature is maintained at 540 K for the sake of clarity. Unfortunately, the hot-spot axial position cannot be controlled by acting on the shell temperature if we are already using it for controlling the hot-spot temperature. An additional manipulated variable and a new objective function are needed. Assuming for the sake of simplicity that HP and HC are constant for both time (450 s and 90 s, respectively) and spatial (7 m and 2 m, respectively) coordinates, 13086

dx.doi.org/10.1021/ie401511e | Ind. Eng. Chem. Res. 2013, 52, 13079−13091

Industrial & Engineering Chemistry Research

Article

Figure 8. Plant section scheme.

Table 5. Step Disturbance at the Inlet Syngas Composition unperturbed stream perturbed stream

yCO

yCO2

yH2

yH2O

yCH3OH

yCH4

0.2052 0.2252

0.084 0.084

0.6853 0.6653

0.005 0.005

0. 0.

0.0205 0.0205

Figure 9. Methanol molar fraction. Optimized variable.

Figure 11. Hot-spot temperature.

Figure 10. Shell side temperature.

Figure 12. Flow rate entering the reactor (Msyngas + Mrec).

Analogously to the one-coordinate case, it is necessary to tune all parameters involved in F2: in particular, a specific closed-loop sensitivity analysis was performed to estimate r2 related to the gas recycle flow rate, which significantly influences the value of the objective function and final results of optimization. Figures 15−19 seem to show contrasting trends with respect to the value of r2: the methanol molar fraction exiting the reactor increases by decreasing the fresh process flow rate, whereas it operates a dilution of methanol molar compositions (recycle stream is richer in inerts). These opposite aspects can compensate each other as the massive methanol production is substantially the same and equal to 3.55 kg/s for each value of r2. On this subject,

relevant in the extended NMPC. Such a test is easily accomplished by imposing r1 = 0. Even though many variables do not change significantly with or without this term, Figure 14 shows the relevant deviation in the shell temperature (manipulated variable). Such a deviation is relevant when some dynamics are occurring, but also for the new steady state achieved by the system (Tshell = 528 K with r1 = 5 and Tshell = 524 K with r1 = 0). Moreover, it is worth noting that the shell side temperature is frequently changed by the control scheme when the antiringing term is removed by the objective function. All of these considerations force us to preserve the antiringing term for the shell side temperature. 13087

dx.doi.org/10.1021/ie401511e | Ind. Eng. Chem. Res. 2013, 52, 13079−13091

Industrial & Engineering Chemistry Research

Article

Figure 17. Tuning of parameter r2. Hot-spot temperature. Figure 13. Hot-spot position.

Figure 18. Tuning of parameter r2. Inlet reactor stream (Msyngas + Mrec). Figure 14. Tuning of parameter r1. Shell side temperature (manipulated variable no. 1).

Figure 19. Tuning of parameter r2. Hot-spot position. Figure 15. Tuning of parameter r2. Outlet methanol molar fraction.

the criteria to tune r2 must refer to (1) the ability to preserve and control the hot-spot position, (2) the ability to satisfy the upper bound of hot-spot temperature, and (3) the ability to smooth the control action by reducing oscillations of the system. Figure 15 shows that the parameter r2 is crucial to monitoring the hot-spot position. The tuned parameters are reported in Table 6. Figures 16 and 17 show that a reduced antiringing term r2 related to the recycle stream (i.e., r2 = 1) increases the relative importance of this manipulated variable with respect to the shell side temperature. This implies a reduced manipulation of the steam drum pressure that could be relevant in terms of equipment mechanical stress. Also the hot-spot value is Table 6. Parameters F2

Figure 16. Tuning of parameter r2. Shell side temperature. F2 13088

q1

q2

r1

r2

p1

1

5

5

2

10

dx.doi.org/10.1021/ie401511e | Ind. Eng. Chem. Res. 2013, 52, 13079−13091

Industrial & Engineering Chemistry Research

Article

remarkably lower for the same reason. In spite of a(n) (approximately) constant methanol production, the significant reduction of the recycle stream imposed by a reduced term r2 (i .e., r2 = 1 in Figure 18) implies an increased stream rich in hydrogen to the flare system. This is important since, although a higher methanol concentration is outflowing the reactor (Figure 19), the hydrogen losses to the flare increase to overcome the process benefits. The optimal value of r2 value is consequently equal to 2. Greater values of it lead to a reduction in the recycle manipulation, but one of the main characteristics of the control system fails to reach its objectives: Figure 19 clearly shows that the set point of the hot-spot position is not anymore satisfied.

■ ■

7. CONCLUSIONS AND FUTURE DEVELOPMENT This work has defined, implemented, tuned, and exploited an extension of nonlinear model predictive control (NMPC) by applying it to an industrial case study and demonstrating its online feasibility and tangible benefits. The traditional theory of NMPC is extended from the classical case that accounts for time evolution to the time and spatial coordinates for the optimization of distributed problems. In practice, it has been possible to optimize the methanol production yield by simultaneously controlling the temperature hot-spot value and position leading to direct benefits in the catalyst life-cycle. The potentialities of the extended NMPC still have to be fully investigated, especially looking forward to the possibility to manage the catalyst deactivation by controlling the hot-spot position while preserving the optimal methanol synthesis condition. This could lead to new designs of the tubular reactor and novel maintenance strategies.





ρcat = density of the catalytic pellet [kg/m3] ρgas = density of the gas phase [kg/m3] ωi = mass fraction of the component i in the gas phase [−]

NONLINEAR MODEL PREDICTIVE CONTROL HC = control horizon HP = prediction horizon m = number of coordinates N = number of sampling times included in the HC set = set point of controlled variable tar = target of manipulated variable u,̅ U̅ = vector and matrix of manipulated variables x,̅ X̅ = vector and matrix of controlled variables CATALYST AND REACTOR SPECIFICATIONS ρcat = 1770 dp = 5.47 × 10−3 av = 627 ε/τ = 0.123 number of tubes = 2962 tube length = 7 outer diameter of the tube = 3.8 × 10−2 thickness of the tube = 1.98 × 10−3 tube conductibility = 19

MATHEMATICAL RELATIONSHIPS (PROPERTY: METHOD) gas mixture viscosity: Lucas72−74 gas mixture density: RKS75 gas mixture conductibility: Stiel and Thodos76 external heat transfer coefficient: Mostinsky77 internal heat transfer coefficient: Chilton−Colburn analogy, Gupta and Thodos relationship78,79 gas mixture diffusivity: Fuller, Schettler, and Giddings80

AUTHOR INFORMATION

Corresponding Author

*Phone: +39 02 2399 3273. Fax: +39 02 7063 8173. E-mail: fl[email protected].



Notes

The authors declare no competing financial interest.



REFERENCES

(1) Benamor, S.; Doyle, F. J.; McFarlane, R. Polymer grade transition control using advanced real-time optimization software. J. Process Control 2004, 14 (4), 349−364. (2) Cervantes, A. M.; Tonelli, S.; Brandolin, A.; Bandoni, J. A.; Biegler, L. T. Large-scale dynamic optimization for grade transitions in a low density polyethylene plant. Comput. Chem. Eng. 2002, 26 (2), 227−237. (3) Chatzidoukas, C.; Perkins, J. D.; Pistikopoulos, E. N.; Kiparissides, C. Optimal grade transition and selection of closed-loop controllers in a gas-phase olefin polymerization fluidized bed reactor. Chem. Eng. Sci. 2003, 58 (16), 3643−3658. (4) Lima, N. M. N.; Filho, R. M.; Embiruqu, M.; Maciel, M. R. W. A cognitive approach to develop dynamic models: Application to polymerization systems. J. Appl. Polym. Sci. 2007, 106 (2), 981−992. (5) Manenti, F.; Rovaglio, M. Integrated multilevel optimization in large-scale polyethylene terephthalate plants. Ind. Eng. Chem. Res. 2008, 47 (1), 92−104. (6) Dubljevic, S.; Christofides, P. D. Predictive control of parabolic PDEs with boundary control actuation. Chem. Eng. Sci. 2006, 61 (18), 6239−6248. (7) Dubljevic, S.; El-Farra, N. H.; Mhaskar, P.; Christofides, P. D. Predictive control of parabolic PDEs with state and control constraints. Int. J. Robust Nonlinear Control 2006, 16 (16), 749−772. (8) Dubljevic, S.; Mhaskar, P.; El-Farra, N. H.; Christofides, P. D. Predictive control of transport-reaction processes. Comput. Chem. Eng. 2005, 29 (11−12), 2335−2345. (9) Dones, I.; Manenti, F.; Preisig, H. A.; Buzzi-Ferraris, G. Nonlinear Model Predictive Control: a Self-Adaptive Approach. Ind. Eng. Chem. Res. 2010, 49 (10), 4782−4791.

SYMBOLS: METHANOL SYNTHESIS REACTOR MODEL Aint = internal area of the tube [m2] a = specific surface area of the catalytic pellet [m2/m3] cpmix = specific heat of gas at constant pressure [J/kg·K] f x = x component activity [−] Kx/kx = x component kinetic constants (Graaf kinetics) [−] M = mass flow rate [kg/s·tube] MW = molar weight of component [kg/kmol] T = temperature of the gas phase [K] Tshell = temperature of the shell side of reactor [K] U = overall heat transfer coefficient [W/m2·K] z = axial coordinate [m] t = time [s] rWGS = reaction rate of water gas shift [mol/s·kg] rCO→CH3OH = reaction rate of methanol from carbon monoxide [mol/s·kg] rCO2→CH3OH = reaction rate of methanol from carbon dioxide [mol/s·kg] + = mass diffusivity [m2/s] k = thermal conductivity [W/m·K] εb = void fraction of catalytic bed [−] ΔHreac = enthalpy of jth reaction [J/mol] j ηj = efficiency of reaction jth [−] 13089

dx.doi.org/10.1021/ie401511e | Ind. Eng. Chem. Res. 2013, 52, 13079−13091

Industrial & Engineering Chemistry Research

Article

(10) Manenti, F.; Dones, I.; Buzzi-Ferraris, G.; Preisig, H. A. Efficient Numerical Solver for Partially Structured Differential and Algebraic Equation Systems. Ind. Eng. Chem. Res. 2009, 48 (22), 9979−9984. (11) Abel, O.; Marquardt, W. Scenario-integrated on-line optimization of batch reactors. J. Process Control 2003, 13 (8), 703−715. (12) Pistikopoulos, E. N. Perspectives in Multiparametric Programming and Explicit Model Predictive Control. AIChE J. 2009, 55 (8), 1918−1925. (13) Pistikopoulos, S. Explicit MPC: Achieving Model Predictive Control on a Chip. Chem. Eng. Prog. 2009, 105 (8), 16−16. (14) Lima, N. M. N.; Manenti, F.; Maciel Filho, R.; Embiruçu, M.; Wolf Maciel, M. R. Fuzzy Model-Based Predictive Hybrid Control of Polymerization Processes. Ind. Eng. Chem. Res. 2009, 48 (18), 8542− 8550. (15) Takagi, T.; Sugeno, M. Fuzzy Identification of Systems and its Applications to Modeling and Control. IEEE Trans. Syst., Man, Cybern. 1985, 15, 116. (16) Dougherty, D.; Cooper, D. A Practical Multiple Model Adaptive Strategy for Multivariable Model Predictive Control. Control Eng. Pract. 2003, 11, 649. (17) Guiamba, I. R. F.; Mulholland, M. Adaptive Linear Dynamic Matrix Control Applied to an Integrating Process. Comput. Chem. Eng. 2004, 28, 2621. (18) McGahey, S. L.; Cameron, I. T. A multi-model repository with manipulation and analysis tools. Comput. Chem. Eng. 2007, 31 (8), 919− 930. (19) Wang, F. Y.; Bahri, P.; Lee, P. L.; Cameron, I. T. A multiple model, state feedback strategy for robust control of non-linear processes. Comput. Chem. Eng. 2007, 31 (5−6), 410−418. (20) Manenti, F. Considerations on nonlinear model predictive control techniques. Comput. Chem. Eng. 2011, 35 (11), 2491−2509. (21) Christofides, P. D.; Scattolini, R.; de la Pena, D. M.; Liu, J. F. Distributed model predictive control: A tutorial review and future research directions. Comput. Chem. Eng. 2013, 51, 21−41. (22) Manenti, F.; Cieri, S.; Restelli, M. Considerations on the steadystate modeling of methanol synthesis fixed-bed reactor. Chem. Eng. Sci. 2011, 66 (2), 152−162. (23) Manenti, F.; Cieri, S.; Restelli, M.; Bozzano, G. Dynamic Modelling of the Methanol Synthesis Fixed-Bed Reactor. Comput. Chem. Eng. 2013, 48, 325−334. (24) Qin, S. J.; Badgwell, T. A. An Overview of Nonlinear Model Predictive Control Applications; Allgower, F.; Zheng, A., Eds.; Birkhauser: Switzerland, 2000. (25) Qin, S. J.; Badgwell, T. A. A survey of industrial model predictive control technology. Control Eng. Pract. 2003, 11 (7), 733−764. (26) Findeisen, R.; Allgower, R. The quasi-infinite horizon approach to nonlinear model predictive control. Nonlinear Adaptive Control 2003, 281, 89−108. (27) Morari, M.; Lee, J. H. Model predictive control: past, present and future. Comput. Chem. Eng. 1999, 23 (4−5), 667−682. (28) Biegler, L. T. Large-scale Nonlinear Programming: an Integrating Framework for Enterprise-wide Dynamic Optimization. Proceedings of ESCAPE-17, Bucharest, Romania, 2007; Elsevier: Amsterdam, 2007; pp 575−582. (29) Varma, V. A.; Reklaitis, G. V.; Blau, G. E.; Pekny, J. F. Enterprisewide modeling & optimization - An overview of emerging research challenges and opportunities. Comput. Chem. Eng. 2007, 31 (5−6), 692−711. (30) Cucek, L.; Klemeš, J. J.; Kravanja, Z. A review of footprint analysis tools for monitoring impacts on sustainability. J. Cleaner Prod. 2012, 34, 9−20. (31) Klemeš, J. J.; Varbanov, P. S.; Pierucci, S.; Huisingh, D. Minimising emissions and energy wastage by improved industrial processes and integration of renewable energy. J. Cleaner Prod. 2010, 18 (9), 843−847. (32) Lam, H. L.; Varbanov, P. S.; Klemeš, J. J. Optimisation of regional energy supply chains utilising renewables: P-graph approach. Comput. Chem. Eng. 2010, 34 (5), 782−792.

(33) Vaccari, G.; Tamburini, E.; Sgualdino, G.; Urbaniec, K.; Klemeš, J. Overview of the environmental problems in beet sugar processing: Possible solutions. J. Cleaner Prod. 2005, 13 (5), 499−507. (34) Manenti, F.; Rovaglio, M. Market-driven Operational Optimization of Industrial Gas Supply Chains. Comput. Chem. Eng. 2013, 56, 128−141. (35) Manenti, F.; Bozzano, G.; D'Isanto, M.; Nascimento Lima, N. M.; Linan, L. Z. Raising the decision-making level to improve the enterprisewide production flexibility. AIChE J. 2013, 59 (5), 1588−1598. (36) Mayne, D. Q.; Rawlings, J. B.; Rao, C. V.; Scokaert, P. O. M. Constrained model predictive control: Stability and optimality. Automatica 2000, 36 (6), 789−814. (37) Seborg, D. E.; Edgar, T. F.; Mellichamp, D. A. Process Dynamics and Control, 2nd ed.; John Wiley and Sons: New York, 2003. (38) Buzzi-Ferraris, G.; Manenti, F. A Combination of Parallel Computing and Object-Oriented Programming to Improve Optimizer Robustness and Efficiency. Comput.-Aided Chem. Eng. 2010, 28, 337− 342. (39) Manenti, F.; Grottoli, M. G.; Pierucci, S. Online Data Reconciliation with Poor Redundancy Systems. Ind. Eng. Chem. Res. 2011, 50 (24), 14105−14114. (40) Manenti, F.; Signor, S.; Grottoli, M. G.; Fabbri, P. Adaptive Data Reconciliation Coupling C++ and PRO/II and On-line Application by the Field. Comput.-Aided Chem. Eng. 2010, 28, 373−378. (41) Signor, S.; Manenti, F.; Grottoli, M. G.; Fabbri, P.; Pierucci, S. Sulfur Recovery Units: Adaptive Simulation and Model Validation on an Industrial Plant. Ind. Eng. Chem. Res. 2010, 49 (12), 5714−5724. (42) Manenti, F.; Buzzi-Ferraris, G. Criteria for Outliers Detection in Nonlinear Regression Problems. Comput.-Aided Chem. Eng. 2009, 26, 913−917. (43) Buzzi-Ferraris, G.; Manenti, F. Outlier Detection in Large Data Sets. Comput. Chem. Eng. 2010, 35, 388−390. (44) Busch, J.; Oldenburg, J.; Santos, M.; Cruse, A.; Marquardt, W. Dynamic Predictive Scheduling of Operational Strategies for Continuous Processes Using Mixed-logic Dynamic Optimization. Comput. Chem. Eng. 2007, 31, 574−587. (45) Binder, T.; Blank, L.; Bock, H. G.; Bulirsch, R.; Dahmen, W.; Diehl, M.; Kronseder, T.; Marquardt, W.; Schloder, J. P.; Van Stryk, O., Introduction to Model Based Optimization of Chemical Processes on Moving Horizons. Online Optimization of Large Scale Systems: State of the Art; Groetschel, M. Krumke, S. O., Rambau, J. Eds.; Springer: New York, 2001. (46) Henson, M. A.; Seborg, D. E. Nonlinear Process Control; Prentice Hall: Upper Saddle River, NJ, 1997. (47) Imsland, L.; Findeisen, R.; Bullinger, E.; Allgower, F.; Foss, B. A. A note on stability, robustness and performance of output feedback nonlinear model predictive control. J. Process Control 2003, 13 (7), 633− 644. (48) Biegler, L. T.; Hughes, R. R. Feasible Path Optimization with Sequential Modular Simulators. Comput. Chem. Eng. 1985, 9, 379−394. (49) Biegler, L. T.; Hughes, R. R. Infeasible Path Optimization with Sequential Modular Simulators. AIChE J. 1982, 28 (6), 994−1002. (50) Manenti, F.; Dones, I.; Buzzi-Ferraris, G.; Preisig, H. A. Efficient Numerical Solver of Partially Structured DAE Systems. Ind. Eng. Chem. Res. 2009, 48, 9979−9984. (51) Bendersky, E.; Christofides, P. D. Optimization of transportreaction processes using nonlinear model reduction. Chem. Eng. Sci. 2000, 55 (19), 4349−4366. (52) Christofides, P. D.; Daoutidis, P. Nonlinear control of diffusionconvection-reaction processes. Comput. Chem. Eng. 1996, 20, S1071− S1076. (53) Liu, J.; de la Pena, D. M.; Christofides, P. D. Distributed Model Predictive Control of Nonlinear Process Systems. AIChE J. 2009, 55 (5), 1171−1184. (54) Lou, Y. M.; Hu, G. S.; Christofides, P. D. Model predictive control of nonlinear stochastic partial differential equations with application to a sputtering process. AIChE J. 2008, 54 (8), 2065−2081. 13090

dx.doi.org/10.1021/ie401511e | Ind. Eng. Chem. Res. 2013, 52, 13079−13091

Industrial & Engineering Chemistry Research

Article

(55) Ni, D.; Christofides, P. D. Multivariable predictive control of thin film deposition using a stochastic PDE model. Ind. Eng. Chem. Res. 2005, 44 (8), 2416−2427. (56) Palavajjhala, S.; Motard, R. L.; Joseph, B. Blocking and Condensing design for quadratic dynamic matrix control using wavelets. Ind. Eng. Chem. Res. 1994, 33 (5), 1159−1173. (57) Cagienard, R.; Grieder, P.; Kerrigan, E. C.; Morari, M. Move blocking strategies in receding horizon control. J. Process Control 2007, 17 (6), 563−570. (58) Lange, J. P. Methanol synthesis: a short review of technology improvements. Catal. Today 2001, 64 (1−2), 3−8. (59) Olah, G. A.; Goeppert, A.; Surya Prakash, G. K. Beyond Oil and Gas: The Methanol Economy; Wiley-VCH: Weinheim, Germany, 2009. (60) Lurgi GmbH, Lurgi MegaMethanol. www.lurgi.com. (61) Hartig, F.; Keil, F. J. Large-Scale Spherical Fixed-Bed Reactors Modeling and Optimization. Ind. Eng. Chem. Res. 1993, 32 (3), 424− 437. (62) Rezaie, N.; Jahanmiri, A.; Moghtaderi, B.; Rahimpour, M. R. A comparison of homogeneous and heterogeneous dynamic models for industrial methanol reactors in the presence of catalyst deactivation. Chem. Eng. Process. 2005, 44 (8), 911−921. (63) Buzzi-Ferraris, G.; Manenti, F. BzzMath: Library Overview and Recent Advances in Numerical Methods. Comput.-Aided Chem. Eng. 2012, 30 (2), 1312−1316. (64) Buzzi-Ferraris, G.; Manenti, F. Improving the selection of interior points for one-dimensional finite element methods. Comput. Chem. Eng. 2012, 40, 41−44. (65) Graaf, G. H.; Sijtsema, P. J.; Stamhuis, E. J.; Joosten, G. E. Chemical Equilibria in Methanol Synthesis. Chem. Eng. Sci. 1986, 41, 2883−2890. (66) Graaf, G. H.; Stamhuis, E. J.; Beenackers, A. Kinetics of Lowpressure Methanol Synthesis. Chem. Eng. Sci. 1988, 43, 3185−3195. (67) Lommerts, B. J.; Graaf, G. H.; Beenackers, A. Mathematical modeling of internal mass transport limitations in methanol synthesis. Chem. Eng. Sci. 2000, 55 (23), 5589−5598. (68) Narasimhan, S.; Jordache, C. Data Reconciliation & Gross Error Detection. An Intelligent Use of Process Data; Gulf Publishing Co.: Houston, TX, 2000. (69) Romagnoli, J. A.; Sanchez, M. C. Data Processing & Reconciliation for Chemical Process Operations. Process Systems Engineering; Elsevier: Amsterdam, 1999; Vol. 2. (70) Logist, F.; Saucez, P.; Van Impe, J.; Vande Wouwer, A. Simulation of (bio)chemical processes with distributed parameters using Matlab. Chem. Eng. J. 2009, 155 (3), 603−616. (71) Ullmann’s Modeling and Simulation; Wiley-VCH: Weinheim, Germany, 2007. (72) Chung, T.-H.; Ajlan, M.; Lee, L. L.; Starling, K. E. Generalized Multiparameter Correlation for Nonpolar and Polar Fluid Transport Properties. Ind. Eng. Chem. Res. 1988, 27, 671. (73) Lucas, K. Phase Equilibria and Fluid Properties in the Chemical Industry; Dechema, Frankfurt, 1980; p 573. (74) Lucas, K. Chem. Ing. Tech. 1981, 53, 959. (75) Soave, G. Chem. Eng. Sci. 1972, 27, 1197. (76) Stiel, L. I.; Thodos, G. AIChE J. 1964, 10, 26. (77) Mostinsky, I. L. Teplenergetika 1963, 4, 66. (78) Gupta, A. S.; Thodos, G. Direct analogy between mass and heat transfer to beds of spheres. AIChE J. 1963, 9, 751. (79) Gupta, A. S.; Thodos, G. Direct analogy between mass and heat transfer to pressure gasifier. Ind. Eng. Chem. Fundam. 1964, 3, 218. (80) Fuller, E. N.; Schettler, P. D.; Giddings, J. C. New method for prediction of binary gas-phase diffusion coefficients. Ind. Eng. Chem. Res. 1966, 5, 18−27.

13091

dx.doi.org/10.1021/ie401511e | Ind. Eng. Chem. Res. 2013, 52, 13079−13091