An optimization model for site-wide scheduling of coupled production plants with an application to the ammonia network of a petrochemical site

This contribution presents the modeling and optimization of the operation of production plants that are coupled via distribution networks and applies it to a part of the petrochemical production site of INEOS in Köln in Germany. The problem is formulated as a mixed-integer linear problem and solved to generate an optimal monthly plan for a set of plants, tanks, and loading/unloading facilities, while respecting various constraints arising from technical limitations, physical couplings between the plants, production targets, and the schedule for import and export across the company borders via ships and trains. The optimization problem takes into account varying energy prices, the influence of the ambient temperature on the processes, and the inventory management for different types of storages. We solve the optimization problem for the particular case and compare the results for a 1 month scenario to recorded data and show that a significant energy saving potential exists. We discuss the current limitations and outline potential improvements in the context of the application of the optimization model to optimal site planning that leads to an improved coordination of the production in the process industries.


Introduction
Companies in the process industries operate under tough competition on the markets and under pressure to improve their energy and resource efficiency and to become "greener" (Baos et al. 2011). The value chain of chemicals with its logistics and complex production processes is a large and tightly coupled system of systems (Engell and Sonntag 2016). The individual production steps ranging from the processing of raw materials, the production of intermediates to the formulation of tailored products and fine chemicals are coupled by flows of material and energy, usually with recycle streams. Often, the processing steps are geographically close and are performed at production complexes or sites such as in the petrochemical industry. The plants are physically coupled by shared resource networks through which for instance steam at different pressure levels or intermediate products are exchanged. These production sites cannot be regarded as autonomous systems that are decoupled from the environment, as complex interactions exist with the logistic network in which these systems are embedded (see Fig. 1). Often, the production is driven by the availability of resources or by customer demands. In addition to the complexity of the logistics, varying material and energy prices give rise to a dynamically changing pricing landscape in which a site has to operate (Mitsos et al. 2018). These circumstances offer a significant potential for joint optimization that results in a well orchestrated overall system in which the available resources are used efficiently and in which the overall system has an increased responsiveness towards the external factors and constraints (Lund et al. 2015).
The scope of optimization in the process industries is often limited to specific issues and domains, but recently there has been an increased interest and research activity in expanding the scope (Castro et al. 2018). Figure 1 illustrates different domains which the scope of optimization can potentially embrace. Looking at the system from a bottom-up perspective leads to the natural choice of starting with the optimization on the plant level. If the focus is widened to the scope of a complete production site and vertically across the different layers of the organization Fig. 1 Conceptual illustration of the scope of the paper (Harjunkoski et al. 2014), the optimal shared resource allocation within chemical production processes by means of distributed market-based coordination schemes is a way to drive the operation of the individual processing plants to a site-wide optimum (Wenzel et al. 2016). Looking at the upstream connections one can optimize the logistic (pipeline) network with the assumption that the processing plants are consumers with fixed inflow rates as Liu et al. (2016) propose. If one is concerned with the on-site storage and buffer capacities and their management coupled with the production, an optimization of the inventory can be done as in (Rodriguez et al. 2018). Considering the external influences such as varying availability of energy from the grid leads to an optimization of the electric power procurement of the site, known as demand side management (Zhang and Grossmann 2016;Hadera et al. 2016;Leo and Engell 2018).
Optimizing the system as a whole, even if the modeling depth is limited, is a step further towards smart factories in the petrochemical industry (Li 2016), which are able to adjust their operation to the situation outside of the site such as low or even negative power prices due to the availability of electric power from renewable resources (Paraschiv et al. 2014).
New or improved solvers, more powerful computers, and efficient modeling frameworks enable formulating and solving more and more complex optimization problems and combining the scopes of the different optimization fields shown in Fig. 1. In this contribution we develop a generic MILP optimization model that is suited to generate an optimal schedule for the operation of physically coupled processing plants under the consideration of external influences such as energy prices and logistic constraints. We motivate the development with a case study of a part of the petrochemical production site of INEOS in Köln in Germany, where an optimal schedule for the operation of the plants within the ammonia distribution network is aimed at.
The rest of the paper is organized as follows. First, the scope of the site-wide optimization and shared resource allocation is laid out and possible choices of optimization models are discussed. Afterwards, the ammonia distribution network at INEOS in Köln is presented as a motivating example. Then a generic model for the purpose of site-wide scheduling of connected processing plants is formulated as an MILP. The model is parametrized for the INEOS in Köln ammonia distribution network and the results are compared to historic production data that were recorded. A summary and an outlook on further developments conclude the paper.

MI(N)LP for site-wide optimization and shared resource allocation
Site-wide optimization (Cheung and Hui 2004) or even enterprise-wide optimization (Grossmann 2005;Varma et al. 2007;Wassick 2009;Quaglia et al. 2013) is a field of research and application that is highly relevant to the process industries and although the idea is not new, recent advances in mathematical programing and solution strategies have enabled the application to real world problems, though the remaining challenges are non-trivial (Grossmann 2012).
One of the challenges is the complexity of the underlying physical and chemical principles which govern the behavior of the individual processing units. Reaction kinetics and thermodynamics lead to significant non-linearities and the presence of different pieces of equipment, redundant production lines, or different operating windows lead to discrete decisions. Consequently, one ends up with the dilemma of either formulating the overall optimization problem as a mixed-integer non-linear program (MINLP) for which often only approximate solutions can be found (Belotti et al. 2013), or one decides to set up and solve a mixed-integer linear program that results from a (piece-wise) linear approximation of the nonlinearities (Grossmann 2012).
Applications of MI(N)LPs in the process industries are numerous and diverse, since the benefits that a system-wide optimization offers are significant. In the crude oil sector, refineries can for instance be modeled as large storage units and the distillation processes as constant level consumers (Más and Pinto 2003) and the resulting MILP can be solved by a decomposition strategy to handle its complexity. Production planning for petroleum refineries was done with a multiperiod optimization model by Neiro and Pinto (2005). Liu et al. (2016) propose a MILP model for the optimal operation of different gas oil separation plants that are operated in a joint distribution network for which the routing is optimized. In Zhao et al. (2017), an enterprise-wide optimization is proposed that incorporates the refinery and an ethylene plant in an integrated topology for a joint optimization. Further, a site-wide optimization is also of interest for a combined optimization of utility consumption such as in heat-exchanger networks [see Kermani et al. (2017) and the references therein]. In this contribution, a generic formulation of an optimization model is proposed for the generation of an optimal operation schedule for processing plants that are physically coupled by shared resource networks and storage tanks. These shared resources can be different streams of material and energy. An application of enterprise-wide optimization to shared resource networks is presented by Lv et al. (2018) who optimize an interplant fresh water network. In Khor et al. (2014) an extensive review on the optimization of water networks can be found. Martí et al. (2013) investigate the optimal sharing of oxygen in a reactor network. The optimal management of hydrogen networks in a petrol refinery is described by Sarabia et al. (2012). In Wenzel et al. (2016) price-based coordination is used to optimally allocate steam on different pressure levels at a petrochemical production site and in Wenzel et al. (2017) the work is extended to incorporate the intermediate products ethylene and propylene. The optimization of a utility network combined with the optimization of the hydrogen production is done based upon a MILP formulation by Hwangbo et al. (2017).
Extending the optimization of shared resource networks, the scope of the model in this contribution is chosen according to Fig. 1 to incorporate decisions from different operational levels (Castro et al. 2018). The level of detail and the choice of the model equations are motivated by the application to the ammonia distribution network of INEOS in Köln, which is explained in more detail in Sect. 3. Since the current planning practice at INEOS in Köln is constraint driven rather than dominated by non-linearities, we propose to use a large-scale MILP formulation to optimize the network including the tanks and the operation of the individual plants on the site, which are modeled by linear and affine functions as done similarly, e.g., by Gao et al. (2015).

The ammonia distribution network at INEOS in Köln
In this contribution, the ammonia production plant (P1) as well as three connected processing plants of INEOS in Köln are considered. The topology of the setup is shown in Fig. 2. The ammonia plant receives hydrogen (H2) from the cracker complex at the site and natural gas (NG) as reactant streams and fuel gas (FG) as the main energy carrier from the shared resource networks. The produced ammonia can be sent to two types of storage. One type are the deep-cooled storage tanks, Tc1 and Tc2, to which the ammonia is sent via two compressors. In the case of high demand, the stored ammonia can be taken out of these tanks. It is warmed up by a heater, used in the downstream plants or sold to the market. The second type of tanks are the three spherical warm buffer tanks Tb1a, Tb1b, and Tb3, which are considerably smaller than the deep-cooled storage tanks and which are used to buffer the production as well as the ammonia delivery and shipping facilities that transfer the ammonia across the borders of the site via ships or train vessels. The ammonia that is further processed on site is sent to three major ammonia consuming plants: The nitric acid plant (P3) and the two acrylonitrile (AN) plants (P4a) and (P4b), which differ in their size and in the number of reactors. The  Table 5 nitric acid plant uses nitrogen from the ambient air as the second main reactant, the two acrylonitrile plants receive propylene (C3) from the cracker complex.
The purpose of the optimization of this network is to operate the plants at the overall optimum, which might not be identical to the optimum for each individual plant, as external influence factors (ambient temperature, prices for electric power and gas) have to be taken into account as well as the operating points of the connected plants. In addition, the logistics within and outside of the ammonia distribution grid play a major role, as they impose important constraints on the operation of the plants. The presence of numerous discrete decisions renders this problem a challenge both from the modeling perspective as well as from the optimization perspective. The integration of external influences into the site-wide ammonia distribution optimization offers a large potential to adjust the operation of the complex in order to achieve an overall energy and resource efficient operation with an increased responsiveness to external influences.

Formulation of a generic MILP for site-wide scheduling
In the following, a generic model that is suited to determine an optimal schedule for the operation of processing plants at a production site where the plants are coupled by several networks is formulated. The model considers the scope illustrated in Fig. 1 and is motivated by the case study shown in Fig. 2.

Nomenclature
The symbols and the naming conventions used throughout this contribution are listed in Tables 1, 2, 3, and 4. The network connects several systems s ∈ S . The systems are tanks t ∈ T , production plants p ∈ P and nodes n ∈ N , which are connected by pipes (r, s, s � ) ∈ . The mass flows through the pipes in time interval j are denoted by ṁ (j) r,s,s � with resource r ∈ R going from system s to s ′ . In every pipe only one resource is transported. A resource can be a pure raw material, a mixture with a fixed composition, or a carrier of energy as, e.g., steam at a certain pressure level. The amount of generation or consumption of electric power in system s in time interval j is denoted by P (j) el.,s .

Nodes
The streams are joined or split in the nodes N. The mass balance of the incoming streams in,n and the outgoing streams out,n has to hold for every time interval j ∈ J and for all nodes n ∈ N . It can be expressed by Many streams on site are bounded by technical limitations. All mass streams are modeled as positive streams, i.e., the flow direction in every pipe is fixed. Thus, the following constraints are formulated:   The stored mass of resource r in system s in interval j The mass flow of resource r from system s to system s ′ in interval j m (j) prod,p The production of the main product of plant p in interval j, equal to the associated flow ṁ s,m,m � for mode transitions are used. In (4) the systems are restricted to be in one mode at a time and in (5) the transitions between the modes are formulated.
Forbidden transitions are prohibited by forcing the respective binaries to zero:  Often, it is required that the plant stays in the new mode for a fixed time K fix,j s,m,m ′ before the next transition, e.g., because shutting down a plant or starting it up cannot be done instantaneously. These constraints are formulated as Note that for a minimum or a maximum stay time instead of a fixed stay time constraint (7) can be reformulated by changing the relation (see Mitra et al. 2013).
For the different transitions between the different modes, the values of K fix,j s,m,m ′ are assumed to be either fixed values, or they are dependent on the ambient temperature. For all transitions with fixed times f the estimated required time is computed as where ̂a is the estimated ambient temperature in time interval j. Obviously, the accuracy of the estimation depends on the length of the horizon, since there is no accurate weather forecast for more than a few days. One possibility is to use seasonal average temperature patterns for night and day times to optimize the schedule. Based on experience, the minimum and maximum times for starting up and shutting down are known for the specific units. At low ambient temperature a the variable part of the time v is and at high ambient temperature a the variable is v = . The variable estimated time v (̂a) is then computed by simple interpolation for temperatures ̂( j) a ∈ [̂,̂] . Outside of this temperature range the required time v (̂( j) a ) is bounded according to

Plant models
Since the derived models are intended to be used for medium to long term scheduling, the plants are modeled as stationary input output mappings and the dynamics of ramping up or ramping down the production rate are modeled by the above constraints. Plants can either be operated as a whole or they can have multiple parallel units, such as reactors with different operating modes. Details will be provided in the following subsections.

General plant model
For the plant models the production rate ṁ (j) prod.,p of plant p at time j is a continuous model input. Further, the ambient temperature (j) a is considered as a non-influenceable input parameter. It is assumed that within the operating range, the remaining streams of material and energy can be determined from ṁ (j) prod.,p by affine equations. The following equations describe a plant that can only be operated as a whole: where the binary decision variable y (j) p,on in the first terms ensures that the plant does not consume or produce resources if it is not on. With this approach, a static consumption or production of resources of a plant can be accounted for, even if no product is being produced, which is relevant, e.g., if cooling water is continuously consumed independent of the load of the plant.

Plants with multiple parallel units
In Fig. 3 the modeling of a plant p with multiple parallel units is illustrated. In a plant with multiple parallel units, the planned overall product stream ṁ (j) prod.,p is split in an internal node into single streams ṁ (j) prod.,p,e u that enter the individual units e u ∈ E U,p . The mass balance is ensured by Similar to (10f.), the model of a single unit is a set of affine functions of the individual product stream with the ambient temperature as external influence prod.,p . The ambient temperature a is a non-influenceable input parameter. Incoming and outgoing streams are split and joint in internal nodes (circles, see (12ff.)) As illustrated in Fig. 3, in a plant with multiple parallel units, there is a load independent block (li) that accounts for the fraction of consumption or production of the plant that is independent of the amount of product ṁ The overall electric power consumption of plant p is formulated accordingly

Operating modes of plants and units
In general, all plants and units can be in the modes on and off (see Sect. 4.3). Further modes can be defined, if needed. Typical additional modes are the modes shutting down and starting up (Mitra et al. 2013). These modes are required if the time required for a shutdown or a start up are modeled with minimum or fixed stay times as in (7). If the modes shutting down and starting up are used, the direct transitions from on to off and vice versa are modeled as forbidden transitions.
For plants with multiple units, relations between the plant mode and the modes of the units have to be established. The constraint (20) guarantees that a plant can only be in mode on if at least one unit is in mode on, i.e., if all units are switched off, than the plant is switched off, too. Constraint (21) guarantees that the plant must be in mode on as soon as one of the units is in mode on.

Constraints on the load changes
To account for the response times of the different plants on site, ramping constraints are introduced which prevent drastic set point changes across the operational window of a unit or plant. Although the models are steady state models and have thus no inherent dynamics, the incorporation of the ramping constraints adds dynamics to the model.
where the latter part ensures that when a mode change from or to on is performed, the change of the mass flow from or to the maximum or minimum value is not restricted. Note that with this modeling approach a plant can switch from the modes off or starting up to the mode of operating at full load (on) instantaneously.

Tanks
The tanks are described by mass balance equations. All incoming streams increase the level, all outgoing streams decrease it. The mass is balanced over the length of the time interval t Each of the tanks has lower and upper bounds on its level. The bounds result from technical limitations or from management decisions, e.g., to keep a certain level of inventory for the case of unforeseen shutdowns: The basic modes for tanks that are considered are fill for filling, disc for discharging and idle, if neither is done. A mode-dependent reformulation of (2f.) for the inflows and outflows of the tanks is These constraints ensure that a tank can only be filled if it is in the filling mode and it can only be discharged if it is in the discharging mode. Depending on the type of tank, uniqueness of the modes can be enforced [see (4)]. This is necessary if a tank can only be filled or discharged at a time.

Logistics of external sinks and sources
Usually, there are several incoming and outgoing streams of material and energy across the boundaries of the system under consideration such as the streams that are depicted in Fig. 2. In our model, there are two types of incoming and outgoing streams. For one type it is assumed that the streams are available at sufficient quantity at all times, i.e., there is no shortage in supply within the predefined operating ranges of the plants. They are provided by adjacent processes and have lower and upper bounds on the quantity per hour that can be acquired. For the second type, discrete imports and exports of resources across the boundaries of the system are considered. These imports and exports can be realized by, e.g., trains, ships, or trucks. Technical limitations for the import and export can either be realized as bounds on the incoming and outgoing streams or, in addition, handled using discrete modes of the arrival nodes. One example is the restriction to only load or unload one type of delivery. Then the receiving node has a unique operating mode from of a set of modes, e.g., { ship, train, truck} [see (4)] and the bounds are modeled according to (26f.).

External influence factors
The prices for electric power and natural gas from the grid as well as the ambient temperature vary during the day and month. In Fig. 4, the electric power prices and the gas prices in an exemplary month are plotted in €/MWh. The electric power price besides the ordinary day-night and weekdays-weekend pattern shows some hours in which the price is very far away from the average. These situations can be, for instance, holidays or weekend days where there is lots of wind, or sunshine, or both, and few consumers operate at their full capacity. During these hours it can be beneficial for the site to operate equipment with a high energy demand such as the large compressors. Although for the period shown in Fig. 4 the gas price fluctuates only little, changes in the gas price happen frequently and can be significant (Hulshof et al. 2016). Then it is beneficial to adjust the operation of major natural gas consumers on the site to save raw material costs during demand peaks or to use it in the case of excess supply, when the prices are low (Hadera et al. 2016). A broader review on demand response in the context of smart grids is provided by Siano (2014).
In Fig. 5 an exemplary ambient temperature plot of Cologne, Germany during a spring month is shown. It can be seen that the temperature differs significantly during the day as well as during the month. The ambient temperature has an influence on the plant and unit model equations (10ff.) as well as on the fixed and variable transition times of the different operating modes that are modeled by (8).

Optimization of the ammonia network of INEOS in Köln
In the following, the scope of the model is defined and the model is parametrized for the ammonia network of INEOS in Köln. The objective function and its degrees of freedom are stated. Afterwards, the model is solved to generate a schedule for one month and its results are compared with historic data of INEOS in Köln to investigate the quality of the solution of the model.

Definition of the scope of the optimization
The aim of the optimization is to generate an optimal production schedule for the processing plants in the ammonia distribution network at INEOS in Köln, which is illustrated in Fig. 2 and described in Sect. 3. The scope of the optimization includes all areas of optimization that are shown in Fig. 1.

Parametrization the optimization model
The model is parametrized by the topology given in Fig. 2 and by the technical, operational, and managerial constraints provided by INEOS in Köln, which will be discussed in the following sections.

The processing plants
The processing plants in the network are listed in Table 5. The ammonia plant is the only ammonia producer on site. The other three plants consume ammonia. The plants all have the discrete modes on or off. Both acrylonitrile plants have individual reactors that additionally have the transition states shutting down and starting up. The transition states have fixed stay times that depend on the ambient temperature as modeled by (7).

The buffer tanks and storage tanks
As described in Sect. 3, there are two different types of ammonia storages. One type are the two deep-cooled storage tanks Tc1 and Tc2. The second are the spherical buffer tanks. All buffer tanks and storage tanks are listed in Table 5. In the following, details for the two storage types are given.
Deep-cooled storage tanks for ammonia T c ⊂ T cannot be filled and discharged at the same time (4). There are no restrictions on changing between the three modes and thus switching can be done after every time step t . The tanks are only used for long term storage to keep a certain level of inventory for the situation of an unforeseen shutdown or to build up an amount of ammonia for sale. Before the ammonia is stored in these tanks, it has to be cooled down by the compressors C1 and C2. When discharging the tank, the ammonia has to be heated up on the way to the warm buffer tanks. This operation requires energy that could be reduced by employing an optimal schedule of the overall network. During one time interval ammonia is only sent to one of the two tanks or only discharged from one of the two tanks. This constraint is based on a heuristic of the plant personnel to minimize the pressure loss due to operating two pipes at the same time. It can be formulated as Ammonia buffer tanks T b ⊂ T can be filled and discharged at all times. However, they have outlets to different nodes. They can discharge towards either of following processes (discp), ships (discs) or train vessels (disct) or stay idle (idle). However, there is no restriction to only one mode being active in each interval, e.g., it is possible to send ammonia to processes and a train vessel simultaneously. As in (27), the mass flows are set to zero if the corresponding discharge mode is inactive. There are no forbidden transitions between the available modes. In this contribution the imported amounts of ammonia via ships and trains are implemented as hard constraints. Thus, the discharging modes for trains and ships are fixed throughout the horizon (see Sect. 5.5). In theory, it is possible to import and export to and from ships and trains via all three ammonia buffer tanks. However, the current practice at INEOS in Köln is to only use Tb3 for import and export. Export suppression from Tb1a and Tb1b is enforced by The restriction of importing ammonia only via the buffer tank Tb3 is formulated as where ext denotes external sources.

The optimization problem
The optimization problem is stated in a condensed form as the following MILP where f ∶ ℝ n m ×ṅm×n P el. ×n y ×n z → ℝ is the objective function and m, ṁ , P el. , y, and z are the decision variables with y and z as binaries. The optimization is formulated for 31 days, with a resolution of t = 1h , which results in a horizon J = {1, 2, … , 744} . The number of variables 1 for the tank levels is n m = 7450 . The number of variables for the pipes is ṅm = 93,744 , the number of variables for the power n P el. = 2232 , the number of the binary mode variables n y = 40176 , and the number of the binary transition variables n z = 138,384 . With this, the total number of variables is 281,986. The choice of the resolution of t = 1 h , results in a large-scale problem which could in principle be avoided by choosing a more coarse resolution. However, the fine resolution is needed in this case study to account for the fixed stay times. In addition, we observed infeasibilities for the buffer tanks, ships, and trains which are caused, e.g., if a discharging mode of a tank is active for too long, which causes the tank to run empty.

Objective function
The objective function is formulated over the scheduling horizon ∀j ∈ J . It contains the following contributions The economic performance of the site over the horizon in terms of profit and additional costs that accompany certain transitions are accounted for. The terms contributing to the objective are discussed in the following in detail.

System profit
For all plants and tanks (systems) the profit is calculated as the difference of the exported and the imported mass weighted by an (internal) price of the specific resource. In addition, the cost for electric power is considered where (j) r,s,s � is the price of resource r for a specific stream and (j) is the time-variant price of electric power. It is important to note that the prices (j) r,s,s � are also time-variant for natural gas and they can be in principle time-variant for other streams, too. When summing up the profit terms for all tanks and plants the internal sales cancel out and only accumulation in tanks and transfers across the site boundary appear in the total balance.

Transition costs
Shutting down and starting up a reactor or plant does not only require a specific amount of time (see Sect. 4.3), but is also associated with costs. The costs for these procedures are taken from past experience and are included in the cost function. For instance shutting down a reactor requires a predefined cleaning procedure that has a fixed associated cost. Accounting for the transition costs is done by the binary transition variables and individual prices s,m,m ′ for all transitions that lead to additional costs We assigned marginally different values s,m,m ′ for identical equipment such as the reactors to break ties. The following relation is an example for the reactors in plant P4a: This formulation improves the solution time, because the time that the optimizer requires to evaluate symmetric solutions (Margot 2010) is reduced. The optimizer chooses the cheapest transition among the transitions of identical equipment and if is small, the structure of the overall solution is not changed. In this case study no extra penalty on the size of the load change is considered. However, such a penalization could be added in a straightforward fashion by a linear reformulation of the absolute change of the load prod.,p |ṁ (j) prod.,p | , where prod.,p is a parameter that penalizes the size of the load change.

Formulation of the optimization scenario
To determine the quality of the solution of the site-wide optimization model, a problem instance of an example scenario provided by INEOS in Köln is formulated and solved. Past recorded production data from one sample month 2 has been taken as a reference. The optimization problem takes the monthly production target and initial conditions for all modes and tanks levels into account. The following sections provide a detailed description of the scenario.

Background
The scheduling horizon is one month in Spring (31 days), where it is assumed that in the subsequent month a plant shut down of the ammonia plant P1 happens. Consequently, all plants have to operate such that for the next month a sufficiently large amount of ammonia is stored on site to ensure an uninterrupted operation of the downstream processes. Given a planned schedule of purchases and sales of ammonia, the tank levels and the operating levels of the plants are optimized in order to reach the defined target amount of ammonia that is stored on site at the end of the month while meeting the predefined monthly production capacities for the products.

Logistic constraints
For the optimization, the schedule of logistics w.r.t. the ships and trains is given in Fig. 6. During the investigated month there are five ships that need to be unloaded. The Fig. 6 Purchases and sales via ships and trains during one month. The loading and unloading is assumed to be done with maximal flow rates. For the sake of illustration, the number of symbols indicates the time that is required for loading or unloading of ships and trains. One train vessel symbol refers to ≈ 5 h of processing time. To unload a ship requires ≈ 12 h of processing time schedule for the trains is illustrated by different numbers of train vessels, where one train vessel requires ≈ 5 h of loading or unloading time. For the processing of these imports and exports the mass flow rates have been fixed, since it is assumed that the processing is done as fast as possible at fixed flow rates.
where the reduced load ũ b (j) NH 3 ,s,s � is necessary to process amounts that cannot be processed during a full number of time intervals n × t, k ∈ ℕ , i.e., if the unloading would require 4.5 h at full load, then during the first four hours the bound is set to full load ub NH 3 ,s,s � . This expresses the same situation as if one would process 4.5 h at full load.

Energy prices and ambient temperature
The energy prices for power and natural gas as well as the ambient temperature vary over time during the scheduling horizon. The time varying prices influence the objective function (33) and need to be estimated for the length of the scheduling horizon J. An extensive review of electric power price forecasting can be found in Weron (2014) and more general for forecasting methods in energy planning models in Debnath and Mourshed (2018). Since in this contribution we simulate the optimization model against a set of historic data of INEOS in Köln, the optimizer knows the prices in advance. In practice, of course only predictions which are subject to uncertainty are given. The prices for power and gas for the considered scenario are shown in Fig. 4 for the scheduling horizon in €/MWh. The ambient temperature varies on different time scales. There are daily and seasonal variations. In the model, the ambient temperature influences the model equations of the resource consumption and production of the plants (10ff.) and it influences the fixed transitions times for shutting down and starting up the reactors listed in Table 5 (see 8). As for the energy prices, in this contribution the ambient temperature is known to the optimizer for the complete scheduling horizon.

Initial conditions and target amounts
For all the systems listed in Table 5, the modes at the beginning of the horizon j = 1 are fixed to the modes that were recorded by INEOS in Köln. Also, the initial levels of all tanks within the network were fixed.
NH 3 ,s,s � if import or export is performed at full load ub (j) NH 3 ,s,s � if import or export is performed at reduced load 0 otherwise ∀j ∈ J, (NH 3 , s, s � ) ∈ in,site ∪ out,site , As explained in Sect. 5.5.1, at the end of the horizon a specific minimum target amount of ammonia that is stored on site had to be realized. Hence, the following constraint is introduced that sums up the amount of ammonia stored in all tanks at is a given value taken from the data of INEOS in Köln.
The current practice for the site-wide scheduling at INEOS in Köln is to fix monthly production targets for the processing plants. These targets can be translated into ammonia consumption targets of the plants that produce subsequent products from ammonia. Thus, the following constraints are added to the problem formulation for the plants P3, P4a, and P4b where ṁ (j) NH 3 ,t,p denotes the amount of ammonia that every plant p (except for the ammonia plant) receives from the buffer tank t during the scheduling horizon (see Fig. 2). The fixed consumption target m NH 3 ,p,⊚ is taken from the actual production capacities in the scenario.

Implementation
The optimization problem was implemented with the help of the mathematical modeling package JuMP (Dunning et al. 2017), which is a package written in Julia (Bezanson et al. 2017). This implementation is independent of the chosen solver, since  Programming Kit) 2018) were tested with different settings for the relative MIP gap (see Table 6).

Results and discussion
This section contains the results for the described scenario. The optimized schedule and its production rates and tank levels are compared to the historic daily operation at INEOS in Köln for that month. For the discussion of the optimized schedule we use the results with the smallest MIP gap, which were obtained within 60 min by the solver Gurobi (see Table 6). Figure 7 shows the optimized hourly production rates of the plants with respect to the ammonia usage, i. e, the ammonia plant P1 produces ammonia and the others consume it. The dashed lines represent the recorded operation at INEOS in Köln. The optimized schedule for the ammonia plant assigns full load towards the end of the month, while at the beginning the load is set to the lower bound, which is about 80% of the NH 3 production capacity. In between, the optimized production trajectory of the plant is ramped up and down twice between the lower and upper bound without a complete shut down of the plant. The ramping is determined by the constraints (22f.). The reasons for the particular pattern can be explained by the pattern of the gas price (see Fig. 4), which is slightly lower at the end of the month. In addition, the ammonia plant reduces its load to not exceed the amount that is taken from the proceeding processes plus what is required at the end of the month as defined in (37) and (38). The ammonia production does not exceed the minimum required amount of ammonia on site, because additional ammonia has to be stored in the deep-cooled storage tanks Tc1 and Tc2, which requires an additional consumption of electric power without the possibility of sales within the horizon.

Operation of the plants and their reactors.
The schedule that is computed for the nitric acid plant P3 has the largest number of set point changes among the plants. The plant is not shut down completely, since for this action a high transition cost is assigned by (34). It is more efficient to reduce its load, which is done by ramping up and down between the upper and the lower bounds of the ammonia takeup capacity. It is visible that there exists a coupling between the operation of plants P1 and P3. At the beginning of the horizon the load of P1 is reduced and the load of P3 is reduced, too. However, the coupling is not straightforward, since the plants are connected via the ammonia storages. The intermediate reduction of the load of P3 is either caused by a shortage in the buffer tanks or by a propagation of the load reduction of P1 during an earlier time interval. Similar to the operation of P1, P3 does not operate at maximum load throughout the horizon in order to not exceed the monthly target. The required operational load is shifted according to the needs of the site.
In the optimized schedule for the AN plant P4a all the reactors are operated either at full load, or switched off completely, which is done in the last quarter of the horizon (see Figs. 7,12). For this plant it is more beneficial to shut down individual equipment and pay the one time transition cost (34) instead of reducing the load of the reactors. A similar situation can be observed for the other AN plant P4b, which switches off one of the four reactors in the first third of the horizon (see Fig. 7, 12). For these reactors it is more economical to operate them at their maximum load or to switch them off, which is related to large constant terms in the affine model equations (13f.). The assignment of the timings for the shut downs is done either because of the prices for gas and electric power or according to the availability of ammonia in the buffer tanks, which is affected by the purchase and sale schedule and which is coupled to the operation of the remaining plants.
When the optimized schedule is compared to the recorded operation one can see that the optimizer makes use of the possibility to switch off reactors or even complete plants. The recorded operation at INEOS in Köln changes the load of the plants slowly rather than performing sudden changes, which is typical for the operation of large continuous processes. In contrast, the optimizer evaluates the tradeoff between shutting down and operating at a less efficient operating point such as the lower bound.  Optimized tank levels. Figures 8, 9, 10, and 11 give insight into the optimization results for the tank levels. Figure 8 shows the optimized tank levels of all individual ammonia tanks within the network. It is visible that the optimizer strives to fill all tanks towards the end of the month, which is in line with the required value defined in the scenario (37). The difference between the two different tank types, namely the deep-cooled storage tanks Tc1 and Tc2 (denoted as TcX) and the buffer tanks Tb1a, Tb1a, and Tb3 (denoted as TbX) becomes visible in Fig. 9 where the levels of the grouped tank types are compared to the recorded data. The integrated amount of ammonia on site matches exactly the amount that was recorded. However, the distribution is different. Although, the optimizer and INEOS in Köln use the tanks TbX to buffer the imported and exported amount (wavy pattern of TbX), the final amount in the buffer tanks is larger compared to the recorded data and the amount in the deep-cooled storage tanks is smaller. This results from minimization of the costs on site, where the optimizer strives at spending as little money as possible for storing ammonia. Sending ammonia towards the deep-cooled storage tanks via the compressors requires liquefaction, which adds extra costs and reduces the resource efficiency on site. Figure 10 gives further details on the filling of the individual deep-cooled storage tanks Tc1 and Tc2. Starting from the initial values, both tank levels increase. The optimizer chooses to first fill Tc1 at a higher rate, then Tc2 is preferred from 400 h on. In contrast to the recorded data, the optimizer decided to alternate between the filling of the two tanks. This is possible because switching between the two tanks can be done after every time interval t without any further cost. Overall it can be seen that especially Tc1 has still some capacity left at the end of the planning horizon in the optimized schedule, which is caused by the minimization of the storage cost on site and thus by a minimization of the amount that is stored in the deepcooled tanks. Figure 11 shows the trajectory of the buffer tank Tb3. This buffer tank processes the flows of ammonia that enter or leave the site via ships and trains (see Sect. 5.5.2). It is clearly visible that the pattern of the ship arrivals is reflected in both the recorded data and the optimized schedule (cf. Fig. 6). The optimized schedule generally has a higher filling level, which is especially of importance at the end of the horizon, where more ammonia is stored in comparison with the recorded data. The optimizer uses this free capacity to store ammonia on site without the necessity to cool it down. Further differences can be found in the emptying procedure. The tank needs to be nearly empty to receive the load of a ship when it arrives. While in the recorded operation the tank is emptied earlier than required, which might be due to experience and implicit uncertainty handling of the site management, the optimizer empties the tanks just before the ship arrives, since there is no safety margin included in the formulation.
Optimized operating modes. Figure 12 illustrates the operating modes of the major equipment and plants on site. Most of the plants stay switched on throughout the horizon, which is plausible for the operation of large petrochemical processes. Reactors are switched off only in the plants P4a and P4b. Plant P4a is completely shut down (both of its reactors are switched off) and in plant P4b one reactor is switched off. It is also visible that the shut down and start up times differ from reactor to reactor and in addition they are temperature dependent [cf. the shut down time of the reactors P4aReac1 and P4bReac1, see (8)].
It is interesting to analyze the switching of the compressors C1 and C2, which is predominantly determined by the need to liquefy ammonia to reach the target capacity at the end of the month in the deep-cooled storage tanks. In addition to these needs, the optimizer considers the price of electric power to run the power consuming devices during times of low prices of electric power as illustrated for the compressors in Fig. 13. In the figure, a selection of regimes of low and high prices for electric power are highlighted. It can be seen that the optimizer chooses to switch the compressors on during times of low prices and decides to turn them off if the price of electric power is elevated.
Energy savings for deep-cooled ammonia storage. There are many aspects that can be considered to evaluate the quality of the optimized schedule. One aspect is to look at the electric power required for the storage of deep-cooled ammonia. There are two factors that lead to savings. First, the alignment of the liquefaction of ammonia with the power price leads to savings. Second, the minimization of the amount Fig. 13 Operating modes of the two compressors plotted jointly with the prices for electric power. The vertical stripes illustrate exemplary regimes of low and high prices for electric power Fig. 14 Relative savings in power because of less liquefaction of ammonia in the deep-cooled storage tanks and alignment of the compressor operation to the power prices. 100% corresponds to the costs that occurred on site for the given price of electric power and the processed amount of ammonia of liquefied ammonia reduces the power consumption. As an example, the relative power expenses integrated over the horizon are shown in Fig. 14. At the beginning, both the recorded and the optimized trajectories start at zero expenses. Since the power price is negative at the beginning (see Fig. 4), the relative expenses become negative. The data is scaled such that at the end of the horizon the recorded operation reaches 100% of the expenses. It can be seen that the optimized schedule leads to less expenses than the recorded data. At the end of the horizon, the optimized schedule saved ≈ 25% of the costs for the liquefaction, which is caused by less liquefied ammonia and the adjustment of the compressors to the power prices.
Solution times. The solution times required to solve the optimization model significantly depend on the scenario and how many logistic constraints and operating modes are prefixed. In addition, we tested the performance of different solvers, two commercial ones and two open source solvers. The results are shown in Table 6 for varying tolerances of the MIP gap of the solvers. Both open source solvers were not able to find a solution to the problem instance within the allowed time of 60 min. In terms of the objective values, CPLEX and Gurobi perform similarly if the allowed time is 60 min and the tolerance for the relative MIP gap is set to 10 −4 . For the most coarse termination tolerance of 10 −2 , Gurobi finds a solution within the tolerance much faster than CPLEX. For a termination tolerance of 10 −3 both solvers require comparable solutions times with similar relative MIP gaps.

Conclusion
From the generic model proposed in Sect. 4 it is possible to generate an optimal schedule for the operation of processing plants coupled by resource networks. The results for the provided scenario of INEOS in Köln could be obtained within a reasonable time if the termination tolerance is not chosen too tight. The results are plausible and they are in good agreement with the recorded data. The optimizer reduces the production cost by shifting the operation periods and loads of the processes and the operation of energy consuming equipment. In the perspective of demand side management the benefits are limited because of the inflexibility of the processes, which is not uncommon in the petrochemical industry. The large production plants are slow in changing their operating point and are directly coupled to the operation of the adjacent processes. Since between the ammonia plant and the compressors there is no intermediate storage, the operation of the compressors is strongly linked to the operation of the plant, which is limited by ramping constraints. In addition, the overall system is heavily driven by logistic constraints, which significantly reduces the degrees of freedom for an optimizer to improve the operation. The logistic constraints of the arrivals of ships and trains are fixed in this contribution. However, with the generic modeling approach presented in Sect. 4 it is possible to vary the arrival times of the ships and trains as well.
One issue that we noted is that the contributions of the individual terms in the objective function vary by orders of magnitude. The predominant terms are the production loads of the plants, since these terms determine the profit of the site and are the main driving force. The "minor" decisions such as when to shut down one reactor do not have a large impact on the objective. Hence, the solution for these "minor" decision can vary depending of the chosen tolerance and solver. If one is not satisfied with the pattern of the solution and the respective setting of the operation modes, one can implement further constraints that are inspired by daily operation such as minimum stay time constraints. This will avoid that equipment like the compressors switch their modes too frequently.
Towards a practical implementation, the efficiency of the optimization model can be further improved if more discrete variables are fixed by experience or external references. In this situation, however, the optimization model might represent more the reality of the current operation than revealing the potential improvements of the overall system, i.e., there is the risk of modeling the site-wide operation with all the constraints that are in practice regarded as constraints, but which potentially prevent the system to operate at an even better operating point.
An important point that needs to be addressed is the prediction of the prices of gas and electric power, as well as of the ambient temperature. If it is intended to take these into account, then precise predictions are necessary to adjust the operations in time. The fact that the optimizer in this contribution is aware of the exogenous influences over the complete horizon has to be kept in mind when interpreting and evaluating the results and potential savings.

Summary and future work
We first presented a generic MILP model for the optimization of the operations of coupled production plants that can have multiple units. The optimization model takes into account the storage tanks and the logistics of the site as well as demand side management. The model was parametrized for the petrochemical production site of INEOS in Köln, where four processing plants are physically coupled by a joint distribution network. For a scenario based on historic data, an optimization problem was formulated to generate a one month schedule for the operation. The solution of the optimization problem is in good agreement with the recorded operation of INEOS in Köln, since the relevant constraints for the operation of the network are captured in the model. We showed that there is a significant saving potential with respect to the power consumption for the deep-cooled storage tanks on site. The optimization model was solved by the two commercial solvers CPLEX and Gurobi, while the two open source alternatives failed. Gurobi found a solution with a slightly smaller MIP gap within an allowed time of one hour.
Future work includes the incorporation of educated guesses or estimations for the energy prices as well as the ambient temperature, since the required planning horizon is supposed to be longer than a reasonable prediction horizon of the energy market or the weather. In addition, to cope with the inherent uncertainty of train and ship schedules as well as unforeseen failures in the production, a moving horizon formulation and daily rescheduling is proposed. This can be extended by stochastic formulations (Sand and Engell 2004;Cui and Engell 2010). Another interesting aspect to study is the extension of the described network by more plants and possibly by connecting it to more shared resource networks such as steam or other base chemicals shared on site. Further, it can be interesting to apply distributed optimization algorithms such as market-based coordination in order to perform the optimization of coupled production plants that belong to different business units or companies that operate at large chemical sites or in industrial clusters as done by Wenzel et al. (2016) for coupled continuous plants.