Microsoft Word - 476hernandez.docx CHEMICAL ENGINEERING TRANSACTIONS VOL. 43, 2015 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Chief Editors: Sauro Pierucci, Jiří J. Klemeš Copyright © 2015, AIDIC Servizi S.r.l., ISBN 978-88-95608-34-1; ISSN 2283-9216 Global Solution Approaches for Biomass to Commodity Chemicals (BTCC) Investment Planning Problem Ismail Fahmi, Selen Cremaschi* Department of Chemical Engineering, The University of Tulsa, 800 South Tucker Drive, Tulsa, Oklahoma 74114, USA selen-cremaschi@utulsa.edu Incorporating biomass, a potential feedstock for chemical process industry, to our existing system will require significant investments on biomass conversion technologies. A recent model, a mixed integer nonlinear program (MINLP), addresses this investment planning problem. However, it has been observed that it becomes intractable quickly with the increase in the number of technologies. This paper investigates specialized solution approaches by exploiting the MINLP structure. Replacing bilinear terms with linearly segmented tight relaxations and nonlinear terms with linearly segmented under and upper estimators, solving the resulting MIP, and initializing DICOPT at the MIP solution solves the original MINLP efficiently. 1. Introduction Biomass has great potential as chemical process industry (CPI) feedstock because it is abundant, locally available, and renewable (Peres et al., 2013). Incorporating this relatively new feedstock to our existing CPI will require significant amounts of investments: (1) for increasing the efficiency of the technologies that can be used to convert biomass to commodity chemicals, and hence make them cost competitive to support the CPI system, and (2) for expanding the production capacities of these technologies to meet the current and future market demands. As such there is great opportunity for investigating how these investments will impact the evolution of the biomass-to-commodity-chemicals (BTCC) system. Recently, a new representation and the corresponding model is introduced to study the investment planning problem for incorporating biomass as a CPI feedstock (Fahmi et al., 2014). The investment planning problem involves a series of yes/no decisions for investing in a technology at a particular time, and associated with these decisions is the amount of investment. The resulting model is a non-convex mixed integer nonlinear program (MINLP). Fahmi et al. (2014) considered a small case study, the evolution of ethylene and propylene production from biomass (corn grain and corn stover) and naphtha for a planning horizon of 50 years. They assumed that ethylene can be produced through three different technologies: biomass fermentation, biomass gasification, and naphtha cracking, and propylene can be produced through two different technologies: naphtha cracking, and metathesis reaction of ethylene and butenes. The resulting MINLP, which had 11,502 constraints and 7051 variables, was implemented in GAMS V23.6.2, and solved using BARON V9.0.6 with a Dual Intel E5405 2.0 GHz processor and 8 GB RAM memory, and it took about 44 CPU hours to solve to a relative gap of 5 %. Although the MINLP was shown to be able to capture the effects of raw material costs, the capacity expansion costs, the maximum annual budget, the demand increase, the inflation rate, and the extraction coefficient of non-renewable materials on the optimum production plan and total cost, it is computationally expensive to solve and becomes intractable quickly as the number of technologies increase. This paper develops and tests specialized solution algorithms that exploit the structure of the BTCC investment planning model. The source of computational complexity in the BTCC investment planning model are nonconvex nonlinear constraints and bilinear terms. Two different approaches are used to convexify the bilinear terms: linearly segmented tight relaxations and radix-based relaxations. The nonlinear terms are replaced with linearly segmented under and upper estimators. Once these terms are relaxed, the solution of the resulting MIP is used as the initialization to solve the original MINLP with two local solvers, DICOPT and DOI: 10.3303/CET1543222 Please cite this article as: Fahmi I., Cremaschi S., 2015, Global solution approaches for biomass to commodity chemicals (btcc) investment planning problem, Chemical Engineering Transactions, 43, 1327-1332 DOI: 10.3303/CET1543222 1327 SBB. We used ethylene production from naphtha and biomass as the case study and compared the performance of both approaches in terms of overall solution time and the quality of the solution. 2. The BTCC investment planning model Given a number of biomass and fossil-based feedstock processing technologies and their characteristics, and the initial commodity-chemicals-production-system market conditions, the MINLP determines the amount and schedule of investments (both R&D and production capacity expansions) for each technology to yield the minimum cost over the planning horizon. The initial commodity-chemicals-production-system market conditions are defined by initial raw material costs (CRv,0) of biomass and fossil based feedstocks, the current and future forecasted demands for the products, the current capacities (CXe,0) and total R&D expenditure (CRDe,0), and current costs of capacity expansions (CCe,0) of technologies that use fossil-based feedstocks and biomass. For each technology e, the relationship between its future cost and its capacity and R&D investments are defined via a two-factor learning curve (Kouvaritakis et al., 2000). These factors are learning- by-doing elasticity (αe) and learning-by-searching elasticity (βe), which describe how the expansion cost changes with changes in the total installed production capacity and R&D expenditure, respectively. The BTCC investment planning model is given in Figure 1. Details of the BTCC investment planning model can be found in Fahmi et al. (2014). 3. Structure of the BTCC investment planning model Most terms in the constraints of the BTCC investment planning model are linear (Figure 1). The nonlinear equations are the objective function, the two-factor learning curve that defines how the technology expansion costs change with changes in the capacity expansions and R&D expenditures (Figure 1, Technology Costs constraint), and the constraint that limits the production to technologies that has reached a certain maturity level (Figure 1, Capacity Constraints). The objective function has two nonlinear terms, one in calculation of the expansion cost and other in the calculation of the raw material costs. We can linearize the objective function by separating the nonlinear terms as shown in Eq(1). This transformation adds two equality constraints, Eq(1b) and Eq(1c), which contain bilinear terms to the MINLP.        eexpenditur D&Rcost material Rawcost Expansion  + − + + + + = − ∈ e t t tete VRv t t tv e t t te i CRDCRD i BLTerm i BLTerm TC )1( )( )1( 2 )1( 1 1000 1,,,, (1a) tetete XCXBLTerm ,,,1 = et,∀ (1b) tvtvte RCRBLTerm ,,,2 = VRvt ∈∀ , (1c) Following a similar logic, the Technology Cost constraint can be replaced with the constraint set given in Eq(2). This separation operation yields one linear equality constraint (Eq(2a)), and exposes three equality constraints, one with a bilinear term (Eq(2b)), and the other two with nonlinear terms (Eq(2c) and Eq(2d)). teete BLTermCXCC ,0,, 3= et,∀ (2a) tetete NLTermNLTermBLTerm ,,, 213 ⋅= et,∀ (2b) e e te te CX CX NLTerm α         = 0, , ,1 et,∀ (2c) e e te te CRD CRD NLTerm β         = 0, , ,2 et,∀ (2d) The Capacity Constraint (Figure 1) includes another bilinear term and is reformulated as shown in Eq(3): tete BLTermP ,, 4≤ et,∀ (3a) tetete CXYBLTerm ,,3,,4 ⋅= et,∀ (3b) With these modifications, it is observed that the BTCC investment planning problem has two groups of nonlinear terms (NLterm1e,t and NLterm2e,t) and four groups of bilinear terms (BLterm1e,t, BLterm2e,t, BLterm3e,t, and BLterm4e,t). Because the values of learning elasticities, αe and βe, are negative, the nonlinear terms monotonically decrease with zero asymptotic value and the resulting curve is concave up. Both 1328 nonlinear terms contribute to the non-convexity of the model as they are in equality constraints. The other computational complexity arises due to the bilinear terms of the problem. Figure 1: The MINLP of the BTCC investment planning problem Objective Function          eexpenditur D&R cost material Raw cost Expansion Min    + − + + + + − = − ∈ − e t t tete VRv t t tvtv e t t tetete i CRDCRD i RCR i CXCXCC TC )1( )( )1( )1( )( 1000 1,, ,, 1,,, (4a) Subject to Technology Costs ee e te e te ete CRD CRD CX CX CCCC βα                 = 0, , 0, , 0,, et,∀ (4b) Raw Material Costs  = += t j jvvvtv RkCRCR 1 ,0,, { }VRRvVRvvt ∉∧∈∀ , (4c) ( )tvtv IRCRCR += 10,, VRRvt ∈∀ , (4d) Product Demands t vvtv DD )1(0,, γ+= VPvt ∈∀ , (4e) Meet Product Demands tvtv DR ,, ≥ VPvt ∈∀ , (4f) = e teevtv PbR ,,, VPvt ∈∀ , (4g) No Accumulation of Intermediates 0,, = e teev Pb { }VRvVPvvt ∉∧∉∀ , (4h) Raw Material Requirements  −= e teevtv PbR ,,, VRvt ∈∀ , (4i) Capacity Constraints tetete CXYP ,,3,, ⋅≤ et,∀ (4j) Capacity, R&D, and Budget Stock Bounds tete CXCX ,1, ≤− et,∀ (4k) tete CRDCRD ,1, ≤− et,∀ (4l) CostmaxCRDCRD RCR CXCXCC e tete VRv tvtv e tetete ≤−+ + −    − ∈ − )( )(1000 1,, ,, 1,,, t∀ (4m) Additional constraints for technologies with multiple inputs and/or outputs ', , '' i i atin atin ii ii MWP MWP λ λ = ∀i,i'∈[1,m], ∀t∈[t,tmax] (4n) ', , '' j j cto cto jj jij MWP MWP ϑ ϑ = ∀j,j'∈[1,n], ∀t∈[t,tmax] (4o) ', , , ' ii tin tin i i P P θ= ∀i,i'∈[1,m], ∀t∈[t,tmax] (4p) ', , , ' jj to to q P P i i = ∀j,j'∈[1,n], ∀t∈[t,tmax] (4q) Stage-gate representation constraints tseste YLOCX ,,, ⋅≥ set ,,∀ (4r)  = −−⋅≤ 4 1 1,,, )( s sstsete HIHIYCX et,∀ (4s) tsetse YY ,1,,, −≤ { }{ }maxSsset ,,3,2|,, ∈∀ (4t) 1,,,, −≥ tsetse YY set ,,∀ (4u) Nomenclature Sets e∈E: technologies or pseudo-arcs v∈V: materials or pseudo-nodes v∈VP: products, VP⊆V v∈VR: raw materials, VR⊆V v∈VRR: renewable raw materials, VRR⊆VR⊆V ai∈Ae: input material nodes for tech e with m inputs, i∈[1,m], Ae⊆V cj∈Ce: output materials for techn e with n outputs, j∈[1,n], Ce⊆V de: delivering node for tech e with multiple input and/or output, de⊆V ini∈Ine: input transfer pseudo-arcs for tech e with m inputs, i∈[1,m], Ine⊆E maine: main arc of tech e with multiple input and/or output, maine⊆E oj∈Oe: output pseudo-arcs for tech e with n outputs, j∈[1,n], Oe⊆E re: receiving node for tech e with multiple input and/or output, re⊆V Parameters γv: annual increase rate for demands, ∀vVP kv: extraction cost coefficient of nonrenewable raw materials (∀vVRvVRR) IR: inflation rate Dv,0: initial product demands (∀vVP) ηe: yield of tech e λi and ϑj: stoic. Coef. of reactants and products θi,i': production ratio btwn input material ai and ai' qi,i': production ratio btwn input material cj and aj' Costmax: maximum annual budget available     − = otherwise 0 ogy by technol produced is material if 1 logy for techno material raw a is material if /1 , ev ev b e ev η Variables CCe,t: expansion cost of tech e at time t CRv,t: cost of nonrenewable material v at time t, ∀vVRvVRR CXe,t: cumulative capacity of tech e at time t CRDe,t: total R&D expenditure of tech e at time t Pe,t: production from tech e at time t Rv,t: amount of material v produced (∀vVP) or consumed (∀vVR) at time t TC: total cost Xe,t: capacity expansion of technology e at time t    = otherwise 0 at time stageat least at is y technologif 1 ,, tse Y tse 1329 4. Linear relaxations for the bilinear and nonlinear terms Two different approaches are used for obtaining linear relaxations of the bilinear terms: (1) linearly segmented tight relaxations (Misener et al., 2011), and (2) radix-based relaxations (Kolodziej et al., 2013). For the bilinear equation z =x.y where x,y∈R, xL ≤ x ≤xU, and yL ≤ y ≤ yU, the linearly segmented tight relaxation approach uses Np number of partitions in one variable, say x, and introduces concave and convex linear constraints above and below the original surface, z, in each partition, np, whose length is equal to a. A binary variable takes the value of one if partition np is selected, and is equal to zero otherwise. A continuous switch variable, Δynp, is used to incorporate the term (y-yL) if partition np is selected in the concave and convex envelope. The main principle behind radix-based relaxations is to approximate the bilinear constraint by representing one of the variables as a discrete value on a specified radix (commonly set to 10) to an arbitrary precision, defined with a range of 0 ×10p to 9×10P. To remove all nonlinearities from the original MINLP, the nonlinear terms defined in Section 3 should also be relaxed. These nonlinear terms are monotonically decreasing functions with zero asymptotes, and we propose using linearly segmented under and upper estimators for replacing these constraints. Here, the basic idea is to build supporting tangent lines at every breakpoint of each segment as an under estimator, and a straight line connecting the bounds and the function evaluation at the bounds as an upper estimator. Let )(xfy = with UL xxx ≤≤ be a monotonically decreasing function with zero asymptote, the formulation for the linearly segmented under estimator is given by Eq(5). Np xx a LU − = (5a)  == +≤≤−+ Np np np L Np np np L bnpaxxbnpax 11 ..).1.( (5b) 1 1 = = Np np npb (5c) ( ) ( )LL xx xfxx x xf y L +− ∂ ∂ ≥ = . )( (5d) ( ) ( )npaxfnpaxx x xf y LL npaxx L ... )( . ++−− ∂ ∂ ≥ += np∀ (5e) where Np is the total number of segments. The upper estimator for the function can then be defined via Eq(6). ( ) ( ) ( ) ( )LL LU LU xfxx xx xfxf y +− − − ≤ . (6) Once the bilinear terms are relaxed using one of the approaches and the nonlinear terms are replaced with under and upper estimators as discussed, the resulting formulation becomes a mixed integer linear program (MIP). We use the solution of this MIP as initialization for the original MINLP, and solve it using local solvers. The local solvers considered in this study are DICOPT and SBB. The solution obtained by the defined approach is compared to the solution of the original MINLP obtained by BARON to assess its quality. 5. Results and Discussion 5.1 Test problem The evolution of ethylene production from biomass (corn grain + corn stover) and naphtha in all US market for a planning horizon of 50 y is used as the test problem. In this simplified commodity chemical production system, biomass can be converted to ethanol via fermentation, which is followed by catalytic dehydration yielding ethylene. Alternatively, biomass may first be gasified to syngas, which is then converted to ethanol via catalytic conversion. Naphtha is cracked to yield ethylene. The details of this test problem can be found in Cremaschi (2011). The equations that contain the bilinear and nonlinear terms are modelled using the equations given in Section 3. The resulting MINLP, implemented in GAMS V23.6.2, has 5,952 constraints and 4,301 decision variables. It is solved globally using BARON with relative optimality gap 0.1 %. The solution is 30,927.0 billion US dollars, and it was obtained in 497.57 CPU seconds on a computer with Windows 7 Professional 64-bit operating system, Dual Intel E5405 2.0 GHz processors, and 8 GB RAM memory. The local solver DICOPT yields a solution of 30,924.8 billion US dollars in 1.84 CPU seconds when initialized at Ye,s,t = 1, CXe,t = 30, Xe,t = 1, RDe,t = 1, CRDe,t = 1, CCe,t = 1, Pe,t = 10, CRv∈VR,t = 1, and Rv,t = 100. 1330 5.2 Results with linearly segmented tight relaxations of bilinear terms (LSTRBT) In this section, the results presented are obtained by first solving the MIP using CPLEX to a relative optimality gap of 0.1 %, and then using this solution as the initialization for DICOPT and SBB. The number of variables and constraints of the MIP are 4,301+1,950Np and 13,152+1,350Np, respectively. The maximum number of segments considered in this study is 15. Three performance metrics measures the quality of the solutions: (1) the percentage difference between the MIP objective function value and the MINLP objective function value obtained using the MIP solution initialization, (2) the comparison of the obtained objective function value to the BARON solution, and (3) the overall computational time necessary to obtain these solutions. Figure 2 shows plots of the performance metrics with increasing number of segments when DICOPT and SBB are used as the local solvers for the MINLP. Figure 2(a) shows how the percentage difference between MIP solution and the corresponding MINLP solution changes with the number of segments. This percentage difference is calculated as (MIP-MINLP)*100/MIP. Figure 2(a) reveals that for both solvers the difference decreases slightly with the increase in number of segments with a decreasing gradient. The plot also suggests that increasing the number of segments beyond five does not impact the difference significantly. Figure 2(b) plots the normalized objective function values (Total Cost) versus the number of segments. Here, the normalized total cost is calculated by dividing the total cost obtained by the proposed approach by the objective of the original MINLP solution by DICOPT (30,924.8 billion US dollars). The lower and upper bounds of the objective function obtained by BARON are also shown in Figure 2(a). The plot suggests that the solution obtained by the proposed approach are not sensitive to the number of segments, and they are relatively close to the solution obtained by DICOPT and within the bounds of the BARON solution. The computational cost in CPU seconds is plotted against the number of segments in Figure 2(c). The solution time for the original MINLP with BARON is included as the dashed line for reference. When DICOPT is used as the local solver, the proposed approach is faster than BARON with a similar solution quality up to 11 segments. For SBB, the corresponding number of segments is nine. (a) (b) (c) Figure 2: The change in (a) difference between MIP and MINLP solutions, (b) normalized cost, and (c) computational cost obtained with LSTRBT and linear under and upper estimators for nonlinear terms 0.00% 0.05% 0.10% 0.15% 0.20% 0.25% 0.30% 0.35% 0.40% 0.00% 0.02% 0.04% 0.06% 0.08% 0.10% 0.12% 0.14% 0.16% 0 5 10 15 P e rc e n ta g e d if fe re n c e b e tw e e n M IP a n d M IN L P , S B B P e rc e n ta g e d if fe re n c e b e tw e e n M IP a n d M IN L P o b je c ti v e s , D IC O P T Number of segments DICOPT SBB 0.9989 0.9991 0.9993 0.9995 0.9997 0.9999 1.0001 1.0003 1.0005 0 5 10 15 T o ta l C o s t (n o rm a li ze d ) Number of segments MIP + DICOPT DICOPT MIP + SBB 1.003 BARON Lower Bound BARON Upper Bound 0 1000 2000 3000 4000 5000 6000 0 5 10 15 S o lu ti o n T im e ( C P U s e c o n d s ) Number of Segments MIP +DICOPT BARON MIP + SBB ~ 65000 1331 5.3 Results with radix-based relaxations of bilinear terms The idea in this approach is to generate the MIP that contains the radix-based relaxed bilinear terms, and linearly segmented under and upper estimators of the nonlinear terms. The solution of this MIP will be used as the initialization for the local solvers. Table 1 lists the range of power (for Radix 10) used for each variable that is discretized for each bilinear term of the BTCC investment planning problem. The resulting number of constraints and variables for the ethylene production case study are 58,252+250Np and 86,301+500Np, respectively. Even with a single segment for obtaining the linear upper and under estimators for the nonlinear terms, the CPLEX did not yield a solution after 48 h due to the large number of constraints and variables. These results suggest that radix-based relaxation approach is not appropriate for this problem. Table 1: Power range for the discretized variables of the bilinear terms Bilinear term Variable p P BLTerm1e,t Xe,t -3 0 BLTerm2e,t Rv,t -3 2 BLTerm3e,t NLTerm2e,t -3 2 BLTerm4e,t CXe,t -3 0 6. Conclusions This paper investigated efficient solution approaches for the BTCC investment planning problem. The two sources of nonconvexities, bilinear terms and nonlinear equality constraints, are relaxed using different approaches. Linearly segmented tight relaxation and radix-based relaxation approaches are considered for bilinear terms, and linearly segmented under and upper estimators are proposed to replace the nonlinear terms. The MIP solutions obtained with different number of segments are used as initializations for solving the original MINLP with DICOPT and SBB, local solvers. The results recommend replacing the bilinear terms with linearly segmented tight relaxations and the nonlinear terms with linearly segmented under and upper estimators with fewer than five segments, solving the resulting MIP to a 0.1 % gap, using DICOPT as the local solver. Acknowledgements The authors greatly acknowledge the financial support provided by the National Science Foundation CAREER Grant No. 1055974. References Cremaschi S., 2011, Modeling next generation feedstock development for chemical process industry, Computer Aided Chemical Engineering, 29, 1040–1044. Fahmi I., Nuchitprasittichai A., Cremaschi S., 2014, A new representation for modeling biomass to commodity chemicalsdevelopment for chemical process industry, Computers and Chemical Engineering, 61, 77-79, DOI: 10.1016/j.compchemeng.2013.10.012. Kolodzieja S.P., Grossmanna I.E., Furman K.C., Sawaya N.W., 2013, A discretization-based approach for the optimization of the multiperiod blend scheduling problem, Computers and Chemical Engineering, 53 (11 June 2013), 122–142, DOI: 10.1016/j.compchemeng.2013.01.016. Kouvaritakis N., Soria A., Isoard S., 2000, Modelling energy technology dynamics: Methodology for adaptive expectations models with learning by doing and learning by searching., International Journal of Global Energy Issues, 14 (1-4),104–115. Misener R., Thompson J.P., Floudas C.A., 2011, APEGEE: Global optimization of standard, generalized, and extended pooling problems via linear and logarithmic partitioning schemes, Computers and Chemical Engineering, 35 (5), 876–892, DOI: 10.1016/j.compchemeng.2011.01.026. Peres A.P.G., Lunelli B.H., Filho R.M., 2013, Application of biomass to hydrogen and syngas production, Chemical Engineering Transactions, 32, 589–594, DOI: 10.3303/CET1332099. 1332