Short Term Planning Framework for Enterprise wide Production

54 secs ago - Download Citation · Email a Colleague · Order Reprints · Rights & Permissions · Citation Alerts · Add to ACS ChemWorx. SciFinder Subscri...
0 downloads 0 Views 4MB Size
Subscriber access provided by Bibliothèque de l'Université Paris-Sud

Process Systems Engineering

Short Term Planning Framework for Enterprise wide Production & Distribution Network of A Cryogenic Air Separation Industry Shamik Misra, Divya Saxena, Mangesh Kapadi, Ravindra D. Gudi, and Rachakonda Srihari Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b05138 • Publication Date (Web): 06 Nov 2018 Downloaded from http://pubs.acs.org on November 6, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 63 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Short Term Planning Framework for Enterprise wide Production & Distribution Network of A Cryogenic Air Separation Industry †

Shamik Misra,

Divya Saxena,





Mangesh Kapadi,

R.Srihari

†Chemical

Ravindra D. Gudi,

∗,†

and



Engineering Department,Indian Institute of Technology, Bombay, Mumbai, India ‡Praxair India Private Limited, Bangalore, India E-mail: [email protected]

Abstract

In this paper, a novel mixed integer linear programming model is proposed for short term planning of an enterprise wide production & distribution network for air separation industries. The objective is to develop a framework for the planning level, that is based on a more realistic abstraction of the scheduling level complexities, so as to evolve optimal implementable targets to the scheduling level. The approach is also oriented towards promoting a tighter coordination between the plants in an enterprise and to make intelligent use of the distribution network so as to minimize the overall production & distribution cost for the entire enterprise. The model involves intelligent formulation techniques to curtail the number of binary variables which in turn makes the model computationally ecient. A non-uniform time discretization framework is also adopted 1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

and the proposed framework is validated on a case study representative of typical industrial size production & supply chain network of air separation enterprises. Instead of rigorous modeling, the planning framework optimizes an industrial size test case to the required accuracy in stipulated time. This paper also outlines a decomposition based framework to achieve enterprise wide optimization in air separation industries. Cryogenic Air Separation, Enterprise wide planning, Integration of Production & Distribution Optimization, Non-uniform Time Discretization. Keyword:

1 Introduction Cryogenic air separation is amongst the most eective and ecient ways to separate the components of air viz. oxygen, nitrogen, argon in their gaseous and liquid form to a high purity. 1 However, the cryogenic air separation process is very energy intensive and to survive in the competitive business environment, air separation industries need to adopt an enterprise wide optimization(EWO) approach, which helps in reducing the global operational costs by integrating various decision making layers involved in manufacturing as well as distribution. 2 Enterprise wide optimization facilitates information exchange between dierent layers of the automation hierarchy and helps in making more optimal production decisions by integrating various aspects of the supply chain. EWO reduces the chances of getting suboptimal or even conicting decisions between manufacturing and distribution layers. There are several facets of enterprise wide optimization ranging from planning, scheduling to real time optimization, 2 and several applications in dierent chemical sectors have been explored in literature. 3,4 However, this paper mainly focuses on the planning layer and its integration with the scheduling layers to achieve enterprise wide optimization of an air separation industry. Manenti and Rovagilo(2013) 5 apply EWO for a case study involving two air separation plants(ASP) and show the ecacy of EWO in increasing the exibility of the whole enterprise, which is unattainable if individual plants are optimized separately. Manenti and 2

ACS Paragon Plus Environment

Page 2 of 63

Page 3 of 63 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Rovagilo(2013) also list certain peculiarities regarding of ASP operational procedures (such as common customers, shared inventory, dierent delivery procedures for gaseous and liquid products, dependence on contractual obligations regarding inventory capacities and power consumption) which weigh down individual plants exibility. However, EWO can address these issues and can help in minimizing the operational costs by reaping benets from some of the contracts. Though the case study presented by Manenti and Rovagilo(2013) only considered two ASUs and mainly focused on various contractual agreement, their approach did not include any distribution aspects and its impact on the overall protability of the enterprise. Their approach however did assess the dierence between plant wide and enterprise wide optimization and motivated the EWO approach towards achieving improved exibility. Zhang

et al.(2016) 6 listed some recent attempts towards achieving EWO for industrial de-

mand side management, and mentioned the following challenges that need to be addressed in the process of implementing EWO in an industrial framework: 1) Proper modelling of the process limitations and exibilities, 2) incorporating constraints related to energy requirement and costs into the production optimization framework, 3) decision making across multiple time and spatial scales etc. While addressing the above challenges is mandatory, an additional requirement is that the production and distribution aspects also need to be jointly considered to eectively realize enterprise wide objectives. Production cost minimization involves proper selection of plant operation modes (including slates or operation points) 79 at the right time of the day to derive benets from time of day power pricing. In other words, production cost optimization aims to minimize the power consumption and also minimizes the power bill at minimum power consumption, while meeting the required product demand. On the other hand, distribution cost optimization involves proper selection of routes, proper assignment of vehicles to selected routes, minimum distribution eort without any stock out or late delivery. It is important to note that if the production optimization problem is solved independently, then the solution obtained reects the minimum power cost. However, it does not guarantee the minimum distribution cost; similarly, optimization of distribution as3

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

pects,separated from production limitations, leads to suboptimal solutions. Hence, to arrive at the global optimal solution, the trade-o between production and distribution costs need to be considered rather than minimizing these two costs separately. To tackle the resulting combinatorial nature of the temporal and spatial complexities a number of decomposition based approaches have been proposed earlier in literature. 10,11 An air separation industry typically has substantial number of plants and large number of associated customers and a decomposition based scheme that promotes a tight coordination between the two layers needs to be adopted to eectively exercise EWO in such an enormous and complex setup .A typical decomposition based scheme that promotes a tight coordination between the two layers is outlined in this paper, and is illustrated in Figure 1. The production targets (such as: production at each plant, demand contribution of each primary and secondary inventories) and distribution targets(such as: liquid pickup targets at each inventories, vehicle and inventory route selection, extra vehicle purchase) are sent to the individual scheduling layer where more rigorous optimization of production and distribution aspects are exercised to meet the abovementioned targets. Inventory Targets

Gas & Liquid Demand

Power Availability & Cost

EWO : Enterprise Wide Production & Distribution Planning(EPDP) Layer

Production Target

Distribution Target Distribution Scheduling & Tour Optimization

Single Plant/Enclave Production Scheduling

EWO : Scheduling Layer Figure 1: Decomposition Framework to Achieve EWO in Air Separation Industry

4

ACS Paragon Plus Environment

Page 4 of 63

Page 5 of 63 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

However, to achieve tight integration via decomposition approaches, the main bottleneck as described by Grossmann(2005), 2 lies in establishing proper framework at the upper level problem so that it provides achievable targets to the lower level. As the upper layer encompasses the entire enterprise, intelligent formulation techniques are needed to investigate every complex trade-os with minimum computational eorts. The lower level frameworks,which involves 1) single plant 79,12,13 or enclave scheduling 14 and 2) distribution scheduling and inventory routing, 15 are already available in literature. So the main objective of this paper is to provide a suitable framework for the planning layer where both production and distributing aspects of an air separation industry are represented & optimized simultaneously. While earlier work of Marchetti

et al.(2014) 16

does seek to achieve simultaneous produc-

tion and distribution optimization for industrial gas operations their proposed framework was substantially simplied and did not represent the detailed production and distribution constraints in a rigorous manner. The relatively simplied abstraction could lead only to a conservative coordination, which could perhaps be sub-optimal or even unattainable at the scheduling level. In the approach presented in our paper, we seek to provide a novel framework that provides for a realistic abstraction of the detailed production and distribution constraints in the planning layer. The framework presented in our paper, owing to the use of intelligent formulation techniques, as would be presented latter, is computationally ecient while being suciently rigorous. The results presented in the following sections highlights that the computational time and achieved optimality gap both are signicantly lower than earlier reported instances(Marchetti et

al.(2014)) on a problem that is signicantly more rigorous than the reported one. Zamarripa et al.(2016) 17 extended the formulation presented by Marchetti et al.(2014) and adopted a rolling horizon framework to facilitate computational eciency. However, their formulation also carry forward the same abstraction & assumptions made by Marchetti et

al.(2014) which may lead to sub optimal or infeasible solutions. Zhang

et al.(2017) 18 highlighted the importance of having an detailed and accurate representation of the critical operational aspects in the planning layer. 5

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

In the paper presented by Zhang

et al.(2017), 18

a multi grid framework was proposed to

address production routing problems and a heuristic based iterative solution method was employed for solution of the large instances. However, some important complexities in production of an ASP such as (i) eect of multiple transition modes, (ii) shared inventories,  (iii) common customers etc. as well as in distribution such as (i) secondary distribution network, (ii) dierent type & size of vehicles, (iii) Intra plant product exchange network,  (iv) pickup from supply sources etc. are not considered. Zhang

et al.(2017), 18

attempted heuristic based product and vehicle routing in their pre-

sented framework.In this paper, we propose a novel concept of time capacity at the planning layer, the inclusion of which provides optimal decisions regarding inventory, customer & vehicle suitability, and also provides realistic abstraction regarding the number of trips the vehicles need to make to fulll all the demands. The proposed formulation also provides an estimate on the number of vehicles of dierent sizes that need to be purchased/hired in a distribution region to successfully quench all the distribution needs. The abovementioned optimal decisions are sought to be cascaded to the lower layer where a detailed routing can be done in a manner as part of the decomposition strategy illustrated in Figure 1. The novel framework presented in this paper, incorporates all the complexities mentioned above and presents a rigorous and computationally ecient model formulation which is able to solve signicantly large scale instances compared to the case studies published in literature 16,18 with relatively smaller optimality gap in comparatively less computational time. The remainder of the paper is organized as follows: Section 2 briey discusses the major production and distribution aspects of a typical a multi site air separation network and enlists the possibilities as well as the associated complexities. Section 3 discusses the proposed MILP framework to achieve enterprise wide production and distribution planning. The ecacy of the framework is highlighted through illustrative and industrial case studies in Section 4, followed by concluding remarks in Section 5.

6

ACS Paragon Plus Environment

Page 6 of 63

Page 7 of 63 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

2 Problem Denition Air separation industry is a very important ancillary to the steel industry and a large enterprise with extensive reach in the market usually has several plants span over the whole country. As mentioned earlier, the cryogenic air separation process separates components of air in their gaseous and liquid form at a very high purity. The aforementioned gaseous and liquid products are two major business streams of air separation industry. These two business lines are tightly coupled with each other from an operational point of view as the products are obtained simultaneously from the same sets of production units. The onsite and merchant liquid business involves production and sale of oxygen, nitrogen, argon in their gaseous and liqueed respectively. The gaseous products are consumed by the onsite customer plants while liquid products are delivered to merchant liquid customers using distribution networks. The merchant liquid distribution network in a typical geographical spread can be categorized into two parts: primary and secondary distribution networks. The primary distribution network involves transportation of liquid products from plant locations to customers as well as to secondary storage locations. The products can also be distributed to customers from secondary storage locations, which is called as secondary distribution. The primary as well as secondary distribution networks use vehicles(trucks) of varying capacities (ranging from 3-30 tonnes) to supply liquid products to customers. The production of gaseous and liquid products is limited by plant capacity. Furthermore, in most of the cases, the demand for gaseous product drives the functioning of production plant and thereby aects the availability of liquid products at a particular plant location. The plants in an enterprise may vary in capacity, number and type of units. The end products may also vary across plants. In addition, the design of the plant inuences various operation modes in which the plant can operate and the power consumptions at each of these modes. The cost of power (which also can vary with location, time of day and contractual agreements) decides the cost of production. Similarly, the cost of distribution is constrained by the selection of routes, vehicle types and 7

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

capacity, transit time and distance etc. Simultaneous optimization of production and distribution aspects using the enterprise wide production and distribution planning tool described below helps constructing the following business decisions, ˆ Production and distribution plan for at least a full week. ˆ Inventory prole for plants, secondary storage and customer locations. ˆ Gaseous demand fulllment. ˆ Global (within and across country) liquid demand distribution and fulllment. ˆ Achieving globally minimum (electricity and distribution) cost reaping benet of all available energy and distribution contracts. ˆ Eective utilization of distribution capacity and routes to optimize distribution costs. ˆ Exploitation of the synergy between plants, which will help to make more informative decisions regarding the plant operations and also will help to build more resilient supply chain. The proposed approach in our paper seeks to achieve all of the above decision making through a tighter coordination between the production and distribution layers. Now, before describing the proposed mathematical framework in details, a brief discussion on the cryogenic air separation process and the associated plant topology is provided here.

2.1 Brief Discussion on the cryogenic air separation process An air separation plant generally contains the following types of units: 1) air separation unit (ASU), 2) liqueer, 3) compressors, 4) driox units and 5) Vent units. The process network for a standard cryogenic air separation plant is illustrated in Figure 2. Air is liqueed in ASU and the components such as, oxygen, nitrogen, argon are separated by cryogenic distillation method. The above-mentioned components are recovered with high purity in their liquid 8

ACS Paragon Plus Environment

Page 8 of 63

Page 9 of 63 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

and gaseous form in this method. Gaseous oxygen(GO2) and gaseous nitrogen (GN2) are directly produced in the ASU. However, Gaseous argon(GAr) is produced by a process called drioxing, 7,9 by which liquid argon are vaporized into its gaseous form to cater to the onsite customers. Other products(such as liquid oxygen(LN2) and liquid oxygen(LO2)) can also be drioxed to compensate for the inadequate production of the gaseous counterparts. A more detailed account regarding drioxing process and its eect on reducing overall production cost can be found in Misra

et al.(2017). 9

Gaseous products can further be compressed in dierent compressor units to meet customer VMPGN2

VGN2

VGO2

VentMPGN2 VentGO2

OFF

ON

ON

VentGN2

OFF

ON

OFF

LMCompGN2 GO2

MPGN2

OFF

ON

GN2

OFF

ON

MHCompGN2

ON

OFF

LiqGN2

Startup

NoAr On

Startup

AS ON

Shutdown

LHCompGN2

Shut down

ON

HPGN2

OFF

OFF

DrioxN2 LN2

OFF

ON

LO2 OFF ON

OFF

DrioxO2

LAr VGAr

ON

DrioxAr

GAr

ON

OFF

VentHPGN2 VHPGN2

Figure 2: Process Network with Operation Modes of an ASP specications. As can be seen from Figure 2, gaseous nitrogen are passed through three dierent compressor units to produce medium pressure GN2 (MPGN2) and high pressure GN2 (HPGN2). If the liquid products available in the inventory is not enough to meet the liquid product demand, gaseous products can also be liqueed using the liqueer unit(the unit LiqGN2 in Figure 2 is one such liqueer unit which transforms GN2 to LN2). Each of these units can operate in dierent operational modes which is also illustrated in Figure 2. The process network shown in Figure 2 constitutes a rigorous representation of the 9

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

actual process, and minute nuances of this process are incorporated in the detailed scheduling frameworks presented by Misra

et al.(2017,2018). 9,14 The main bottleneck in extending

the abovementioned rigorous single plant/enclave formulations to simultaneously optimize productions in a large number plants leads to combinatorial explosion. With increasing number of constraints and variables, the branch and bound problem size increases exponentially making the problem computationally very expensive, if not intractable. To make the planning framework computationally ecient while capturing the intricacies of production and distribution aspects with sucient rigor and accuracy, certain assumptions and abstractions are made which are enlisted in below with proper justication.

2.2 Important Model Characteristics & Assumptions Choice of units 1. Gaseous product compressors are accessory units and are only required to further compress the product as per customers specication. Post production processing units such as gaseous product compressors can be excluded from the planning layer as they do not aect the overall production. 2. On the other hand, driox units need to be a part of the planning layer as drioxing decisions provides an important exibility to the optimizer to nd the minimum power cost by reducing the production at peak power cost hours.

Operation modes 1. In the planning layer the transition modes such as startup and shutdown, have been lumped with the unit downtime(i.e. `OFF' mode). Since the product ow during the above mentioned transition modes are not of desirable purity, the calculated production target will still be realistic & will not be aected due to the above mentioned approximation. 10

ACS Paragon Plus Environment

Page 10 of 63

Page 11 of 63 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

2. However, the power consumption during startup and shutdown cannot be ignored and therefore are represented in the formulation. 3. Another transition mode `NoArOn' represents the excess time needed for liquid argon to reach its required purity limits. During `NoArOn' mode production of LO2 and LN2 are as similar as of `ON' mode but there is no liquid argon production. This operation mode needs to be meticulously captured in the planning layer as the production and power consumption during this operation mode are signicant. More details about 'NoArOn' mode is presented in the following section.

Gaseous product demand In the planning layer, the compressed gaseous products are not considered as separate products. For instance, if the customer quotes demands for HPGN2 and MPGN2, then in the planning layer proposed in our paper, these two demands will be added & will be considered as GN2 demand.

Liquid Product Demand 1. Theoretically, the liquid product demands can be treated as global, i.e. any plant in the enterprise can contribute to fulll the liquid demand quoted by any of the customer. Though the above scenario is theoretically possible, it is not feasible in reality as such liquid product exchange will not be optimal if the customer is far away from the plant. However, the optimizer will still consider this as a possible scenario, though not feasible in reality, and engage computational resources to nd out its feasibility, which will make the algorithm computationally expensive. 2. The organization can easily decide the compatibility between the customer and plant depending on (i) the feasible distribution radius, based on the experience and an (ii) an indicative calculation of the probable distribution costs based on experience. Based on

11

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the above two factors, a customer plant suitability matrix CMc,m can be envisaged and provided as an input to the algorithm, which indicates the feasible mapping between the customer & plant.

Distribution region 1. The entire distribution network is divided in multiple regions r depending on liquid source and customer locations. The idea behind creating this segregation is to guide the optimizer to select vehicles which are available nearby the sources. 2. The organization can choose these regions based on geographical location or operational practices. Vehicle availability will be quoted based on the region with the assumption that each and every liquid source in that region can use any of those vehicles available in that region. However, the vehicles used to transport the products, can travel from a region to another depending on the customer plant suitability.

Vehicle types 1. Vehicles used in air separation industry for intra country distribution are generally delivery trucks. In this paper, the delivery trucks are divided into two groups based on their size, i) small & ii) large. One can also categorize the trucks into several groups based on their size, transfer mechanisms required during loading & unloading operations. 2. The reason behind this segregation is that not all the customers can receive liquid products from all types of trucks. Customer site constraints such as availability of proper roads, parking space & other constraints related to the availability of suitable unloading mechanisms may preclude some truck options. 3. The same suitability also applies to the other supply sources. In our paper it is assumed that all the plants and secondary storages are suitable to all kinds of available vehicles 12

ACS Paragon Plus Environment

Page 12 of 63

Page 13 of 63 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

considering standardization. However, if needed one can easily implement a plant and vehicle suitability in a similar way.

Transit time & costs 1. The travel time and cost between the inventories or among liquid sources and customer are normally calculated based on the distance and average speed of a truck with certain capacity. 2. However in our formulation, the actual transit time is calculated considering loading & unloading time, waiting time, curfew hours in addition to the travel time to realistically represent the overall transit time. Transit time between one inventory to another, customers and supply sources are calculated, respecting vehicle compatibility with the source and the destination. 3. Furthermore, in our work the transit cost is taken as twice the actual cost to account for the to & fro visit.

Product distribution/transfer mechanisms 1. There are three types of possible liquid product movement in a region: 1) Plant or secondary storage inventory to customer, 2) other supply source to plant or secondary storage inventory & 3) product exchange amongst the inventories. 2. we look at the possible distribution aspects i) single drop & ii) multi drop. 3. If we just account for the single drop scenario then we will grossly overestimate the distribution cost & it is impossible to point the exact multiple drop scenario without delving into detailed routing calculation. Incorporating detailed routing decisions in a enterprise wide production and distribution framework makes the problem intractable.

13

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

4. To account for both the situation in a more simplistic manner, we introduce a novel concept of time capacity, which is explained in detail in the following section. It needs to be noted that all the above-mentioned assumptions and abstractions are mainly utilizes the inherent exibility of the process operations that can be exploited to make a computationally ecient framework, without compromising with the important production and distribution decisions, the idea of which can easily be extended and adopted for other manufacturing processes. Some of the aforementioned characteristics and assumptions are reiterated in the formulation section to give the constraints more contextualize representation.

3 Enterprise wide Production & Distribution Planning Model Formulation As mentioned earlier, model formulation at the planning level needs to be suciently rigorous to generate realistic targets to the lower layers. In the sequel we pose these rigorous constraints in a MILP based framework to represent the problem. It is to be noted that all the variables presented in the following MILP formulation should be assumed as positive variables unless specied otherwise. Details of all the indices, sets, variables and parameters are provided in the nomenclature section along with the descriptions presented in the text.

3.1 Non-uniform Time Discretization In the following a non-uniform discrete time framework is adopted to represent the process constraints and limitations of production and distribution of air separation industry. Two dierent time slots are adopted. The ner grids are used to represent the detailed nuances of the production aspects such as 1) unit operational modes and their transitions, 2) power consumption 3) production of liquid and gaseous products 4) Gaseous product demand ful14

ACS Paragon Plus Environment

Page 14 of 63

Page 15 of 63 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

llment etc, whereas, the coarser time grids are used to represent calculation regarding liquid product inventory management and other distribution aspects. Several factors determine the choice of the discretization time & the possibility of designing a non uniform time discretization framework. The choice of the discretization units demands an in depth knowledge of the process so that it does not overestimate or underestimate any operational nuances. For example, if the peak time for the grid power source spans a minimum duration of 4, the smallest unit of the time scale will be 4 hours for those plants which use grid power source as an option. The choice of the smallest unit also depends on the demand periods. For example, the onsite demands for gaseous products are foretasted hourly based on the historical data and generally does not have a signicant change in every hour and can be averaged over 2-4 hour time window. On the other hand, the liquid demands are predened and the demand window spans over a day, and therefore for all the calculations regarding liquid inventory management and distribution aspects day based discretization can be adopted. Another signicant factor in choosing the smallest time discretization unit is the minimum uptime or downtime of the units in the plants. The time discretization unit should reect the minimum uptime or downtime of the units as closely as possible to decrease the chance of over or under approximation of power costs. Since the objective of the planning layer is to generate an achievable production target for the local scheduling and distribution network, the formulation level at this layer needs to represent suciently detailed but computationally ecient abstraction of the lower level nuances. In any case, a detailed scheduling will be done locally which will consider a more detailed description of the plants and will provide optimum schedules based on the production target. The eect of non-uniform discretization in reducing the computational eort was rst introduced by Velez & Maravelias (2013). 19 Misra

et al.(2017,2018)

applied the same non

uniform time discretization techniques for formulating single ASP 9 & enclave scheduling 14 frameworks to reap benets from dierent demand windows of gaseous and liquid products. For the case study discussed in the following, we have considered the smallest time unit as 15

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 63

4 hours. This could be expected to reduce the dimensionality of the variables as instead of 168(24X7) time slots, every variable will be calculated for 42(6x7) time slots over the same week. For plants where time of the day power pricing is not applicable, the time slot duration selection can be guided by the minimum uptime or downtime required for the safe operation of the units. It is preferred to have 4 hour time slot or at max 6 hour time slots for air separation plants to curb the chances of over/under approximation of the energy requirements. Inventory balance calculation and liquid demand distribution calculations are done on a day basis. The entire time horizon discretized equally on the basis of smallest grids are illustrated in Figure 3. In Figure 3, N T represents the length of the current planning horizon whereas,

-NPT

0

1

NT

Planning Horizon

Historical Horizon

Figure 3: Discrete Time representation of the entire horizon

−N P T..0 indicates the past time periods. The past data is needed to determine the status of the production units at the start of the current horizon. This data helps to determine the mode transitions, minimum up and downtime constraints for the production units at the onset of the current planning horizon. In case of distribution, past information regarding the distribution vehicles helps to calculate the exact quantity and time of availability of the distribution capacities. A similar equisized daily discretized framework is used for the calculations regarding liquid product inventory and distribution.

3.2 Allocation Constraints Throughout this paper index m ∈ M is assumed to represent the installations (plant or secondary storages). These installations are of two types, 1) plant (PL), 2) secondary storage 16

ACS Paragon Plus Environment

Page 17 of 63 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

(SS). Some plants may have larger inventory and those inventories can be used as secondary storage for other plants. These condition can easily be handled using inter plant liquid exchange framework, that is also a part of the presented framework. The production unit allocation constraints are only applicable to PL type installations. Not all the plants have all the units (e.g. non availability of Liqueer, at certain plants), and the availability of certain units(u) in a plant (m) are denoted by the parameter M Um,u . A binary variable yAlm,u,t is considered which assumes a value 1 if the unit is `ON'. If single/multiple units of a plant would be going under maintenance at a time period t then that condition is denoted by the parameter M Fm,u,t = 1. The constraint that represents a maintenance scenario is shown in Eq. 1.

yAlm,u,t = 0 ∀m(: M T ym 6= ‘SS 0 ), u(: M Um,u = 1), t(: t = 1..N T and M Fm,u,t = 1)

(1)

The most accurate estimate of the production and associated power consumption can be found by solving the rst principle models. However, the detailed rst principle based mechanistic models usually pose complex non-linearity and incorporating those models in a enterprise wide production and distribution planning framework makes the entire problem intractable. A simpler way is to nd a feasible operating region for every unit for each operation modes from the historical data or by solving the mechanistic models separately. The operation points enveloping the feasible regions are called as slates, and it is assumed that any point inside these feasible regions can be represented by a linear combination of the pre-calculated slates. Multiple feasible regions can also be calculated for a unit running on certain operation modes. A detailed description of the surrogate model and feasible region calculation can be found in Zhang

et al.(2016b). 20 Zhang et al.(2016a) 7 and Misra et

al.(2017) 9 applied convex region surrogate models for estimating production & power consumptions of every unit of an air separation plant. In Eq. 2, the formulation presented by Misra

et al.(2017) 9 is extended to calculate the coecients of all the slatesl selected at a 17

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 63

time period t for a unit u running in a certain operation mode in a plant m. L P

xSlCm,u,l,t = yAlm,u,t

l=1(:M U Lm,u,l =1)

(2)

0

∀m(: M T ym 6= ‘SS ), u(: M Um,u = 1), t(: t = 1..N T )

3.3 Transition Constraints As mentioned earlier, the transition modes such as `startup' or `shutdown' are not considered in the planning layer as separate operation modes. The power consumption during startup and shutdown modes are not been ignored and are accurately represented in later part of this formulation. Though the abovementioned transition modes are not represented as separate operation modes in this formulation, unit uptime and downtime constraints cannot be neglected. To indicate the transition event from `ON' to `OFF', the continuous variable

xShDwm,u,t will assume a value of 1. Similarly, the variable xStU pm,u,t takes on a value of one when a certain unit goes `OFF' to `ON'. The transition constraints are adopted from Hedman et

al.(2009). 21 The relation between the unit assignment variable and the transition

variables are captured using the Eqs. 3-4.

xStU pm,u,t − xShDwm,u,t = yAlm,u,t − yAlm,u,t−1 ∀m(: M T ym 6= ‘SS 0 ), u(: M Um,u = 1), t(: t = 2..N T )

ini xStU pm,u,t − xShDwm,u,t = yAlm,u,t − Alm,u,t−1 0

(3)

(4)

∀m(: M T ym 6= ‘SS ), u(: M Um,u = 1), t(: t = 1)

ini The parameter Alm,u,t indicates the operating status of the unit (`OFF' or `ON') at the start

of the current horizon. After a change in operation mode, as stated above, a unit has to stay in that mode for the minimum uptime period, to avoid any physical damage. Similarly, every unit also has a downtime period which indicates the time periods the unit needs to be in `OFF' mode. The minimum uptime and downtime of each unit u part of a plant m is 18

ACS Paragon Plus Environment

Page 19 of 63 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

min min stored in the parameter U T rm,u and DT rm,u respectively. The constraints to represent the

aforementioned limitations are modeled in Eqs. 5-6. t P

xStU pm,u,k ≤ yAlm,u,t

min +1) k=(t−U T rm,u

(5)

∀m(: M T ym 6= ‘SS 0 ), u(: M Um,u = 1), t(: t = 1..N T ) t P

xShDwm,u,k ≤ 1 − yAlm,u,t

min +1) k=(t−DT rm,u

(6)

0

∀m(: M T ym 6= ‘SS ), u(: M Um,u = 1), t(: t = 1..N T ) The above mentioned equations aptly represent the shutdown and startup transitions for every units. However, there is another transition mode `NoArOn' that needs to be accurately represented to properly estimate the liquid argon production. In a typical air separation unit(ASU), the time needed to reach deliverable purity for LAr is signicantly greater than the time needed for LO2 and LN2. To model this phenomena, Misra

et al.(2016,2017) 8,9

proposed addition of a new transition mode `NoArOn' for air

separation units during which the production of LO2 and LN2 are as similar as of `ON' mode but there is no liquid argon production. Though in the proposed EPDP formulation, transition modes are not considered in the planning level, the eect of of `NoArOn' condition cannot be ignored. The reasons are listed in the following: ˆ Reaching purity for LAr takes signicantly longer time than LO2 or LN2. ˆ During this transition mode though there is no production of LAr of desirable purity, LO2 and LN2 are produced with the same purity and quantity as in the `ON' mode. So, if we lump this transition mode with the 'OFF' mode, the resultant LO2 and LN2 production would be under approximated. ˆ Similarly, if the extended transition times required for LAr are ignored and only `ON' mode is considered then the resulting LAr production targets calculated in the planning 19

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 63

layer, would be unachievable by the single plant /enclave scheduling layers. To incorporate the `NoArOn' phenomena, without an introduction of additional operation modes, the following constraint sets(Eqs. 7-10) are formulated. t P

xN oArm,t −

xStU pm,u,td = 0

(7)

td=(t−N oArTm,u +1) 0

∀m(: M T ym 6= ‘SS ), u(: M Um,u = 1 and Asum,u = 1), t(: t = 1..N T )

xF Pm,u,p,t − xN oArm,t ∗ ArPmmax ) ≤ 0 ∀m(: M T ym 6= ‘SS 0 ), u(: M Um,u = 1 and Asum,u = 1), p(: LArFp = 1), t(: t = 1..N T ) (8)

xF Pm,u,p,t − xPm,u,p,t ≤ 0 ∀m(: M T ym 6= ‘SS 0 ), u(: M Um,u = 1 and Asum,u = 1), p(: LArFp = 1), t(: t = 1..N T ) 

xF Pm,u,p,t − xPm,u,p,t + ((xN oArm,t − 1) ∗ ArPmmax )

(9)

 ≥0

∀m(: M T ym 6= ‘SS 0 ), u(: M Um,u = 1 and Asum,u = 1), p(: LArFp = 1), t(: t = 1..N T ) (10) The variable xN oArm,t basically indicates the time needed for liquid argon to reach purity after ASU is turned `ON'. Eqs.(8-10) dictates that when variable xN oArm,t is 1, then for product LAr, the variable xF Pm,u,p,t assumes the same production quantity as indicated by the variable xPm,u,p,t .The variable xPm,u,p,t represents the quantity of product p produced in plant m ∈ M in unit u ∈ U at time period t. The dierence between the variables xPm,u,p,t and xF Pm,u,p,t is added as total produced quantity of liquid argon at time period t in the inventory balance equation. As at the time of `NoArOn' periods, xF Pm,u,p,t assume the same value as xPm,u,p,t , there will be no production of argon during those time periods which is evident from Eq.10. It should be noted that, the variable xF Pm,u,p,t is dened and calculated for liquid argon only, to capture the eect of prolonged purity reaching time needed for liquid argon compared to other liquid products.

20

ACS Paragon Plus Environment

Page 21 of 63 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

3.4 Mass Balance Constraints In an ASP, the components of air mainly oxygen, nitrogen and argon are produced in their gaseous and liquid forms by cryogenic air separation method. The production or consumption in unit u during time slot t are expressed as a combination of the production/consumption rates of the selected operation slates. The variable xPm,u,p,t indicates the total product/state produced/consumed in unit u during time t at plant m and are mathematically expressed in Eq. 11. L P

xPm,u,p,t =

xSlCm,u,l,t ∗ SlP rm,u,l,p ∗ (T Et − T St )

l=1(:M U Lm,u,l =1)

(11)

0

∀m(: M T ym 6= ‘SS ), u(: M Um,u = 1), p(: M U Pm,u,p = 1), t(: t = 1..N T ) Where, T St and T Et represents the start and end time of a time slot respectively.

3.4.1 Gaseous Product Demands and Exchange Onsite gaseous product demands are fullled via pipelines. One or more plants which are connected to the customer via pipelines can quench the gaseous product demand of that customer. In case of an enclave of plants, plants which are connected with each other using pipelines, can also exchange gaseous products among themselves( if protable). Now, the onsite gaseous demands needs to be fullled completely. As shown in Eq. 12, one or more plants can contribute in quenching the demands provided those plants are connected with the customer via pipelines. M P

ODc,p,t −

xODm,c,p,t = 0

m=1(:M T ym 6=‘SS 0 and M Cm,c =1) 0

0

0

(12)

0

∀c(: CT yc = G ), p(: P T yp = G and CPc,p = 1), t(: t = 1..N T )

Generally, gaseous product demands drives air separation plants operations, however, if the production of certain gaseous products is inadequate to meet the demand, liquid products 21

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 63

can be drioxed and sent to the customer. More details about the drioxing process and related constraints can be found in Misra

et al.(2017 & 2018). 9,14 The plants in an enterprise may

dier in their capacity, choice of units and even in the products. Some plant may only produce liquid products. The suitability between the plant and deliverable products are captured in the parameter named F P Fm,p . The mass balance for deliverable gaseous products are expressed in Eq.13, while Eq.14 denes the mass balance constraints for intermediate gaseous products. U P

xPm,u,p,t −

u=1(:M Um,u =1)

P

xODm,c,p,t

c=1(:M Cm,c =1)

P

+

xGExm0 ,m,p,t

md=1(:M M Pmd,m,p =1)

P



(13)

xGExm,m0 ,p,t = 0

md(:M M Pm,md,p =1)

∀m(: M T ym 6= ‘SS 0 ), p(: P T yp = ‘G0 and F P Fm,p = 1), t(: t = 1..N T ) U P u=1(:M Um,u =1)

P

xPm,u,p,t −

xGExm,m0 ,p,t

md(:M M Pm,md,p =1)

P

+

xGExm0 ,m,p,t = 0

(14)

md=1(:M M Pmd,m,p =1)

∀m(: M T ym 6= ‘SS 0 ), p(: P T yp = ‘G0 and F P Fm,p = 0), t(: t = 1..N T ) The parameter M U Pm,u,p = 1 indicates the product p that can be produced in unit u of plant m. Similarly, the parameter M U Pm,u,p = −1 denotes a set products that are consumed in unit u of plant m. P T yp is a parameter which indicates the type of the product. One has to note that the constraints representing unit allocations, mode transitions and the gaseous product mass balance (Eqs. 1-14), are only written for the installations which have associated production units. It is important to note that all these constraints are formulated at the smallest time slots in the horizon.

22

ACS Paragon Plus Environment

Page 23 of 63 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

3.4.2 Liquid Product Inventory Management As the demand for liquid products are quoted day wise, liquid inventory balance equations can also be expressed on a daily basis.Furthermore, liquid product balance equations are formulated for all types of installations as can be seen from Eq. 15 and Eq. 16. Eq. 15 describes the material balance for the inventories which are at plant location(PSI) and associated with the production units, and has product formation and consumption terms. Secondary storage inventories(SSI) as represented in Eq. 16 do not have any production / consumption components and receive products from plants. Secondary storages can also take part in fullling customer demands using secondary distribution networks. M P

xQi,p,d = xQi,p,d−1 +

N PT

U P

(xPm,u,p,t − xF Pm,u,p,t )

m=1(:IMi,m =1) t=1(T Dt =d) u=1(:M U Pm,u,p =1) M N U P PT P



xPm,u,p,t

m=1(:IMi,m =1) t=1(T Dt =d) u=1(:M U Pm,u,p =−1) S P

xLPs,i,p,d − xWi,p,d − xDFi,p,d

+

s=1(:SIPs,i,p =1) N D P

+

i0 =1(:d−LExT Ti0 ,i,p >0)

xLExi0 ,i,p,d−LExT Ti0 ,i,p −

I P

xLExi,i0 ,p,d

i0 =1(:IIPi,i0 ,p =1)

∀i(: IT yi 6= ‘SSI 0 ), p(: P T yp =0 L0 and IPi,p = 1), d(: d = 1..N D) (15) S P

xQi,p,d = xQi,p,d−1 +

xLPs,i,p,d − xWi,p,d − xDFi,p,d

s=1(:SIPs,i,p =1)

+

N D P

xLExi0 ,i,p,d−LExT Ti0 ,i,p −

id=1(:d−LExT Ti0 ,i,p >0)

I P

xLExi,i0 ,p,d

id=1(:IIPi,i0 ,p =1)

∀i(: IT yi = ‘SSI 0 ), p(: P T yp =0 L0 and IPi,p = 1), d(: d = 1..N D) (16) Where,

xQi,p,0 = Qini i,p

∀i, p(: P T yp =0 L0 and IPi,p = 1)

(17)

where xQi,p,d indicates the inventory level for a product p (where,p ∈ P (: P T yp =0 L0 )) at day d . The Variable xDFi,p,d in Eq. 15, indicates the amount of liquid product p decanted from inventory i at day d to satisfy demand from all the customers associated with that plant

23

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 63

or secondary storage. If required, liquid products can also be purchased from supply sources to fulll inventory target requirements, and the purchased quantities are indicated here by the variable xLPs,i,p,d . The variable xLExi,i0 ,p,d indicates the amount of product exchanged between inventory i and inventory i0 . The inventory quantity at the end of each time period should adhere to the maximum and minimum capacity limits, and is enforced by constraints described in Eqs.18-19.

min QCi,p ≤ xQi,p,d ∀i, p(: P T yp =0 L0 and IPi,p = 1), d(: d = 1..N D)

(18)

max ≥ xQi,p,d ∀i, p(: P T yp =0 L0 and IPi,p = 1), d(: d = 1..N D) QCi,p

(19)

The product inventory are also bound by the targets that are set based on contractual agreements and business practices such as onsite customer contracts, upcoming maintenance activity etc. The target constraints, as presented in Eq. 20,however are soft constraints and can be violated incurring penalty costs. Though positive and negative both violations are indicated by the variables xP Qpi,p,d and xP Qni,p,d respectively, only the negative violations(i.e. when the product inventory is less than the inventory target) are penalized in the objective function. Furthermore, it should be noted that among xP Qpi,p,d and xP Qni,p,d , only one variable can assume a nite value at any point of time and that is equal to the violated amount.

xQi,p,d + xP Qni,p,d − xP Qpi,p,d − QTi,p,d = 0 ∀i, p(: P T yp =0 L0 and T V Fi,p = 1), d(: d = 1..N D)

(20)

However, if in some cases violation of a target is not permissible according to the business contract, then the inventory must be above the target level at the specied time periods.

xQi,p,d ≤ QTi,p,d 0

0

∀i, p(: P T yp = L and T V Fi,p = 0), d(: d = 1..N D)

24

ACS Paragon Plus Environment

(21)

Page 25 of 63 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Equation 22 indicates that if there is some contractual obligation regarding the minimum liquid requirement in the inventory, the amount of liquid product in that said inventory should be more than the minimum requirement.

QM Ci,p ≤ xQi,p,d ∀i, p(: P T yp =0 L0 and IPi,p = 1), d(: d = 1..N D)

(22)

3.4.3 Liquid Demand Constraints In air separation industry, liquid product demands are generally of two types: 1) Regular demands and 2) Spot Demands. Regular demands are mandatory to be fullled whereas, spot demands are optional. Since spot demands earn a higher revenue, maximum fulllment of spot demands without aecting other limitations are desirable. One liquid customer can be catered to by multiple plants and also one plant or one secondary storage unit can cater to multiple customers. Eq. 23 indicates that the entire liquid product decanted from an inventory to meet demand on any day should equal to the total demand(regular as well as spot) from all the associated customers. If the quantity of a certain product, that is decanted from an inventory without violating other contractual constraints, is not adequate to meet the total demand assigned to that inventory, the residual can be fullled by purchasing from competitors. The mass balance around liquid demand fulllment are shown in Eq.23. C P

(xRDc,i,p,d + xSDc,i,p,d ) = (xDFi,p,d +

c=1

S P

xLP Ds,i,p,d−ST Ts,i,p )

s=1(:SIPs,i,p ) 0

0

∀i, p(: P T yp = L and IPi,p = 1), d(: d = 1..N D)

25

ACS Paragon Plus Environment

(23)

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 63

As mentioned earlier, multiple inventories can cater to a customer, which is captured in the regular demand fulllment constraints, that are illustrated in Eqs. 24-25. min Dc,p,y,q ≤

N D P

I P

 xRDc,i,p,d



i=1(:ICPi,c,p =1)

s e :d+CT Ti,c,p ≥CDTc,p,y,q and d+CT Ti,c,p ≤CDTc,p,y,q

d=1

(24)

∀c, p(: P T yp =0 L0 and CPc,p = 1), y(: Y T yy = ‘REGULAR0 , q(: q = 1..CCnc,p,y )

max ≥ Dc,p,y,q

I P

N D P

 xRDc,i,p,d



i=1(:ICPi,c,p =1)

e s and d+CT Ti,c,p ≤CDTc,p,y,q :d+CT Ti,c,p ≥CDTc,p,y,q

d=1

(25)

∀c, p(: P T yp =0 L0 and CPc,p = 1), y(: Y T yy = ‘REGULAR0 , q(: q = 1..CCnc,p,y ) As EPDP also includes the distribution aspects, the transit time between the plant and customer as well as plant and supply sources along with the distribution cost, also needs to be considered. The parameter CT Ti,c,p indicates the maximum transit time needed to transfer a certain liquid product from inventory i to customer c. Similar constraints can be written for spot demand fulllment. As spot sales are optional, the optimizer has additional exibility in the decision of spot demand fulllment. However, if a decision is taken in favour of spot demand fulllment and some minimum spot demand is mentioned then it is necessary to fulll the minimum demand at the least. A binary variable

ySSc,p,y,q is introduced which turns 1 if spot sale of product p is fullled. Constraints regarding spot demand, when the minimum limit is mentioned, are expressed in 26 and 27. N D P

min ySSc,p,y,q ∗ Dc,p,y,q ≤

 d=1

 xSDc,i,p,d

s e :d+CT Ti,c,p ≥CDTc,p,y,q and d+CT Ti,c,p ≤CDTc,p,y,q

∀c, p(: P T yp =0 L0 and CPc,p = 1), y(: Y T yy = ‘SPOT0 , q(: q = 1..CCnc,p,y )

26

ACS Paragon Plus Environment

(26)

Page 27 of 63 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

N D P

max ySSc,p,y,q ∗ Dc,p,y,q ≥

 xSDc,i,p,d

 d=1

s e :d+CT Ti,c,p ≥CDTc,p,y,q and d+CT Ti,c,p ≤CDTc,p,y,q

(27)

∀c, p(: P T yp =0 L0 and CPc,p = 1), y(: Y T yy = ‘SPOT0 , q(: q = 1..CCnc,p,y ) If there is no minimum spot demands are specied, then the binary variable ySSc,p,y,q is not required and the aforementioned constraints can be replaced by Eq. 28. N D P

max ≥ Dc,p,y,q

 xSDc,i,p,d

 d=1

(28)

e s and d+CT Ti,c,p ≤CDTc,p,y,q :d+CT Ti,c,p ≥CDTc,p,y,q

∀c, p(: P T yp =0 L0 and CPc,p = 1), y(: Y T yy = ‘SPOT0 , q(: q = 1..CCnc,p,y ) Liquid product can also be purchased from other sources (purchased from competitor plants, if protable) to meet demands and inventory targets. The variable xLP Ds,i,p,d indicates the liquid product p purchased from supplier s at d to quench liquid demand at inventory i.

3.4.4 Liquid Product Supply Source Limitation As mentioned earlier, if the production of liquid products at a certain time periods is not sufcient to meet liquid demands and inventory targets, the shortage quantity can be purchased from supply sources (competitor plants). Products received from supply sources generally have an upper limit on their availability as shown in Eq. 29. Furthermore, if the lower limit is also mentioned then the purchased quantity need to be more than the minimum amount. N D P s d=1(:d≥STs,p,q

&

!

I P e d≤STs,p,q )

xLP Ds,i,p,d + xLPs,i,p,d

≤ SLmax s,p,q

i=1(:SIPs,i,p =1)

(29)

∀s, p(: P T yp =0 L0 ), q(: q = 1..SCns,p and SLmin s,p,q = 0) It is obvious that purchasing liquid products for fullling liquid demands (including spot demands) or inventory targets is costly. Hence, liquid product purchase from competitor plants is generally discouraged. However, this option is necessary as it helps in maintaining

27

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 63

an undisrupted supply chain. There are several other important production aspects that need to be incorporated in the formulation such as 1) contractual limitations on liquid decantation from plant associated inventories(onsite contracts) 2) constraints related to drioxing and free liquid etc. These constraints are essential to include in the planning framework for an adequate representation of the process dynamics which helps in providing achievable production targets to the local optimization framework of individual plants & enclaves. However, as these constraints are discussed and adequately represented in Misra et

al.(2018), 14 those are not reiterated in this

paper to maintain brevity.

3.5 Power Consumption & Cost Calculation As in cryogenic air separation method, the main raw material air is free, the major element of the production cost is power cost. So, proper estimation of the power cost is imperative in the EPDP framework. The power consumed by the unit during 'ON' mode are represented as a linear combination of the power requirements of the selected slates. Though the production during the transition modes such as `startup' and `shutdown' are not considered, the power consumption in those transition periods are signicant enough and cannot be ignored. Eq. 30 estimates the power requirement of unit u in plant m at time slot t.

xEm,u,t =

L P

xSlCm,u,l,t ∗ SlP wm,u,l ∗ (T Et − T St )

l=1(:M U Lm,u,l =1)

+xSUm,u,t ∗ SuP wm,u + xSDm,u,t ∗ SdP wm,u

(30)

∀m(: M T ym 6= ‘SS 0 ), u(: M Um,u = 1), t(: t = 1..N T ) The total power required to run all the units should be less than the power available. This criteria is handled in Eq. 31. The parameter U tP wm,t indicates the amount of power needed to run the oce utilities at each installation locations m. Now, Eq. 30 is written only for those installation which has production units associated, however, Eq. 31 is for all the

28

ACS Paragon Plus Environment

Page 29 of 63 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

installations as secondary storage locations also needs utility power. U P

E P

xEm,u,t + U tP wm,t ≤

u=1

xEPm,e,t

∀m, t(: t = 1..N T )

(31)

e=1(:M Em,e =1)

Now, there may be more than one power sources associated with an installation with dierent availability limits. The power purchased from a suitable power source should less than the power available from that power source. At any point of time, not purchasing power from one or more sources is an option; however, if minimum purchasing limits are mentioned, then if purchased, the purchased power quantity needs to be greater than minimum limit of that power source to maintain supply network stability. The availability limits may be specied for a particular hour or for a period of time. When the availability limits are available on an hourly basis, they are accordingly converted into time slot limits taking into consideration the duration of the time slots. Hence,

yEm,e,t ∗ ELmin m,e,t ≤ xEPm,e,t

∀m, e, t(: ELmin m,e,t > 0 and t = 1..N T )

(32)

xEPm,e,t ≤ yEm,e,t ∗ ELmax m,e,t

∀m, e, t(: ELmax m,e,t > 0 and t = 1..N T )

(33)

If the minimum power availability limit for a power source associated with an air separation plant, is not quoted, then above set of Equations can be replaced with Eq. 34.

xEPm,e,t ≤ ELmax m,e,t

∀m, e, t(: ELmax m,e,t > 0 and t = 1..N T )

(34)

3.5.1 Power Cost In the proposed formulation, time of the day power pricing is assumed, irrespective of the power source. If power price is not varying with time of day, constant price is assumed for all time slots. The total power purchase cost for each plant m from a power source e can be

29

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 63

calculated using Eq. 35.

xECm,e =

N PT

xEPm,e,t ∗ EP Pm,e,t ∀m, e(: M Em,e = 1)

(35)

t=1

Several dierent power contracts such as 1) minimum take or pay, 2) discount on block purchase may be available at each installations and considering these contracts in the planning layer are utmost important to correctly estimate the production cost at each installations. However, constraints related to this contracts are already available in literatures 7,9,14 and is not reiterated in this paper to maintain brevity.

3.6 Distribution Capacity Constraints The distribution optimization in the planning level generates an optimum inventory routing solution and pick-up schedule indicating the liquid product pick-up target at each and every liquid source (such as plant, secondary storage location, or other supply sources), that respects the delivery date and transit time to the destination locations. As mentioned earlier, in the presented formulation, the entire distribution network is divided in multiple regions

r depending on liquid source and customer locations. Vehicle availability are quoted based on the region with the assumption that each and every liquid source in that region can use any of those vehicles available in that region. However, the vehicles used to transport the products, can travel from a region to another depending on the customer plant suitability. There are three types of possible liquid product exchange: 1) Plant or secondary storage inventory to customer, 2) other supply source to plant or secondary storage inventory & 3) product exchange amongst the inventories. Equation 36 calculates the distribution of product p using various compatible vehicle v that are sent to customer c from inventory i to quench the total demand for that product. Similar equations which indicates product quantity distribution among carrier vehicles are also written for supply to inventory transfer (Eq. 37) and inter plant/secondary location product

30

ACS Paragon Plus Environment

Page 31 of 63 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

exchange(Eq. 38). The parameter CP kT yi,c,p = ‘Y 0 indicates that the customer takes care of the product transfer aspects, so there will be no distribution cost. If CP kT yi,c,p = ‘N 0 then the vendor needs to track and rell customer inventory based on their consumption pattern(vendor managed inventory).



N D P

 xRDc,i,p,d + xSDc,i,p,d

d=1

V P



xCDisi,c,p,v = 0

v=1(:V CT yi,c,p =V T yv )

(36)

∀c, i, p(: ICPi,c,p = 1 and CP kT yi,c,p = ‘N 0 ) N D P

  xLPs,i,p,d + xLP Ds,i,p,d −

d=1

V P

xSDiss,i,p,v = 0

v=1(:V ST ys,i,p =V T yv )

(37)

∀s, i, p(: SIPs,i,p = 1 and SP kT ys,i,p = ‘Y 0 ) N D P

xLExi,i0 ,p,d =

V P

xIDisi,i0 ,p,v

v=1

d=1 0

(38)

0

∀i, i (: i 6= i), p(: IIPi,i0 ,p = 1) In the planning layer, since we seek to integrate the production and distribution, it is important to realistically represent the transit times and evaluate transit cost to achieve an overall optimum solutions. The travel time and cost between the inventories or among liquid sources and customer are normally calculated based on the distance and average speed of a truck with certain capacity. However in our formulation, the actual transit time is calculated considering loading and unloading time, waiting time, curfew hours in addition to the travel time to realistically represent the overall transit time. Transit time between one inventory to another, customers and supply sources are calculated, respecting vehicle compatibility with the source and the destination. Furthermore, in our work the transit cost is taken as twice the actual cost to account for the to and fro visit. The total transit time(travel + loading/unloading time) needed for a vehicle v for delivering product p at customer location

c and come back to inventory i are calculated in Eq. 39. In a similar way the total transit time needed for v to transfer product from supplier or exchange product amongst inventories

31

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 63

can be calculated and shown in Eq.40 and Eq.41 respectively. C = 2 ∗ (V CT Ti,c,v,p + P LuTv,p ) T rTi,c,v,p

∀c, i, p(: ICPi,c,p = 1 and CP kT yi,c,p = ‘N 0 ), v(: V CT yi,c,p = V T yv )

S T rTs,i,v,p = 2 ∗ (V ST Ts,i,v,p + P LuTv,p ) 0

(39)

(40)

∀c, i, p(: ICPi,c,p = 1 and CP kT yi,c,p = ‘N ), v(: V CT yi,c,p = V T yv ) I T rTi,i 0 ,v,p = 2 ∗ (V LExT Ti,i0 ,v,p + P LuTv,p )

∀c, i, p(: ICPi,c,p = 1 and CP kT yi,c,p = ‘N 0 ), v(: V CT yi,c,p = V T yv )

(41)

In the above framework the loading/unloading time are assumed only to be dependent on the type of vehicle and the product. However, it can also vary with the source and destination location can easily be incorporated if needed. As is well known, the optimal distribution routing from inventory/plant to customer location via variety of routing option is a very complex problem. In our work, we approach this complexity by introducing a novel concept termed as time capacity.

3.6.1 Overall Time Capacity Now in the above the total transit time needed to distribute a truckload of specic product

p from a liquid source i to a customer c in a vehicle v is calculated(Eq. 42). By dividing the abovementioned transit time with vehicle capacity for that specic product, the unit quantity distribution time can be calculated. If we multiply the above-mentioned unit quantity distribution with the quantity of the total amount of liquid product that needed to be transferred we get the total time that is needed to complete a specic liquid product distribution. the same calculation can be done for all kinds of liquid product transfers. Now, the time needed to complete a specic liquid product distribution in a region of all kinds (inventory-customer, supply source-inventory, among inventories) should be less than the total transit time available in the planning horizon. This limitation can be captured 32

ACS Paragon Plus Environment

Page 33 of 63 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

using the Eq. 42. S P

I P

S T rTs,i,v,p

s=1 i=1(:IRi,r =1) (SIPs,i,p =1)

+

V Qv,p

∗ xSDiss,i,p,v

+

C P

I P

C T rTi,c,v,p

c=1 i=1(:IRi,r =1) (ICPi,c,p =1) I P

I P i=1(:IRi,r =1)

!

i0 =1 (IIPi,i0 ,p =1)

I T rTi,i 0 ,v,p

V Qv,p

V Qv,p

! ∗ xCDisi,c,p,v

! ∗ xIDisi,i0 ,p,v

≤ 24 ∗ N D ∗ V Ar,v,p

(42)

∀r, p(: P T yp = ‘L0 ), v(: RV Pr,v,p = 1) The total transit quantity available for distribution of a product is a function of availability of compatible vehicles in that region and their capacity. It is to be noted that the aforementioned constraints take care of both quantity as well as time limitations for distribution.

3.6.2 Daily Distribution capacity Though the abovementioned constraints aptly captures the overall distribution time and capacity limitation for the entire planning horizon, it might overestimate the distribution capacity available in a region on a particular day. To calculate the daily distribution capacity balance, rst we need to calculate the amount of distribution capacity that are getting consumed in a day for transferring a certain products among the liquid source(inventory or supplier) to the consumer(customer or other inventories). The distribution of liquid products from inventory i by a certain vehicle v on a specic day d to a customer c is captured by the variable xV CDisDayc,i,p,v,d and calculated using Eqs. 43-44. Similar set of equations can be written for product exchange and purchase scenarios are represented using Eqs. 45-46 and Eqs. 47-48 respectively. V P

xV CDisDayc,i,p,v,d = xRDc,i,p,d + xSDc,i,p,d

v=1(:V CT yi,c,p =V T yv ) 0

∀c, i, p(: ICPi,c,p = 1 and CP kT yi,c,p = ‘N ), d(: d = 1..N D)

33

ACS Paragon Plus Environment

(43)

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

N D P

xV CDisDayc,i,p,v,d = xCDisi,c,p,v

d=1

Page 34 of 63

(44)

0

∀c, i, p(: ICPi,c,p = 1 and CP kT yi,c,p = ‘N ), v(: V CT yi,c,p = V T yv ) V P

xV IDisDayi,i0 ,p,v,d = xLExi,i0 ,p,d

v=1

(45)

0

∀i, i , p(: IIPi,i0 ,p = 1), d(: d = 1..N D) N D P

xV IDisDayi,i0 ,p,v,d = xIDisi,i0 ,p,v

d=1

(46)

∀c, i, p(: IIPi,i0 ,p = 1), v V P

xV SDisDays,i,p,v,d = xSDs,i,p,d

v=1(:V ST ys,i,p =V T yv )

(47)

∀s, i, p(: SIPs,i,p = 1 and SP kT ys,i,p = ‘Y 0 ), d(: d = 1..N D) N D P

xV SDisDays,i,p,v,d = xSDiss,i,p,v

d=1

(48)

0

∀s, i, p(: SIPs,i,p = 1 and SP kT ys,i,p = ‘Y ), v(: V ST ys,i,p = V T yv ) Equation 49 captures the daily distribution capacity balance and the variable xV Qr,p,v,d indicates the distribution capacity of vehicle v available at day d in a region r for transferring liquid product p. It is assumed that all the deliveries start at the starting of the day, and as the vehicle comes back to the source completing all the intended delivery that newly

34

ACS Paragon Plus Environment

Page 35 of 63 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

available distribution capacity are added to the distribution capacity pool of the next day.

xV Qr,p,v,d = xV Qr,p,v,d−1 − −

I P

i,c,p =1) c I(:ICP P P

xV CDisDayc,i,p,v,d

c=1 i=1(:IRi,r =1) I(:IIPi,i0 ,p =1)

P

xV IDisDayi,i0 ,p,v,d −

i=1(:IRi,r =1) i0 =1(:IRi0 ,r =1)

+

i,c,p =1) c I(:ICP P P

s,i,p =1) s I(:SIP P P

xV SDisDays,i,p,v,d

s=1 i=1(:IRi,r =1)

C d−1(:DEd−1 >(DSdd +T rTi,c,v,p ))

P

xV CDisDayc,i,p,v,d

c=1 i=1(:IRi,r =1 dd=1(:DSd−1 (DSdd +T rTs,i,v,p ))

P

xV SDisDays,i,p,v,d

s=1 i=1(:IRi,r =1 dd=1(:DSd−1 (DSdd +T rTi,i0 ,v,p ))

+

P

P

P

i=1(:IRi,r =1)

i0 =1

I dd=1(:DSd−1