Simulation and optimization of cryogenic air separation units using a homotopy-based backtracking method

Separation and Purification Technology 67 (2009) 262–270

Contents lists available at ScienceDirect

Separation and Purification Technology journal homepage: www.elsevier.com/locate/seppur

Simulation and optimization of cryogenic air separation units using a homotopy-based backtracking method Lingyu Zhu a,b , Zhiqiang Chen a , Xi Chen a,∗ , Zhijiang Shao a , Jixin Qian a a b

State Key Lab of Industrial Control Technology, Department of Control Science and Engineering, Zhejiang University, Hangzhou 310027, China College of Chemical Engineering and Materials Science, Zhejiang University of Technology, Hangzhou 310014, China

a r t i c l e

i n f o

Article history: Received 28 November 2008 Received in revised form 10 February 2009 Accepted 11 March 2009 Keywords: Cryogenic air separation Homotopy method Backtracking Automatic load change

a b s t r a c t In response to frequently changing demands on gaseous products of cryogenic air separation units, an automatic load change (ALC) system is desired to optimally change the products and recycle rates under widely variable load conditions. However, due to the complex heat integration in the process, simulation and optimization of cryogenic air separation units (ASUs) often fails to converge with traditional Newtonbased algorithms. A homotopy-based backtracking method (HBM) is presented and applied to the process operation of a cryogenic ASU under wide changes in load conditions. Different from traditional process simulation, a load vector is introduced to denote the variation of the operation point due to the load change. The HBM solves the process in a way that heads toward the target and backtracks at failure. The results show that a large number of operating points that failed to converge with traditional algorithms can be successfully optimized with the proposed method. © 2009 Elsevier B.V. All rights reserved.

1. Introduction A cryogenic air separation unit (ASU) produces large volumes of oxygen, nitrogen, and argon at high purity. It is always connected to a manufacturing process such as production of primary metals or chemicals, or gasification. If the ASU is supplying gaseous products to a batch process like a steel furnace, it must automatically and rapidly respond to changing demands from the terminal user. Accurate predictions of the steady states at various operating points during load changes are essential for process operation. This requires a simulation and optimization technique that can always converge quickly and successfully. A number of research reports on the process simulation and control of cryogenic air separation units have been published. Vinson [1] discussed air separation control technologies and considered the installation of model predictive control (MPC) as the best industry practice. Zhu et al. [2] developed a reduced order model for nitrogen purification using nonlinear wave theory. As an extension of this study, Bian et al. [3] established a dynamic simulation of the column and the integrated condenser/reboiler using Aspen Dynamics® and applied it to study the dynamic characteristics of the nitrogen plant. Bian et al. [4] also proposed compartmental modeling of high purity air separation columns by simplifying the stage-by-stage balance equations to compartment-by-compartment balance equations. These studies normally focused on columns. However, the sepa-

∗ Corresponding author. Tel.: +86 571 87953068; fax: +86 571 87951068. E-mail address: [email protected] (X. Chen). 1383-5866/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.seppur.2009.03.032

ration of air into its components is energy-intensive and most units used in this process are highly thermally coupled. Therefore, the whole plant needs to be considered to reveal its true characteristics. Sirdeshpande et al. [5] described the mass and energy balances of the air separation cycle with a simplified algebraic model. However, the rigorous equilibrium-based method for a multi-components distillation column is much more effective for providing basic data to design and optimize this application. The traditional sequential modular (SM) method with iterating and tearing techniques often fails to converge quickly for this complex process with sensitive heat-coupling flowsheeting features under such widely varied operating conditions. The equation-oriented (EO) method based on Newton’s algorithm, on the other hand, typically converges quickly. However, the EO method requires a good initial guess close to the solution, especially for complex process simulations. Therefore, Newton’s method has difficulties in cases with wide and frequent load change requirements. Homotopy methods [6] are a class of robust and globally convergent algorithms. The idea of these algorithms is to continuously deform a simple (easy) problem into the given (hard) problem with a homotopy index changing from zero to one. If the family of deformed problems converges easily, the method is particularly suitable for highly nonlinear problem for which good initial solution estimates are difficult to obtain. This kind of method has been successfully applied to chemical engineering process simulations [7–15]. Unlike the traditional homotopy methods for static process simulations, we propose in this paper a homotopy-based backtracking method (HBM) for real time optimization of a cryogenic ASU under widely varying load conditions. This method uses the opera-

L. Zhu et al. / Separation and Purification Technology 67 (2009) 262–270

263

Fig. 1. Flow diagram of an internally compressed cryogenic air separation unit.

tional difference between the base and target points during a load change to define the homotopy parameter. Moreover, a backtracking technique is applied during the transition from the base point to the target point to improve efficiency. 2. Problem description In this project, the flowsheeting study focuses on an internally compressed cryogenic air separation plant with a nominal capacity of 20,000 Nm3 /h gaseous oxygen. As illustrated in Fig. 1, the flowsheet consists of cooling and separating sections in a cold box. The air compression and CO2 removal sections are not under discussion in this project. In the plant, the compressed and cooled air streams are distilled in an integrated four-column distillation system, which consists of a high-pressure column (HPC), a low-pressure column (LPC), an argon sidearm column (ASC), and an argon distillation column (ADC). The cleaned air is split into three fractions: high-pressure air (HPA), main air (MA), and turbine air (TA). All the feeds are sent to a multistream heat exchanger, E1, to exchange heat with exit products and the waste nitrogen. HPA is mostly liquefied through a throttle valve and sent to an intermediate location of the HPC. MA and TA are directly fed into the bottom of HPC as saturated vapors. In HPC, the air is separated into high purity liquid nitrogen, nitrogen-rich liquid, and oxygen-rich liquid streams. LPC produces high purity gaseous nitrogen (GAN) at the top and liquid oxygen (LOX) at the bottom. After part of the LPC bottom stream is withdrawn as the liquid oxygen product, the remaining portion is pumped with an elevated pressure and vaporized in E1 to produce gaseous oxygen (GOX). Waste nitrogen (WN) is a side product of LPC drawn below the first section of the packing. A sub-cooler, E2, is designed to cool the HPC product streams against WN and GAN. The crude argon (AR) in this flowsheet is produced at the top of the ASC. High purity gaseous argon (GAR) and liquid argon (LAR) products are distilled in the ADC equipment. The process is a highly energy-integrated system. All columns and heat exchangers are thermally coupled without any utility inputs. HPC and LPC are designed to share a common condenser/reboiler. The overhead of the ASC is condensed by the oxygen-rich liquid air. ADC is reboiled by nitrogen vapor, called

the pressure nitrogen (PRGAN) and condensed by liquid nitrogen. Gaseous products are evaporated and warmed by recycles and feeds in the heat exchangers E2 and E1. The ASU is expected to optimally change its operating point in response to demand changes with an advanced control system. To maximize profit, a process simulation and optimization system is developed in this project for setting variables during load changes. In this system, the liquefied nitrogen and both the gaseous and liquefied argon products are set to constant flow rates. The product mix, therefore, is determined by three other variables, {FGOX , FLOX , FGAN }, which represent the flow rates of gaseous oxygen, liquefied oxygen and gaseous nitrogen, respectively. Each product mix represents the online gas and liquid demand of a customer. Besides the demands of the product mix, the economic objective of minimizing the compression cost on feeds is also desirable for the process operation. With a steady state process model (f) described by the balance equations of the material, equilibrium, summation, and enthalpy (MESH) for all units and an additional work equation for the expander, the following process optimization problem can be built to locate the optimal operating point at various loads. mink1 FHPA + k2 FMA + k3 FTA Xopt

s.t. : f (Xprd , Xopt , Xfix , Xexp , Xcal ) = 0 Xprd = {FGOX , FLOX , FGAN } U L ≤X Xcal cal ≤ Xcal U L ≤X Xopt opt ≤ Xopt

(1)

where k1 , k2 , and k3 are cost coefficients of compressing HPA, MA and TA, respectively; Xprd is the desired product mix at a changing load; Xopt represents the optimized variables; Xfix represents the fixed parameters; Xexp represents the empirical variables; Xcal L and X U are the bounds on represents the calculated variables; Xcal cal the calculated variables, such as the product purity specifications, L and X U operating pressure and temperature limits, etc.; and Xopt opt are the bounds on the optimized variables. The mathematical model of the process is built in Aspen Plus by a set of nonlinear algebraic equations with a total of 7447 variables. Among these variables, the fixed parameters include equipment structure parameters and hypothetically constant variables such

264

L. Zhu et al. / Separation and Purification Technology 67 (2009) 262–270

Table 1 Description of constant variables.

simulation. In this project, we define the homotopy parameter as

Variable

Description

Dimension

Xprd Xopt Xfix

{FGOX , FLOX , FGAN } {FTA , FMA , FHPA , FPRGAN , FRELIN , FLWN , FAR , FVENT , FAR RECYCLE } Feed composition Equipment structure parameters Temperature Pressure and pressure drop Efficiency Heat loss

3 9 9 227 12 40 10 15

Xexp

h(x, t) = f (x, ˛(t)) = 0

3. Homotopy-based backtracking method When the demand on an ASU is changed, the process is required to operate at a new point obtained by the process optimization. The traditional Newton-based methods, however, have difficulties in this task, because their convergence relies heavily on a good initial guess. The high heat-coupling feature of the process results in a very small feasible region. The traditional Newton-based methods always fail to converge without a good initial guess. The homotopy (or continuation) method [6–8] is a type of method for process simulation that overcomes the dependence on a good initial guess. The general concept of the homotopy method is to transform the original problem of solving algebraic equations, f(x) = 0, into solving a new equation set, h(x, t) = 0. The parameter t is the homotopy parameter. The homotopy function, h, is defined as (2)

It is designed so the problem can be easily solved when t equals 0, i.e., for h(x, 0) = 0. On the other hand, when t equals 1, h(x, 1) = 0 is equivalent to f(x) = 0. Therefore, the solution of f(x) = 0 is equivalent to the solution of h(x, 1) = 0. When t varies from 0 to 1, a solution path, x*(t) is created for h(x, t) = 0. The strategy of the homotopy method is to track the solution path from x*(0) and gradually reach solution x*(1), which is also the solution to the original problem. In this project, we propose a method for process simulation and optimization with various load changes and apply it to the operation of the cryogenic ASU described in the previous section. Different from traditional process simulation, an extra load vector, ˛, is introduced to denote the variation of the operation point due to the load change. The process simulation at operation point ˛ is indicated as follows: f (x, ˛) = 0

˛ − ˛bp ˛tp − ˛bp

(4)

which represents the ratio of the distance between an intermediate operation point and the target operation point over the distance between the base operation point and the target operation point. The homotopy function is defined as

as the heat duty of the throttle valves. These parameters maintain the same values as those at the base operating point during the load change. The empirical variables are equipment performance parameters affected by the load. Pressure drop rate and efficiency rate are two good examples. These parameters are changed by empirical correlations. The optimized variable set includes nine controllable flow rates of feeds, recycles and vents, which need to be adjusted in conjunction with the load. Information on each type of variables is summarized in Table 1.

h(x, t) = tf (x) + (1 − t)g(x)

t=

(3)

where the solution to the above problem, x*, is a function of the load vector, ˛. Two load vectors, ˛bp and ˛tp , are used to denote the base point and the target point, respectively. Given the solution at the base point, x*(˛bp ), we can develop a homotopy method to reach the solution at the target point, x*(˛tp ). There are two key issues in this method. First, how should the homotopy function and parameters be defined? Second, what strategy should be used to efficiently control the search to track the solution path? The construction of the homotopy function and the selection of the homotopy parameter are always big challenges in process

(5)

When t equals 0, the solution is known as x*(˛bp ). When t equals 1, the solution is identical to the expected solution at the target operating point, x*(˛tp ). Unlike the traditional homotopy method that changes the homotopy parameter gradually from 0 to 1, a backtracking strategy is proposed in this project to deal with the second key issue mentioned above. The proposed homotopy-based backtracking method (HBM) solves the process in a way that heads toward the target and backtracks at failure. When the process is adjusted from a base operating point, ˛bp , to a new target point, ˛tp , h(x, 1) = f (x, ˛(1)) = f (x, ˛tp ) = 0 is first solved with the solution x*(˛bp ) as an initial guess, where t equals 1 in this step. If the load change is significant, the process may not converge with the initial solution. Then the simulation is conducted at an intermediate operating point by backtracking the homotopy parameter, t, with ratio d (d > 1). This parameter d is called the backtracking constant. It is set to 2 in this project unless otherwise specified. The new operating point represented by the load parameters is as follows: t=

t d

˛(t) = ˛bp + t ∗ (˛tp − ˛bp )

(6) (7)

If the process successfully converges at this new point, this point will be set as the base point; the simulation will then test the target point again with the solution at the new location as the initial guess. Whenever the process fails to converge, backtracking will be conducted to a point closer to the base point. This process will continue until the target point is successfully solved, or it reaches a boundary which cannot be surpassed. Fig. 2 presents the flowchart of the proposed algorithm. It should be noted that the HBM also requires another algorithm such as the Newton-based method as a kernel. It can be proved that the HBM with Newton’s method as the core algorithm can converge to the target with any backtracking constant d > 1 under appropriate assumptions. Due to the limited scope of this paper, the proof is omitted. In this project, HBM is implemented via the Open Object Model Framework (OOMF) [16], which is a kernel component of Aspen Plus. The structure of the implementation is illustrated in Fig. 3. HBM is built in an OOMF script in Aspen Plus to solve NLA or NLP problems with algorithm LSSQP. OpenSolver DLLs based on the AOS interface are used to link to external data exchange. It should be noted that this method is not limited to process simulation; it can also be extended to process optimization as shown later in this paper. 4. Results and discussion The industrial ASU illustrated in Fig. 1 consists of 33 pieces of equipment, covering 6 types including distillation column, heat exchanger, valve, mixer, splitter and expander. A mathematical model was constructed in Aspen Plus from rigorous material and energy balances for the plant. Stream results at the base operating point are listed in Table 2. The product mix Xprd at the base operating point is {20,000, 1600, 40,193.32}. The molar fractions in bold type list those key product compositions with specifications. The upper limit of the nitrogen residual in the oxygen product is 2.5 ppm.

L. Zhu et al. / Separation and Purification Technology 67 (2009) 262–270

265

Fig. 2. Flowchart of HBM.

in the system. A feasibility study around the designed base operating point helped to understand the special features of this system. 4.1. Feasibility study around base operating point

Fig. 3. Implementation structure of HBM in AspenPlus® .

The purity of the oxygen product is specified to be higher than 99.8%. Oxygen and nitrogen residuals in the argon product need to be limited to 1.5 ppm and 3.5 ppm, respectively. These specifications are all satisfied in the simulation result at the base operating point. When we conducted process optimization with the algorithms in Aspen Plus, it failed to converge. This failure was mainly due to the narrow feasible region caused by the high heat-coupling

As shown in Eq. (1), the feasible region is constrained by the MESH model and the product specifications. Since it is difficult to visualize the feasible region in a high-dimensional space, we first explored it by studying the effect of single variable changes. Keeping the other variables at their fixed nominal values, we changed one variable of either Xprd or Xopt at a time to see if the product quality still met the specifications. The test for each variable was conducted by gradually increasing the change ratio until a violation of the constraints was observed, or the variable had changed by ±10%. The feasible regions of each variable are described by the shadow columns shown in Fig. 4. The results are divided into two parts; the upper part of the figure includes the variables with a relatively large fluctuation ranges; the lower part includes the other variables with very limited tolerance of change. Among the variables with large fluctuation range, FPRGAN denotes the flow rate of pressure nitrogen, which is a small recycle heat exchange stream with very few influences on the flowsheet. The rates of gaseous nitrogen product (FGAN ) could be decreased down to −10% due to the discharge of nitrogen in waste nitrogen. The rate of the ADC vent could also be reduced down to −10% due to the increase of argon discharge. The results show that most variables, including the flow rates of the air feeds and

Table 2 Stream results at base operating point. Stream

HPA

MA

TA

GAN

LIN

GOX

LOX

GAR

LAR

WN

VENT

Phase

Mixed

Vapor

Vapor

Vapor

Liquid

Vapor

Liquid

Vapor

Liquid

Vapor

Vapor

Molar fraction O2 Ar N2 P, MPa T, K F, Nm3 /h

0.2095 0.0093 0.7812 0.56 97.70 28,600

0.2095 0.0093 0.7812 0.57 103.40 45,980

0.2095 0.0093 0.7812 0.57 100.00 28,720

1.6 ppm 0.0001 0.9999 1.15 306.60 40,193.32

2 ppm 0.0002 0.9998 1.317 82 0

0.9990 0.0010 0 ppm 31.00 306.60 20,000

0.9990 0.0010 0 ppm 1.41 93.26 1600

0.6 ppm 0.999997 2.2 ppm 31.00 306.60 400

0.6 ppm 0.999997 2.2 ppm 1.50 91.13 383.7034

0.0008 0.0043 0.9949 1.36 80.14 40,772.5

4.39E−07 0.595800 4.04E−01 1.28 85.01 8

266

L. Zhu et al. / Separation and Purification Technology 67 (2009) 262–270

Fig. 4. Feasibility study of single variable.

oxygen product, could only fluctuate within a small range of ±1%, which indicates that the feasible region is very narrow along many directions. A feasibility study with multiple variables was then conducted on the product mix variables, {FGOX , FLOX , FGAN }. The range of load change of each product was within −5% to +5% of the base point. Fig. 5 shows the three-dimensional results. Each axis represents the change ratio of one variable. When the other variables were fixed at the base operating point, the product mix could change in a small range without violating the constraints. In the figure, all feasible points are denoted by dotted symbols, except the base point denoted as a star. Compared with the single variable results, this three-dimensional analysis discovered a larger feasible region as there were more degrees of freedom. However, this feasible region was still limited with the flow rate of GOX in a range of −0.06% to +0.07% and LOX in a range of −0.8% to +0.6%. Another three-dimensional feasibility analysis was conducted on the flow rates of the three air feeds, FHPA , FMA , and FTA . By maintaining the values of other variables at the base operating point, the three feeds could slightly change their flow rates while still fitting the model and meeting the product specifications as shown in Fig. 6. However, the largest feasible change along each direction was no more than −0.4% to 0.4%, indicating a narrow feasible region.

Fig. 5. Feasibility study of Xprd .

4.2. Optimization under the base operating condition The simulation result at the base operating point is not the optimal one for the following reasons. First of all, the operating point was obtained through calculation in the design mode rather than in the simulation mode. Second, the objective considered in the real process operation is different from that in the design stage. Therefore, we first tested the optimization problem in Eq. (1) at the base operating point. When we applied it in Aspen Plus using DMO or LSSQP algorithms, the optimization always failed to converge. Starting from the initial point based on the simulation result, optimization with the defined objective quickly led the search out of the feasible region. The algorithm failed to re-enter the narrow feasible region; thus it failed to reach the optimal solution. We know that the simulation result presented in the previous section is a feasible solution of the optimization. Therefore, the failure was not due to the feasibility; it was because these traditional gradient-based algorithms are incapable of efficiently dealing with this problem with a narrow feasible region. Therefore, we propose using HBM to implement the optimization at the base operating point. Unlike the original process simulation, the optimized variables of the optimization process were allowed to change within a bound represented by the constraints. This method used the inequality bounds of the three feed flow rates included in the objective function to

Fig. 6. Feasibility study of Xopt .

L. Zhu et al. / Separation and Purification Technology 67 (2009) 262–270

267

Fig. 7. Solving process of HBM.

set up the homotopy parameter. When the optimization failed, we reduced the bound to build a new optimization problem with a solution closer to the original simulation problem. From the convergence properties of the Newton-based algorithm, we know that when the bound is narrow enough, the optimization problem can be solved. Therefore, following the strategy of HBM described in Section 3, the optimal solution can be reached eventually from the initial simulation point by gradually adjusting the bounds. In this problem, the homotopy parameter t was defined to relate the lower and upper bounds of the three optimized variables, FTA , FHPA , and FMA . Denote Xopt1 as the set of these three optimized variables, and Xopt2 for the other optimized variables. We further use U L X¯ opt1 and X¯ opt1 to stand for the specified lower and upper bounds of bp

the optimization problem. Denote Xopt1 as the value of these three variables of the simulation result at the base operating point. We can thus build the homotopy-based optimization as min

Xopt1 ,Xopt2

k1 FHPA + k2 FMA + k3 FTA

s.t. : f (Xprd , Xopt , Xfix , Xexp , Xcal ) = 0 Xprd = {F bp , F bp , F bp } U L ≤X X¯ cal cal ≤ X cal ¯X L ≤ Xopt2 ≤ X¯ U

opt2 L Xopt1 (t) U (t) Xopt1 L Xopt1 (t)

≤ = =

(9)

opt2 U (t) Xopt1 ≤ Xopt1 bp bp U Xopt1 + t(X¯ opt1 − Xopt1 ) bp bp L Xopt1 + t(X¯ opt1 − Xopt1 )

When t equals one, the above problem is the same as the anticipated optimization problem described in Eq. (1); when t equals zero, the three optimized variables involved in the objective function are all fixed by the same bound constraints, which means that the previous simulation result at the base operating point is the optimal solution. Then the HBM was used to implement the optimization by changing t from 0 to 1 with the backtracking strategy. Fig. 7 illustrates the converging procedure with the change of the homotopy parameter. A solid square denotes a convergent optimization and a shadowed square stands for an unsuccessful optimization. The height of the square represents the time cost of the optimization; its center represents the value of the homotopy parameter; the width of the square is meaningless. The figure shows that when the solution at the target bound was first conducted using the base point simulation as an initial guess, it failed to converge; then the backtracking strategy was activated and the ratio of the homotopy parameters was halved. The optimization did not converge until the homotopy parameter had backtracked to 1/28 in the 9th test run. Then it headed toward the target again using the new solution as an initial guess. The whole process converged at the target in the 10th test run. The total time cost of this whole

Fig. 8. Homotopy path of variable bounds.

optimization procedure was about 210 s. The homotopy paths of U L and Xopt1 are illustrated in Fig. 8. The optimization problem Xopt1 converged only when the upper and lower bounds of the variable were very close to each other; then the optimization with the new initial guess successfully and quickly converged to the optimal point. 4.3. Varying load optimizations After conducting the optimization at the base operating point, we further studied process optimizations under varying loads. In real operations, the anticipated load range of each product is between 70% and 110% of the base point. If we discretize the load change range at 5% for each product, a total of 9 × 9 × 9, i.e., 729 operating points must be tested. All the points together constitute a cube. Experience tells us that not all combinations are feasible as the compositions of oxygen, nitrogen and argon in the feed air are fixed. The feasible region is a part of the cube. A robust algorithm is expected to quickly and successfully solve every feasible operating point during the load change. Optimization at each point of the above-mentioned cube was solved with initialization at the optimal solution of the base point. It was first conducted by using Aspen Plus directly. After a comparison of the algorithms Aspen Plus provides, LSSQP was selected for the optimization problems because of its better performance. As illustrated in Fig. 9, the dot symbols represent the convergent oper-

Fig. 9. Convergence results of operating points in the cube.

Converge status

LSSQP/HBM LSSQP/HBM HBM HBM HBM HBM Failed Failed Failed 0.356 0.355 0.353 0.351 0.350 0.348

L. Zhu et al. / Separation and Purification Technology 67 (2009) 262–270

FWN /FAIR

268

ating points that can be successfully solved by Aspen Plus, which cover 404 of the total 729 points in the cube. We then applied the proposed HBM to the same problem. The homotopy parameter was defined similarly to Eq. (4), which denotes the ratio of the distance between an intermediate operation point and the target operation point to the distance between the base operation point and the target operation point. Then the optimization with HBM under load change was formulated as mink1 FHPA + k2 FMA + k3 FTA

0.491 0.489 0.487 0.485 0.484 0.482 0.480 0.478 0.476

Product ratio

Xopt

s.t. : f (Xprd , Xopt , Xfix , Xexp , Xcal ) = 0 tp Xprd = {FGOX , FLOX , FGAN }

bp

tp

bp

Xprd = Xprd + t(Xprd − Xprd )

0.06 0.06 0.07 0.07 0.07 0.07

yN2,AR (ppm) yO2,AR (ppm)

1.50 1.50 1.50 1.50 1.50 1.50

0.02 0.04 0.21 3.50 3.50 3.50

tp

yO2,WN (%)

(10)

U L ≤X Xcal cal ≤ Xcal U L Xopt ≤ Xopt ≤ Xopt

bp

where Xprd is the specified target product flow rate, Xprd is the product flow rate at the base point, and t is the homotopy parameter. When t varies from 0 to 1, the solution also transits from the optimal base operating point to the optimal point under a new load. By applying HBM to various load optimizations, there were 127 additional convergent operating points in addition to those 404 convergent points already identified. As illustrated in Fig. 9, we can clearly see that HBM increased the convergent region. The lower success rate of the traditional method was mainly caused by the large differences between the new operating points and the base operating point, which were beyond the convergence region of Newton’s method. The proposed HBM, on the other hand, overcame the dependency on a good initial guess because it had a gradual approximation feature with backtracking technique in the case of failure.

yO2,N2 (ppm) yO2 (ppm)

0.998 0.998 0.998 0.998 0.998 0.998 21,738.3 21,364.1 20,989.8 20,613.5 24,091.4 24,492.8 31,591.7 31,590.1 31,591.3 31,591 27,735.4 26,946.3 50,578 50,578 50,578 50,578 50,578 50,577.9

FHPA (Nm3 /h) FTA (Nm3 /h) FMA (Nm3 /h)

As discussed in Section 3, the HBM can guarantee convergence unless there is an unsurpassable physical limit in the problem. In other words, the HBM is also useful to determine the feasibility of a proposed optimization problem during load change operations. Because of the fixed composition in the air feed and specifications on purity of each product, we cannot require an arbitrary product mix Xprd ; there is a reasonable range of each product flow rate. The HBM helps to locate the feasible boundary. The following analysis is presented to explain the reasonable relation between the oxygen and nitrogen product rates. The mass balance for nitrogen in the ASU can be established as follows: FAIR yN2,AIR = (FLOX + FGOX )yN2,O2 + (FLIN + FGAN )yN2 + FWN yN2,WN + (FLAR + FGAR )yN2,AR

(11)

where yN2,AIR is the fraction of nitrogen in the air, yN2,O2 is the residual fraction of nitrogen in the oxygen products, yN2 is the purity of the nitrogen products, yN2,WN is the residual fraction of nitrogen in the waste nitrogen, and yN2,AR is the fraction of nitrogen in the argon products. Since the oxygen, argon and nitrogen products are high purity products, we have yN2,O2 ≈ 0

(12)

yN2,AR ≈ 0

(13)

yN2 ≈ 1

(14)

The waste nitrogen has a purity specification: yN2,WN ≥ 0.995

(15)

Substituting Eqs. (12)–(15) into Eq. (11), we have 1760 1680 1600 1520 1440 1360 1280 1200 1120

FLOX (Nm3 /h)

Table 3 Optimization results for FGOX = 20,000 Nm3 /h, FGAN = 44,326 Nm3 /h.

2.50 2.50 2.50 2.50 2.50 2.50

4.4. Location of physical boundary

(FLIN + FGAN ) ≤ FAIR yN2,AIR − 0.995FWN

(16)

L. Zhu et al. / Separation and Purification Technology 67 (2009) 262–270

269

Similarly, the oxygen mass balance of the ASU can be established as: FAIR yO2,AIR = (FLOX + FGOX )yO2 + (FLIN + FGAN )yO2,N2 + FWN yO2,WN + (FLAR + FGAR )yO2,AR

(17)

where yO2,AIR is the fraction of oxygen in the air, yO2 is the purity of the oxygen products, yO2,N2 is the residual fraction of oxygen in the nitrogen products, yO2,WN is the residual fraction of oxygen in the waste nitrogen, and yO2,AR is the residual fraction of oxygen in the argon products. Again due to the high purity specifications in the oxygen, argon and nitrogen products, we have yO2,N2 ≈ 0

(18)

yO2,AR ≈ 0

(19)

yO2 ≈ 1

(20)

The oxygen residue specification in the waste nitrogen is: yO2,WN ≤ 0.2%

(21)

Fig. 11. Results of product ratio without HBM.

Substituting Eqs. (18)–(21) into Eq. (17) results in (FLOX + FGOX ) ≥ FAIR yO2,AIR − 0.2%FWN

(22)

Therefore, the ratio of oxygen products to nitrogen products satisfies: FAIR yO2,AIR − 0.2%FWN yO2,AIR − 0.2%FWN /FAIR FLOX + FGOX ≥ = FLIN + FGAN FAIR yN2,AIR − 0.995FWN yN2,AIR − 0.995FWN /FAIR (23) Eq. (23) means that there is a lower bound for the product ratio. We calculated the ratios for all points in the “cube”, including the points that converged and failed using HBM. The result of the convergence status is shown in Fig. 10. The horizontal axis represents the number of points. The vertical axis is the value of the ratio at each point. We can clearly identify a separating line between two types of points. The upper half is the result of the convergent points, and the lower half is the result of the failed points. This indicates that HBM reached the limit of (FLOX + FGOX )/(FLIN + FGAN ) at about 0.48. As a contrast, Fig. 11 presents the convergence status for all product mixes without HBM; there is no separating line between the convergent and failed points. We can further compare the performance of the algorithms with and without HBM on approaching the physical boundary as follows. As shown in the three-dimensional plot in Fig. 9, if we fix

Fig. 10. Results of product ratio with HBM.

FGOX and FGAN and vary FLOX , we obtain a group of data points vertical to the horizontal surface. For the following example case, where FGOX = 20,000 Nm3 /h, FGAN = 44,326 Nm3 /h, and FLOX decreased from 1760 Nm3 /h to 1120 Nm3 /h, the optimization results are listed in Table 3. The 12th column shows the convergence results. We can see that the LSSQP algorithm alone converged for the first and second operating points and failed for the others. The HBM, in contrast, could converge for the first six operating conditions. The bold type numbers in the table indicate that the specification bounds were reached. When the given flow rate of liquefied oxygen decreased, the residual nitrogen in the argon product, yN2,AR , increased. In the first two conditions, this variable changed slightly and the optimization problem could converge using LSSQP. When the conditions further changed, the calculated nitrogen residue in the argon product dramatically increased to its upper bound, and the traditional algorithm failed to converge. The HBM, however, could still converge until it reached an unsurpassable boundary. The physical boundary in terms of the product ratio (FLOX + FGOX )/(FLIN + FGAN ), can be observed in the 10th column of the table, which agrees with the result in Fig. 10. 5. Conclusion Frequent changes in demand for the oxygen product require a robust simulation program for the air separation unit. Traditional process simulation algorithms based on Newton’s method, however, can only converge in a small feasible region due to their heavy reliance on a good initial guess. A homotopy-based backtracking method is proposed in this paper for the simulation and optimization of air separation units. The method ensures convergence even when an initial guess is far away from the target operation point. The difference between the base and target operating points in the load change is used to define the homotopy parameter. Additionally, a backtracking technique is applied during the transition from the base point to the target point to improve efficiency. The method was successfully applied to process operations of an ASU under a wide range of load conditions. It helped the process optimization under the base operating point by gradually transforming the problem from the designed process simulation. Moreover, the varying load optimization result shows that a large number of operating points that failed to converge with traditional algorithms could be successfully optimized with the HBM. Further analysis showed that this method also helped to locate the boundary of physical feasibility for product mix.

270

L. Zhu et al. / Separation and Purification Technology 67 (2009) 262–270

Acknowledgements We gratefully acknowledge the financial support of National Basic Research Program of China (No. 2009CB320603), National Natural Science Foundation of China (No. 60704029), and 863 Program of China (No. 2007AA04Z192). References [1] D.R. Vinson, Comput. Chem. Eng. 30 (2006) 1436–1466. [2] G.Y. Zhu, M.A. Henson, L. Megan, Sep. Purif. Technol. 24 (2001) 467–487. [3] S. Bian, M.A. Henson, P. Belanger, L. Megan, Ind. Eng. Chem. Res. 44 (2005) 153–167.

[4] S. Bian, et al., Comput. Chem. Eng. 29 (2005) 2096–2109. [5] A.R. Sirdeshpande, et al., AIChE J. 51 (2005) 1190–1200. [6] S. Eilenberg, N.E. Steenrod, Foundations of Algebraic Topology, Princeton Univ. Press, 1952. [7] T. Wayburn, J.D. Seader, Comput. Chem. Eng. 11 (1987) 7–25. [8] J.D. Seader, et al., Comput. Chem. Eng. 14 (1990) 71–85. [9] K.S. Gritton, J.D. Seader, W.J. Lin, Comput. Chem. Eng. 25 (2001) 1003–1019. [10] F. Jalali, J.D. Seader, S. Khaleghi, Comput. Chem. Eng. 32 (2008) 2333–2345. [11] A.R. Ciric, P. Miao, Ind. Eng. Chem. Res. 33 (1994) 2738–2748. [12] K.T. Hana, K.J. Lee, E.S. Yoon, Comput. Chem. Eng. 17 (Suppl. 1) (1993) 479–484. [13] B.P. Singh, et al., Chem. Eng. Res. Des. 83 (2005) 959–968. [14] J.R. Paloschi, Comput. Chem. Eng. 19 (1995) 1243–1254. [15] J.R. Paloschi, Comput. Chem. Eng. 21 (1997) 531–541. [16] Aspen Technology Inc., OOMF 2004.1 OOMF Script Language Reference Manual, 2005.