**Commemorative Issue to Celebrate the Life and Work of Prof. Roger W.H. Sargent**

Special Issue Editors

**Rafiqul Gani Ian Cameron**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Special Issue Editors* Rafiqul Gani Ian Cameron PSE for SPEED Company Ltd. 294/65 RK Office Park Thailand

The University of Queensland Australia

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Processes* (ISSN 2227-9717) (available at: https://www.mdpi.com/journal/processes/special issues/ Roger Sargent).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Article Number*, Page Range.

**ISBN 978-3-03928-134-2 (Pbk) ISBN 978-3-03928-135-9 (PDF)**

Cover image courtesy of shutterstock.com.

c 2020 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**



## **About the Special Issue Editors**

**Rafiqul Gani** is the co-founder and CEO of 'PSE for SPEED', a company located in Bangkok, Thailand. He was formerly Professor of Systems Design at The Technical University of Denmark and Head of the Computer Aided Design Center (CAPEC). He has the position of Visiting Professor at Texas A&M University as well as Zhejiang University, China. He received his PhD from the University of London (Imperial College), having also received numerous other academic and professional awards. Rafiqul was also Editor-in-Chief of the Computers in Chemical Engineering journal from 2009–2015 and on the editorial advisory boards of other journals. He served as President of the European Federation of Chemical Engineering for two terms between 2014 and 2018. He continues to pursue the research and application of computer-aided methods and tools for modelling, property estimation, process–product synthesis and design, and process–tools integration.

**Ian Cameron** is Professor Emeritus at The University of Queensland and a Fellow of the Australian Academy of Technology and Engineering. He is also co-founder, director and principal consultant with Daesim Technologies, Brisbane. He worked for 10 years with the CSR Group in industry sectors such as sugar, building materials and petrochemicals, having roles in process and control system design, plant commissioning, operations, risk management and environmental protection. After completing a PhD at Imperial College London in the area of Process Systems Engineering (PSE), he worked for 9 years as a United Nations (UNIDO) engineering consultant in Argentina and Turkey. His research has focused on intelligent systems for risk management and process diagnosis, advanced control and system modelling, and the development and application of virtual reality systems in industry and innovation in engineering education design and delivery. He has received numerous professional awards that include the Australian Prime Minister's Award for university teaching.

## *Editorial* **In Memoriam of Professor Roger W.H. Sargent, the Founder of "Process Systems Engineering"**

**Rafiqul Gani 1,\* and Ian Cameron 2, \***


Received: 20 March 2020; Accepted: 20 March 2020; Published: 30 March 2020

#### **In Memoriam**

**Roger W.H. Sargent**

#### **1926–2018**

In September 2018, the global chemical engineering community lost a true pioneer in the field. This Special Issue is in celebration of Roger Sargent's influence, contributions and legacy as a founder of the field of 'Process Systems Engineering' (PSE). His many years, initially at Air Liquide in Paris, and subsequently at Imperial College London laid a foundation that has been built upon by many who did not personally know him, but have benefited through his insights and sustained vision.

Besides his research and industry focus, he was also influential within education and the engineering profession. He held leadership positions such as head of the Department of Chemical Engineering Imperial College, Dean of the City and Guilds College and President of The Institution of Chemical Engineers (UK) and director of the Centre for Process Systems Engineering.

The papers in this issue represent areas that, in many ways, were pioneered by Roger as he and his many students developed new approaches to the integrated design, mathematical modelling and the optimization of complex industrial processes. These developments were facilitated by computer-based tools and a wide range of advanced numerical methods. A number of his direct PhD students and some of their academic descendants are authors of these contributions.

In remembering his impact, there follows a personal tribute from Arthur Westerberg, one of Roger's PhD students in the early 1960s. Professor Westerberg himself continued ground-breaking contributions in PSE at Carnegie Mellon University (USA).

In the late 1970s, both of us were privileged to complete our PhD degrees under Roger's expert and constructive guidance. We were thankful for his time, patience and encouragement during that period and his ongoing friendship.

He will be missed by family, friends and a host of those within his academic family tree. However, he will not be forgotten.

Rafiqul Gani and Ian Cameron.

#### **Tribute to Roger Sargent**

by

#### Arthur W. Westerberg

I am going to steal part of the following tribute to Roger Sargent from a message I sent to his son, Philip, just after I learned, sadly, that Roger had passed away.

Many have and will summarize Roger's many contributions to the Process Systems Engineering area. This tribute will be more personal. I was Roger's sixth PhD student, from 1961 to 1964, and, of course, he was a lifelong friend.

Having done a two month summer tour organized for college students in 1959 with Shel Thompson, my roommate at the University of Minnesota, and having Neal Amundson as my undergraduate advisor, who knew just about everyone in academia in the UK, pushed me to apply to go to the UK for my PhD. Negotiations with Kenneth Denbigh, then chairman of the teaching side of the department at Imperial College, garnered me an Assistant Lectureship offer to provide support. After spending a year at Princeton to obtain first an MS degree, I stepped off the ship (on which, incidentally, I met my wife Barbara) knowing no one in London in September of 1961. I remember first visiting the main office of the department at Imperial College and meeting Kenneth Denbigh. After first presenting me his topics, he sent me around to the rest of the department. I entered Roger Sargent's office shortly thereafter, and, of course, Roger had a project for someone who had on his resume experience in using the computer, on projects both at Minnesota and at Princeton. No one in those days would likely have had any computer experience. Furthermore, as we now know, Roger was into computers and their use in process design and control, with his several years at Air Liquide in Paris in the 1950s. Mind you the computers anyone would have had then were toys by today's standards, but he could see the potential and wanted to be in the business of developing that area as a research area. As I came with an Assistant Lectureship, an entry level faculty position, he insisted I join the faculty for morning coffee, lunch and coffee, and afternoon tea. Thus I saw him at least three times almost every day and had lots of time to interact with him on my project. Our paper from that work was really two parts: a core dump of his ideas on how computers would contribute to design and control in our field and a presentation of the ideas in my work on how to shoehorn a flowsheeting calculation onto a vastly undersized computer. Many people in industry told me they were in the modeling departments of their companies because they mentioned they had read that paper in their interviews.

I worked for two years with Control Data in San Diego before joining the Chemical Engineering Faculty at the University of Florida in 1967. By then Roger had given a well appreciated talk at the invitation of the AIChE and then published as a paper in CEP. This talk and the resulting paper again presented his views on PSE. If memory serves, it had about 50 references to relevant papers from many different fields. I remember going on a hunt to collect and review all those papers. It got me off to a running start on my career when I joined the faculty in Florida.

A bit of humor: In 1990 Roger received the CAST Division Computing in Chemical Engineering Award and gave his acceptance speech at the November AIChE meeting. As all know who have heard Roger give a presentation, his talks were seldom if ever done within the prescribed time limit. We were all sitting at the tables listening and, by this time, consuming our desserts and drinking coffee. I remember clearly lifting my coffee cup for a sip. When I looked down, the wait staff has taken away the saucer. Our table lit up with laughter. The wait staff were ready to go home.

Over the years I had many occasions to interact with Roger. We were at many meetings together. I visited Imperial College to give the annual Sargent lecture, Roger came to Carnegie Mellon University to give a departmental seminar (we have three former and current faculty members who are his students: the oldest: me, Ignacio Grossmann and Erik Ydstie). In the late 1980s, as we were just starting our Engineering Design Research Center at CMU, Roger led a group at Imperial College to create a related organization aimed at process systems engineering. As director of our center, I reviewed his proposal and others, and clearly this was the project to fund among those I saw. After, I had many occasions to visit his centre. I have to admit, he seemed to tolerate administration (certainly much more eagerly that I did) and held many significant such positions at Imperial College and in the Institute of Chemical Engineers, as well as running this centre. So his legacy is much more than his research.

Finally, he is the top node of a huge academic family tree (http://titan.engr.tamu.edu/Sargent\_tree/ index.php), that is "the tree" for all of us lucky enough to have been on his.

We will all miss him.

**Author Contributions:** R.G. and I.C. wrote the In Memoriam preface. A.W. wrote the personal tribute to Professor Roger Sargent. All authors have read and agreed to the published version of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Towards the Grand Unification of Process Design, Scheduling, and Control—Utopia or Reality?**

#### **Baris Burnak 1,2, Nikolaos A. Diangelakis 1,2 and Efstratios N. Pistikopoulos 1,2,\***


Received: 3 June 2019; Accepted: 17 July 2019; Published: 18 July 2019

**Abstract:** As a founder of the Process Systems Engineering (PSE) discipline, Professor Roger W.H. Sargent had set ambitious goals for a systematic new generation of a process design paradigm based on optimization techniques with the consideration of future uncertainties and operational decisions. In this paper, we present a historical perspective on the milestones in model-based design optimization techniques and the developed tools to solve the resulting complex problems. We examine the progress spanning more than five decades, from the early flexibility analysis and optimal process design under uncertainty to more recent developments on the simultaneous consideration of process design, scheduling, and control. This formidable target towards the grand unification poses unique challenges due to multiple time scales and conflicting objectives. Here, we review the recent progress and propose future research directions.

**Keywords:** process design; scheduling; process control; integration

#### **1. Introduction**

It has been over half a century since Professor Roger W.H. Sargent envisioned a paradigm shift in chemical process design methodologies, from ad hoc engineering judgment for specific problems to fully computerized systematic approaches based on complex mathematical models [1]. He conceived the notions of "explicitly formulating the techniques" and "precisely defining the objectives" for engineering design problems, which used to be considered to be "an activity not worthy of higher minds" because of the lack of scientific and systematic tools and methodologies. With the advent of computers, he further emphasized the opportunity to expand the process design problem to account for foreseeable variations in the plant environment over its life cycle to achieve more reliable and robust operations. "(During the process design phase) Many parameters are left available for adjustment during plant operation, such as flow rates, tank levels, operating pressures, etc., but here also *the design places limits on the range of variation possible*." stated Professor Sargent to underpin the interdependence between the design and the uncertainty of future operational decisions.

The Process Systems Engineering (PSE) community has been accumulating formidable knowledge and know-how on mathematical modeling techniques in the fields of process design and operations, and developed efficient tools to solve these advanced models since Professor Sargent had outlined the future of PSE in his 1967 perspective article [1]. Moreover, it has been long established that the early design problem should be studied simultaneously with the operational time-variant decisions to improve the operability and flexibility of the process under variable internal and external plant conditions, and consequently to achieve more reliable, economically more favorable, and inherently

safer processes. The most recent efforts towards simultaneous consideration of design and operational decisions explore effective methodologies to integrate the short-term process regulatory decisions (process control) and longer-term economical decisions (scheduling) through mixed-integer dynamic optimization (MIDO) formulations. The proposed solution tools and techniques for this class of integrated problems include (i) discretizing the dynamic high-fidelity representation of the process through orthogonal collocation on finite elements followed by solving a mixed-integer nonlinear programming problem [2], (ii) "back-off" approach to ensure constraint satisfaction under some assumed worst-case scenario [3–5], and (iii) multiparametric programming to explicitly represent the operational strategies to derive tractable and equivalent MIDO formulations [6].

In this paper, we present a historical perspective on the development and progress of modern process design techniques that account for the dynamic variability introduced by the process control and scheduling decisions. In retrospect, we observe the evolution of methodologies from fundamental analyses on design and process uncertainty at steady state to dynamic complex models that explicitly encapsulate the scheduling and control decisions, as illustrated in Figure 1, and summarized as follows.


Clearly, it would be inaccurate and redundant trying to reduce the individual research efforts to a single category. The literature is noticeably diverse in this field with numerous different approaches. However, we find it useful to classify into certain schools of thought that are also in alignment with the historical progress of the field. In Section 5, we further seek to pose the pivotal questions on future challenges and opportunities for the seamless integration of the design, scheduling, and control problems based on the cumulative knowledge of the PSE community and the current trends in the academia.

**Figure 1.** A Venn diagram representation of major operability indices and their relationship with process economics. It is interesting to note that the design optimization approaches started from the outermost layer, and with the advance of modeling techniques, they have progressed towards the center for guaranteed operability, which delivers the optimal process economics.

#### **2. Early Efforts in Design Optimization under Uncertainty**

The ongoing collective efforts towards the grand unification of design, scheduling, and control was inaugurated through steady-state design under uncertainty in plant conditions. Takamatsu et al. (1970) [7] estimated the undesirable effects of variations in system parameters, measured process disturbances, and manipulated variables on plant performance by sensitivity analysis on a linearized model. Nishida et al. (1974) [8] adopted the notion of sensitivity analysis to structure a min-max problem for design optimization, presented by Equation (1).

$$\begin{aligned} \min\_{\text{des}} & \max\_{\theta} & \quad \mathbb{C}(\mathbf{x}, \text{des}, \theta) \\ & \text{s.t.} & \quad h(\mathbf{x}, \text{des}, \theta) = 0 \\ & & \quad \text{g}(\mathbf{x}, \text{des}, \theta) \le 0 \\ & & \underline{\theta} \le \theta \le \overline{\theta} \end{aligned} \tag{1}$$

where *x* is the vector of states of the system, *des* is the vector of design variables including the steady-state manipulated variables, *θ* is the vector of parameters that agglomerates the system uncertainties and process disturbances. Equation (1) is one of the first notable attempts to systematically assess the trade-off between minimizing the investment cost and improving the flexibility of the process design. However, this strategy yields conservative solutions since it does not distinguish the time-invariant design variables and time-variant manipulated variables. Grossmann and Sargent (1978) [9] remedied this issue by treating the time-sensitive variables (i.e., manipulated actions and design variables that can be modified in the future) and fixed design variables separately. They further adopted the parametric optimal design problem proposed by Kwak and Haug (1976) [10], and formulated an objective function to minimize the average cost over the expected range of parametric uncertainty, as presented by Equation (2).

$$\begin{aligned} \min\_{\boldsymbol{\mu}, \boldsymbol{\mu}, \boldsymbol{\Theta} \boldsymbol{\varepsilon}} &\quad \mathbb{E} \left\{ \mathbb{C} (\mathbf{x}, \boldsymbol{\mu}, \boldsymbol{\epsilon} \boldsymbol{\epsilon}, \boldsymbol{\theta}) \right\} \\ \text{s.t.} &\quad \max\_{\boldsymbol{\theta} \in \Theta} \quad \mathbb{g}\_i (\mathbf{x}, \boldsymbol{\mu}, \boldsymbol{\epsilon} \boldsymbol{\epsilon}, \boldsymbol{\theta}) \leq 0, \quad i = 1, 2, \dots, t \end{aligned} \tag{2}$$

where the expected cost function is defined the joint probability distribution of the parameter set *θ*. Equation (2) requires solving infinite nonlinear programming (NLP) problems. Grossmann and Sargent (1978) [9] proposed an efficient solution procedure for a special case of Equation (2), where each constraint *gi* is monotonic in *θ*, through discretization of the problem over the parameter space. However, solving the NLP problem at a finite number of *θ* realizations does not ensure the feasibility of the design. This issue was addressed by Halemane and Grossmann (1983) [11] through reformulating an equivalent design feasibility constraint as presented by Equation (3).

$$\max\_{\theta \in \Theta} \min\_{u \in U} \max\_{i \in I} \mathbb{g}\_i(\mathbf{x}, u, d\mathbf{es}, \theta) \le 0 \tag{3}$$

The max-min-max problem in Equation (3) mathematically expresses the feasibility question "For all the uncertainty realizations Θ, does there exist a control action *u* such that the constraint set *g* is feasible?". Equation (3) was employed in a multiperiod design optimization problem, where the deterministic uncertain parameter *θ* was allowed to vary within a prespecified range [11]. The feasibility constraint then laid the foundation for the concept of feasibility index, *F*, proposed by Swaney and Grossmann (1985) [12], as given by Equation (4).

$$\begin{aligned} F &= \max \quad \delta \\ \text{s.t.} & \max\_{\theta \in \Theta} \min\_{u \in \mathcal{U}} \max\_{i \in I} g\_i(\mathbf{x}\_i u, d\mathbf{s}\_i \theta) \le 0 \\ T(\delta) &= \{ \theta \mid (\theta^{\text{num}} - \delta \Delta \theta^-) \le \theta \le \theta \mid (\theta^{\text{num}} + \delta \Delta \theta^+) \} \end{aligned} \tag{4}$$

where *T* is the hyperrectangle for the uncertain parameters, *δ* is the scaled parameter deviation, and the superscript *nom* denotes nominal conditions. Equation (4) is the first significant attempt to quantify the degree of flexibility of a process design, and has been exploited by numerous studies on design optimization and process operability. However, Equation (4) constitutes a nondifferentiable global optimization problem and is still quite challenging to solve. Therefore, it requires simplifying assumptions and approximations to maintain a tractable problem. Swaney and Grossmann (1985) [13] introduced a heuristic vertex search method and an implicit enumeration scheme for the special case where the critical uncertainty realizations are assumed to lie at the vertices of the hyperrectangle *T*(*δ*). Clearly, this assumption fails to hold when the feasible space of the design problem is non-convex. Grossmann and Floudas (1987) [14] relaxed this assumption by developing a mixed-integer nonlinear programming (MINLP) problem for the feasibility test presented by Equation (3). They further proposed an active constraint strategy for the solution of the resulting MINLP. The mixed-integer formulation also provides a systematic approach to consider all possible critical uncertainty realizations without exhaustive enumeration. The proposed formulation was used for synthesis of a heat exchanger network with uncertain stream flow rates and temperatures [15]. The case of linear constraints reduces to an MILP problem, for which global solution is attainable by standard branch and bound enumeration techniques [14,16,17]. Bansal et al. (2000) [18] developed a computationally efficient theory and algorithm based on multiparametric programming techniques for this special case of flexibility analysis problems. The authors derived explicit expressions for the flexibility index as explicit functions of the continuous design variables. Pistikopoulos and Grossmann (1988a, 1988b, 1988c) used the flexibility test with linear constraints for optimal retrofit design [19–22] and redesign under infeasible nominal uncertainties [23]. Although these approaches are effective and promising to handle the design uncertainty, they require solving nested optimization problems, which poses a major challenge to solve complex and large-scale problems in a reasonable time. Raspanti et al. (2000) [24] proposed replacing the complementarity conditions of the lower level optimization problems with a well-behaved, smoothed nonlinear equality constraints, namely Kreisselmeier and Steinhauser function [25] and Chen and Mangasarian smoothing function [26].

One of the common assumptions in these approaches is the known bounds of the uncertainties, which is rarely the case in real world industrial applications. Pistikopoulos and Mazzuchi (1990) [27] and Straub and Grossmann (1990, 1993) [28,29] extended the flexibility test by assuming a probability distribution model for the parameter uncertainty, which improved the economic performance of the design optimization problem by addressing the "conservativeness" of the solution.

Another common assumption of these approaches is the steady-state operation of the plant design, which creates a significant limitation on the applicability of the methodologies. Although steady-state assumption holds true for the dominant life cycle of the plant operation, design optimization problem may fail to ensure the operability under transient behaviors such as startup or shutdown and transitions between different operating conditions. Dimitriadis and Pistikopoulos (1995) [30] proposed a dynamic feasibility index for the systems that are described by differential algebraic equations (DAE) subject to time-varying constraints. However, the time-dependent uncertainty in their formulation dictates to solve infinitely many dynamic optimization problems. Therefore, the authors assumed that the critical scenarios of uncertainties are known and lie on the vertices of the time-varying uncertainty space, similar to Swaney

and Grossmann (1985) [12]. The simplifying assumption reduced the problem to the form given by Equation (5).

$$\begin{aligned} DF(\text{des}) &= \max\_{\delta, \mu(t), t} \quad \delta \\ &\text{s.t.} \quad \dot{\mathbf{x}} = f(\mathbf{x}(t), \boldsymbol{\mu}(t), \text{des}, \theta(t), t), \quad \mathbf{x}(0) = \mathbf{x}\_0 \\ &\quad \quad g(\mathbf{x}(t), \boldsymbol{\mu}(t), \text{des}, \theta(t), t) \le 0 \\ &\quad \quad \quad \theta(t) = \theta^N(t) + \delta \Delta \theta^\varepsilon(t) \\ &\quad \quad \delta \ge 0, \quad \underline{\boldsymbol{\mu}}(t) \le \boldsymbol{u}(t) \le \overline{\boldsymbol{u}}(t) \end{aligned} \tag{5}$$

where the time dependence of the variables constitutes a dynamic optimization problem, and the solution was determined by control vector parameterization techniques [30]. Dynamic flexibility has been widely used in numerous design optimization applications including batch processes [31], separation systems [32–36], reaction systems [37], and heat exchanger network synthesis [38–40].

The dynamic assessment of the plant feasibility under uncertainty has been also studied through exploiting the multiperiod design optimization formulation proposed by Halemane and Grossmann (1983) [11]. Varvarezos et al. (1992) [41] implemented an outer-approximation approach to solve the multiperiod multiproduct batch plant problems operating with single product campaigns, which was formulated as an MINLP. Pistikopoulos and Ierapetritou (1995) [42] considered stochastic process uncertainty and proposed a two-stage decomposition that can handle convex nonlinear problems.

As presented in this section, the early studies on integrated design optimization have primarily focused on (i) investigating the range of operation (flexibility) of a nominal design configuration under foreseeable conditions, and (ii) determining the "best" possible trade-off between the investment cost and the capability of handling variations in the internal and external operating conditions. These studies mostly considered open loop processes, under the traditional assumption that controller design is a sequential task to process design. However, most processes in industry are operated in closed loop, and the controller schemes inherently alter the process dynamics, rendering the open loop flexibility analyses of lesser relevance. In other words, an "attainable" operating point according to open loop flexibility analysis may actually be an infeasible point in closed loop. Realizing the shortcomings of open loop flexibility analyses, researchers began investigating the "controllability" of process systems, and the interdependence of process control and design decisions. In the following section, we present a retrospective background on the integration of process control in the design optimization problem.

#### **3. Integration of Process Control in Design Optimization**

The initial efforts towards the integration of process control and design problems established a fundamental understanding on the interdependence of the two decision making mechanisms. The most pronounced school of thought in the early years to evaluate the controllability of the process design is "dynamic resilience", as conceptually defined by Morari (1983a, 1983b) [43,44].

Morari (1983) [43] described dynamic resilience as "the ability of the plant to move fast and smoothly from one operating condition to another and to deal effectively with disturbances". This depiction implies that there is not a clear-cut distinction between flexibility, which was discussed in Section 2, and resilience. However, Grossmann and Morari (1983) [45] pointed out the main difference as "resiliency refers to the maintenance of satisfactory performance despite adverse conditions while flexibility is the ability to handle alternate (desirable) operating conditions". This distinction is the primary motive for most of the flexibility analyses to study steady-state operations, while the resilience deals with the dynamic operations, as we will discuss in this section.

Dynamic resilience, as described by Morari (1983) [43], aims to find the "perfect controller" that is allowed by the physical limitations of the system to assess the controllability of the process by using the

internal model control (IMC) structure. The proposed technique decomposes the system transfer function *<sup>G</sup>*˜ into (i) a non-singular matrix *<sup>G</sup>*˜<sup>−</sup> to design the perfect controller *<sup>G</sup>*˜ <sup>−</sup><sup>1</sup> <sup>−</sup> , and (ii) a singular matrix *<sup>G</sup>*˜<sup>+</sup> to generate dynamic resilience indices based on (i) bounds on control variables, (ii) presence of right half plane transmission zeroes, (iii) presence of time delays, and (iv) plant-model mismatch. The proposed indices were used to improve the operability of numerous process, including heat integrated reactor networks [46–48], separation systems [49], heat exchanger networks [50].

Among the four aforementioned resilience indices, Perkins and Wong (1985) [51] studied the last two by adapting the "functional controllability" theorem proposed by Rosenbrock (1970) [52]. The authors further define a system to be functionally controllable if there exists a manipulated action *u*(*t*) that can generate any process output *y*(*t*) at any time *t*. Psarris and Floudas studied the dynamic resilience and functional controllability of multiple input multiple output (MIMO) closed-loop systems with time delays [53–55], and transmission zeroes [54,55]. Barton et al. (1991) [56] investigated the open loop process indicators, namely minimum singular value and right half plane zeroes, to assess the interactions between different design configurations and their operability with the best possible control configurations.

In the context of simultaneously assessing the process controllability in process design, one of the first significant contributions is the "back-off approach" introduced by Narraway et al. (1991) [57]. Narraway and Perkins (1994) [58] used this approach to systematically assess the trade-offs between all possible controlled and manipulated variable pairs in a mixed-integer formulation. Bahri et al. (1995) [59] employed the back-off approach to handle process uncertainties in an optimal control problem. The proposed approach is applicable to design linear and mildly nonlinear processes, and relies on three key steps, namely (i) perform a steady-state nonlinear process optimization, (ii) linearize the process at the optimum point, and (iii) "back-off" from the optimal solution by some distance to ensure the feasibility of the operation under some structured disturbance profile. The proposed approach was shown to be effective to select between alternative flowsheets as well as alternative control structures.

With the burgeoning interest in exploring the simultaneous design and control problem, the International Federation of Automatic Control (IFAC) organized the first workshop on "Interactions between Process Design and Process Control" in the Center for Process Systems Engineering at Imperial College London in 1992. The workshop laid the groundwork for a plethora of approaches with a wide range of diversity. Walsh and Perkins (1992) [60] implemented a PI loop in the flexibility analysis, where the input–output loop is selected by an exhaustive screening procedure. Luyben and Floudas (1992) [61] formulated a multiobjective MINLP problem to simultaneously consider the disturbance rejection capacity of the control loop through disturbance condition number and relative gain array to evaluate the interactions between the inputs and outputs of a MIMO system, while designing the process. Shah et al. (1992) [62] used the State-Task Network (STN) representation [63] to simultaneously consider the scheduling and design problems in a batch plant. Thomaidis and Pistikopoulos (1992) [64] introduced a framework to consider the design problem simultaneously with (i) the process flexibility through stochastic flexibility index, (ii) the effect of equipment failures to the overall performance by combined flexibility-reliability index, and (iii) the impact of equipment availability by combined flexibility-reliability index. These aforementioned novel approaches were shown to be promising concepts and techniques to address multiple facets of operational decisions simultaneously with the process design problem. As a result, succeeding studies after this workshop expanded these techniques and branched out to explore further opportunities.

Integrating PI controllers in the design optimization problem was one of the prominent outcomes of the workshop and became the most attractive option for the following research. The literature on PI controllers was already abundant and well-established by the time. Moreover, the explicit form of the controller structure made the integration relatively easy and intuitive, which significantly accelerated the research in closed-loop design optimization. Walsh and Perkins (1994) [65] presented an integrated

PI control scheme and process design for wastewater neutralization. Although the proposed approach was effective for the SISO process, it was reported that it entails further challenges for more complex processes. One major drawback of PI control is its inability to tackle MIMO systems without any advanced modifications in the feedback loop structure. Narraway and Perkins (1993, 1994) [58,66] developed an MILP-based formulation to systematically evaluate the economic performance of every input–output pair combination. Luyben and Floudas (1994a, 1994b) [67,68] adapted a similar approach in a multiobjective framework to determine the best performing input–output pair based on the controllability indices introduced by them, earlier (1992) [61]. The proposed framework was showcased on the design of a heat integrated distillation system [67] and a reactor-separator-recycle system [68]. Mohideen et al. (1996) [32] formulated a multiperiod design and control problem to account for the dynamic variations in the operation, while including the input–output pairing superstructure in the problem. Moreover, the authors used the flexibility index to account for the uncertain parameters in the model and presented a decomposition algorithm for the resulting complex problem. Bansal et al. (2000) [69] constructed a similar formulation as a mixed-integer dynamic optimization (MIDO) problem, which was solved by a Generalized Benders Decomposition (GBD)-based algorithm. The MIDO formulation was presented as follows.

min *<sup>u</sup>*,*des* ∑ *i*∈*NS wiC* - *x*˙ *i* (*t*), *x<sup>i</sup>* (*t*), *u<sup>i</sup>* (*t*), *des<sup>i</sup> s*.*t*. *x*˙ *i* (*t*) = *hd* - *xi* (*t*), *u<sup>i</sup>* (*t*), *des<sup>i</sup>* , *θi* , *t* , *x*(*t*) = *x*<sup>0</sup> *yi* (*t*) = *ha* - *xi* (*t*), *u<sup>i</sup>* (*t*), *des<sup>i</sup>* , *θi* , *t g* - *x*˙ *i* (*t*), *x<sup>i</sup>* (*t*), *y<sup>i</sup>* (*t*), *u<sup>i</sup>* (*t*), *des<sup>i</sup>* , *θi* , *t* ≤ 0 (6)

where *wi* is the discrete probability of a scenario *i* and *NS* is the discretized set of scenarios. The discretization of uncertainty in the process was first proposed by Grossmann and Sargent (1978) [9].

Although the aforementioned PI-based design and control frameworks are applicable on nonlinear processes, the range of operability is usually limited due to the mismatch between the nonlinear process model and the linearized control model. Ricardez-Sandoval et al. (2008, 2009) [70,71] used robust control tools and the back-off approach to integrate PI control and ensure its stability while solving the design optimization problem. The proposed approach was also tested against the Tennessee Eastman Process [72]. The back-off approach was later generalized for control structure selection in nonlinear processes by Kookos and Perkins (2016) [73]. Ricardez-Sandoval & co-workers have extensively studied back-off approach for simultaneous process design and control under uncertainty [74–76].

One main limitation of integrating PI control in the design optimization in a dynamic formulation is the increasing problem size and complexity. Kookos and Perkins (2001) [77] developed an algorithm for the integrated PI control and design optimization problem, where the size of the search space was reduced systematically in each successive iteration. Malcolm et al. (2007) [78] proposed an "embedded control optimization" procedure, where the authors introduced a two-stage decomposition scheme that approximates the complete integrated problem. The proposed approach reduced the problem size and complexity, and was showcased on larger scale problems including a reactor-separator system [79].

Apart from the inability to naturally handle MIMO systems, PI controllers do not explicitly account for any process constraints stemming from operational, environmental, and safety limitations. Model predictive control (MPC) overcomes these shortcomings by postulating a constrained dynamic optimization problem subject to an explicit model of the process [80]. One of the first remarkable efforts to integrate an MPC scheme in a nonlinear design problem was published by Brengel and Seider (1992) [81].

Here, the authors postulate a bi-level optimization problem, where the *leader* has an economic objective, while the *follower* is the MPC formulation, as presented by Equation (7).

min *des Cdes*- *des* + *κC<sup>κ</sup>* - *x*(*t*), *y*(*t*), *u*(*t*), *des*, *θ*(*t*) *s*.*t*. *fdes*- *des*, *θ*(*t*) = 0 *gdes*- *des*, *θ*(*t*) ≤ 0 min *u*(*t*) *Cu* - *x*(*t*), *y*(*t*), *u*(*t*), *des*, *θ*(*t*) *s*.*t*. *x*˙ = *fu* - *x*(*t*), *y*(*t*), *u*(*t*), *des*, *θ*(*t*) *gu* - *x*(*t*), *y*(*t*), *u*(*t*), *des*, *θ*(*t*) = 0 *hu* - *x*(*t*), *y*(*t*), *u*(*t*), *des*, *θ*(*t*) ≤ 0 (7)

where *κ* is the design and control integration parameter that scales the trade-off between the controllability of the system and the investment cost. The bi-level problem presented in Equation (7) is challenging to solve without appealing to simplifications. Therefore, the authors proposed replacing the follower problem by complementary slackness equations. However, the solution strategy was still intractable for more complex systems due to the numerical calculation of the second derivatives [81]. As a consequence, integration of the MPC scheme in the design optimization had been rather limited in the literature for almost a decade, until the invention of multiparametric MPC (mpMPC/explicit MPC).

Bemporad et al. (2002) [82] proposed formulating the MPC problem as an explicit function of the initial conditions of the system. This novel strategy allowed for deriving piecewise affine explicit control laws by treating the initial conditions as parameters. The proposed approach formulated the explicit MPC problem as presented by Equation (8).

$$\begin{aligned} \left\| u\_{t}(\theta) \right\| &= \operatorname\*{arg\,min}\_{u\_{t}} \left\| x\_{N} \right\|\_{P}^{2} + \sum\_{t=1}^{N-1} \left\| x\_{t} \right\|\_{Q}^{2} + \sum\_{t=1}^{N-1} \left\| y\_{t} - y\_{t}^{sp} \right\|\_{QR}^{2} + \sum\_{t=0}^{M-1} \left\| u\_{t} - u\_{t}^{sp} \right\|\_{R}^{2} + \sum\_{t=0}^{M-1} \left\| \Delta u\_{t} \right\|\_{R1}^{2} \\ \text{s.t.} \quad \left\| x\_{t+1} = A\mathbf{x}\_{t} + Bu\_{t} + \mathbb{C}d\_{t}, \quad y\_{t} = D\mathbf{x}\_{t} + Eu\_{t} + Fd\_{t} \\ \underline{\mathbf{x}}\_{t} \le \mathbf{x}\_{t} &\leq \overline{\mathbf{x}}\_{t}, \quad \underline{y}\_{t} \le \mathbf{y}\_{t} \le \overline{\mathbf{y}}\_{t}, \quad \underline{u}\_{t} \le u\_{t} \le \overline{u}\_{t}, \quad \underline{\Delta u}\_{t} \le \Delta u\_{t} \le \overline{\Delta u}\_{t}, \quad \underline{d}\_{t} \le d\_{t} \le \overline{d}\_{t} \\ \theta &= [\mathbf{x}\_{t+0}, u\_{t-1}, d\_{t}, y\_{t}^{sp}, u\_{t}^{sp}]^{T} \end{aligned} \tag{8}$$

where *N* is the prediction horizon, *M* is the output horizon, superscript *sp* denotes set point, *Q*, *QR*, *R*, and *R*1 are the corresponding weight matrices determined by tuning, *P* is calculated by discrete algebraic Riccati equation, and ·*<sup>ψ</sup>* denotes weighted vector norm with a weight matrix *ψ*. Different than conventional MPC, Equation (8) formulates the optimal control problem exactly and completely offline as a function of the set of parameters *θ*. The solution of this problem can be determined by multiparametric programming techniques, which express the solution space as a piecewise affine function, as presented by Equation (9).

$$\begin{aligned} \mu\_t(\theta) &= \mathbb{K}\_{\mathbb{R}} \theta + r\_{\mathbb{R}}, \quad \forall \theta \in \mathbb{C} \mathbb{R}\_{\mathbb{R}}\\ \mathbb{C} \mathbb{R}\_{\mathbb{R}} &:= \{ \theta \in \Theta \mid \mathbb{C} \mathbb{R}^A \theta \le \mathbb{C} \mathbb{R}^b \}, \quad \forall n \in \{ 1, 2, \dots, N \mathbb{C} \} \end{aligned} \tag{9}$$

where *CRn* is referred as a critical region and it is the active polyhedral partition of the feasible parameter space, Θ is a closed and bounded set, and *NC* is the number of critical regions.

The control law given by Equation (9) reduces the complexity of solving an online optimization problem to a simple look-up table algorithm (also known as point location problem) and function

evaluation, all of which are affine operations. Hence, the complexity of implementing an MPC scheme is similar to that of a PI controller.

Sakizlis et al. (2003) [83] exploited the explicit nature of the mpMPC solution in the context of design and control integration. The authors formulated a bi-level mixed-integer dynamic optimization problem similar to Equation (7), where the leader accounted for the investment and operating costs in the objective function subject to the dynamic high-fidelity model, and the follower MPC problem was substituted by the affine control law Equation (9). The proposed formulation offered an elegant and systematic methodology to reduce the complexity of the bi-level Equation (7) into a single level dynamic optimization problem. However, the solution strategy still required repetitive linearizations and solving a multiparametric programming problem at every iteration, which can be restrictive for large-scale complex problems. Diangelakis et al. (2017) [84] alleviated that limitation by deriving a "design dependent offline controller", which allowed for solving a single MIDO problem after integrating the control law in the high-fidelity model. Eliminating the linearization step and formulating a single synergistic design and control problem also improved the economic performance of the resulting process compared to the approach proposed by Sakizlis et al. (2003) [83]. The proposed formulation was also showcased on a tank, a continuous stirred tank reactor, and a residential scale combined heat and power unit. The cost effectiveness of the MPC integrated optimal design was also reported to be superior than PI integrated approaches in the literature. Diangelakis and Pistikopoulos (2017) [85] reported that the mpMPC integrated optimal combined heat and power unit operated more fuel efficient in closed loop than PI integrated design. Similarly, Sanchez-Sanchez and Ricardez-Sandoval (2013) [86] showcased a system of CSTRs, where the MPC integrated framework reduced both the operating and the investment costs compared to the PI control integrated approach.

One common aspect of the studies on simultaneous design and control optimization is the assumption that the process will be operated around the same steady-state point throughout the entire life cycle of the plant. However, the external plant conditions, such as market conditions, may dictate a considerably wider operating region with multiple steady-state points [1]. The increasing competition among the businesses impacts the volatility of the market, which creates rapid fluctuations in the energy and raw material prices as well as their availability. Moreover, the demand rate on the product is also subject to considerable variations during the plant operation. Therefore, it is clear that there exists a "best" operating strategy under the knowledge available to the operator, which necessitates the operability of the plant across a wider range. For example, high production rates may be less profitable during the night time because of increased energy prices and hence, operating the energy intensive processes during the daytime may reduce the operating costs. This indicates that the operating level of a processing unit might vary drastically by the choice of the operator. However, the integrated design and control frameworks discussed in this section usually assume a single operating point around which a controllability and flexibility analysis is conducted. Consequently, these frameworks do not attempt to provide any means of guaranteeing the operability of the process at different regions. In the next section, we will discuss several approaches that account for multiple operating regions in a plant, and their scheduling during the operational optimization.

#### **4. Towards the Grand Unification of Process Design, Scheduling, and Control**

Process design, scheduling, and control problems are traditionally constructed to address different objectives and they span widely different time scales. In a nutshell, the plant design problem dictates the capacity of processing and it usually comprises the most uncertainty due to its years long lifecycle. The scheduling problem addresses the allocation of the resources and time, as well as the operating level of processing units and their maintenance based on some economic criteria over days/months long horizons. Lastly, the control problem maintains the performance of the plant, while satisfying any physical

limitations such as the environmental and safety constraints. The discrepancy in the objectives and time scales creates a challenging problem to systematically evaluate and determine the optimal trade-off between different decision makers.

Process scheduling is more critical in batch operations than continuous operations, as the former are inherently dynamically operated. Accordingly, the initial efforts focused primarily on the batch processes for the integration of the operational optimization and design problems. Birewar and Grossmann (1989) [87] formulated NLP models to incorporate the scheduling decisions in the batch sizing and timing problem in a multiproduct plant for unlimited intermediate storage and zero wait policies. Shah et al. (1992) [62] tackled a similar problem by using the STN representation. White et al. (1996) [88] investigated the switchability of continuous processes between different operating points through formulating an optimal control problem that accounts for the terminal criteria and path constraints within a range of design parameters. Bhatia and Biegler (1996, 1997) [89,90] formulated a dynamic optimization problem, where an economic objective function was subject to a dynamic high-fidelity model of the process described by differential algebraic system of equations. The authors proposed a solution strategy based on discretizing the process model by orthogonal collocation over finite elements, followed by solving the resulting NLP by using a standard solver. The proposed modeling and solution strategy was shown to be promising to satisfy the path constraints, which is a crucial benefit for dynamic systems. Terrazas-Moreno et al. (2008) [2] extended this integration approach to account for the binary decisions in the scheduling problem by formulating a MIDO. Similar to Bhatia and Biegler (1996, 1997) [89,90], the authors first discretized the problem by orthogonal collocation, followed by solving the resulting MINLP.

The early studies that explore the interactions between the scheduling and process control decisions have a significant role in shaping today's approaches for the integrated design optimization problem. In their excellent review article, Baldea and Harjunkoski (2014) [91] classified these attempts to integrate the scheduling and control decisions as (i) "top-down approaches", where the process dynamics and control elements are incorporated in a scheduling skeleton, and (ii) "bottom-up approaches", where the process economics are implemented in the plant-wide control decisions.

In terms of characterizing the transitions between different products in a single operating unit, Mahadevan et al. (2002) [92] introduced a unique "top-down" perspective on the operational optimization problem, revealing that a simultaneous approach on the scheduling and control problem can identify and eliminate the fundamental limiting behavior during the transitions, as showcased on a polymer grade transition process. However, the presented approach requires case specific heuristic decisions to select the "best" fitting scheduling and control configuration and hence, it is not suitable for different applications in the general sense. Chatzidoukas et al. (2003) [93] studied a similar polymerization reactor, and formulated a MIDO problem to determine the time optimal transition between different polymer grades and best performing control structure simultaneously. Flores-Tlacuahuac and Grossmann (2006) [94] introduced a monolithic approach on a multiproduct cyclic CSTR, where the profit was maximized by manipulating the production sequence, transition times, production rates, length of processing times, and amounts manufactured of each product. In contrast to the earlier studies [92,93], the authors focused on the manipulated actions rather than the optimal control configuration. They formulated a MIDO problem, which was solved by discretization of the differential algebraic equations by orthogonal collocation on finite elements followed by solving the resulting MINLP. The presented approach has been extensively studied in the following years to broaden its scope and effectiveness. Terrazas-Moreno et al. (2007) [95] applied this approach on two industrial polymerization reactors. Terrazas-Moreno et al. (2008) [2] formulated a design optimization problem accounting for the scheduling and open loop control trajectories using this approach. Flores-Tlacuahuac and Grossmann (2010, 2011) extended the formulation to partial differential equation systems, and showcased on tubular reactors with single [96] and multiple production lines [97].

This monolithic approach usually generates open loop control trajectories, i.e., no feedback loop is assumed to develop the input and output profiles. However, the processing units are subject to internal process disturbances, and the mismatch between the process and the model leads to deviations in the targeted operations. Zhuge and Ierapetritou (2012) [98] implemented the monolithic approach in closed loop, where the authors initiate a readjustment procedure to solve the integrated problem online if the states deviate from their reference trajectories. This approach does not completely resolve the issue of handling the process disturbances or the process/model mismatch; however, it was shown to mitigate these concerns to a great extent. Gutiérrez-Limón et al. (2014) [99] also implemented a similar closed-loop strategy with a nonlinear model predictive control scheme, while extending the scope of the problem statement to account for an extended horizon production policy. However, both approaches require solving a complex and large-scale MINLP problem at the time steps of the controller, which makes it unsuitable for the processes with fast dynamics.

Low-order representation of fast process dynamics in the scheduling problem has been an effective approach to reduce the computational burden of solving complex optimization problems. Du et al. (2015) [100] proposed a time scale-bridging model that describes the closed-loop input–output behavior of a process in the scheduling formulation, postulated as a MIDO problem. The low-order representation also maintains the stability of the process in the existence of process/model mismatch and handles disturbances. Baldea et al. (2015) [101] extended this approach to MPC governed systems.

Burnak et al. (2018) [102] also addressed the online computational burden of "top-down" approaches by developing a multiparametric programming-based approach, where the authors explicitly mapped (i) the closed-loop dynamic process behavior in a "control-aware" scheduling problem, and (ii) the continuous and binary scheduling level decisions such as the operating level and operational mode of the system in a "schedule-aware" MPC scheme (iii) to yield the optimal operational decisions. The offline nature of the integrated scheduling and control scheme allows for determining the feasible operating space prior to actualizing the operation. Furthermore, reducing the problem complexity from solving online optimization problems to a simple look-up table and affine function evaluation, the framework is well-suited for fast process dynamics. Charitopoulos et al. (2019) [103] employed a similar multiparametric programming approach to include the planning decisions in their framework.

In the "bottom-up" approaches, on the other hand, incorporating the economic objectives in the plant control structures has been perceived as the key for seamless integration of scheduling and control. For this purpose, MPC formulations provide the flexibility to account for a spectrum of objectives in the control level due to their optimization-based structures. Loeblein and Perkins (1999) [104] presented an economic analysis of unconstrained MPC scheme operating under constrained systems. The authors determined the most cost-effective model predictive regulatory control structure by using the back-off approach to satisfy the constraints. Zanin et al. (2002) [105] addressed the discrepancy between the real-time optimization (RTO) and control layers by incorporating the economic optimization problem in the controller and feeding the same piece of information in both layers. The proposed formulation diminishes the discrepancy between the decision layers to yield more economical operations, but the resulting control scheme does not guarantee the stability of the process for the entirety of operations. Rawlings and Amrit (2009) [106] developed asymptotic stability criteria by formulating the so-called "economic MPC" (or EMPC), where the objective function of the MPC is designed to minimize the operational costs instead of maintaining the steady state of the process. This approach aims to replace the conventional two-layer structure with RTO and dynamic regulatory control by a single control layer, where the economic optimization and process regulation are conducted simultaneously. Amrit et al. (2011) [107] further extended the stability criteria by (i) imposing a region constraint on the terminal state instead of a point constraint, and (ii) adding a penalty on the terminal state to the regulator cost.

Similar to the monolithic "top-down" scheduling and control approach, EMPC has been shown to be too complex to be solved in the control time steps. This limitation has led the researchers to develop decomposition algorithms for faster computational times. Würth et al. (2011) [108] proposed a decomposition framework for the single layer dynamic RTO formulation, where the slow trends and process uncertainty is handled in the upper layer, while the lower layer accounts for the fast disturbances acting on the process. Ellis and Christofides (2014) [109] focused on selecting a suitable input configuration for such two-layered dynamic RTO structures such that the asymptotic stability is guaranteed. Jamaludin and Swartz (2017) [110] and Li and Swartz (2019) [111] employed a convex MPC problem in the lower level regulatory control, which enabled its exact substitution with KKT optimality conditions. Simkoff and Baldea (2019) [112] used the same substitution strategy on a production scheduling problem.

Design optimization accounting for the scheduling and control decisions with closed-loop implementation is relatively recent in the literature. Patil et al. (2015) [3] modeled the product transitions in design optimization, while maintaining the stability of the closed-loop system governed by a PI control scheme. The authors formulated an MINLP similar to Equation (6) with the contribution of the criterion, *eig*(*A<sup>z</sup> <sup>i</sup>* (*xlin*)) < 0, which enforces the stability of the linearized states for all products *i* in a multiproduct unit under all critical scenarios *z*. Due to the linearization of the controllers around the operating point, this approach requires repetitive identification of the states at every optimization iteration. Koller and Ricardez-Sandoval (2017) [4] improved this approach by applying orthogonal collocation on finite elements on the integrated problem, and Koller et al. (2018) [5] employed the back-off method to satisfy the constraints under uncertainty by using Monte Carlo sampling techniques to determine the back-off terms.

Recently, Burnak et al. (2019) [6] introduced a multiparametric programming-based theory and framework for the integration of process design, scheduling, and control. The authors derived offline design dependent control and scheduling schemes that can be incorporated in a MIDO formulation in a multi-level fashion, as presented by Equation (10).

min *<sup>u</sup>*,*s*,*des <sup>τ</sup>* 0 *C*(*x*(*t*), *y*(*t*), *u*(*t*),*s*(*t*), *des*, *d*(*t*))*dt s*.*t*. *x*˙(*t*) = *f*(*x*(*t*), *y*(*t*), *u*(*t*),*s*(*t*), *des*, *d*(*t*), *t*) *y* ≤ *y*(*t*) = *g*(*x*(*t*), *y*(*t*), *u*(*t*),*s*(*t*), *des*, *d*(*t*), *t*) ≤ *y x* ≤ *x*(*t*) ≤ *x*, *des* ≤ *des* ≤ *des*, *d* ≤ *d*(*t*) ≤ *d st*(*θs*) = arg min*<sup>s</sup>* ∑ *ts*∈*Ns Cs*(*xts*, *yts*,*sts*, *des*, *dts* ) *s*.*t*. *xts* ≤ *xts*<sup>+</sup><sup>1</sup> = *Ats xts* + *Bts sts* + *Ctsdts* ≤ *xts yts* ≤ *yts* = *Dts xts* + *Ets sts* + *Ftsdts* ≤ *yts sts* ≤ *sts* ≤ *sts*, *dts* ≤ *dts* ≤ *dts <sup>θ</sup><sup>s</sup>* <sup>≤</sup> *<sup>θ</sup><sup>s</sup>* = [*x<sup>T</sup> ts*<sup>=</sup>0, *<sup>y</sup><sup>T</sup> ts*<sup>=</sup>0, *dts*, *des*] *<sup>T</sup>* <sup>≤</sup> *<sup>θ</sup><sup>s</sup> ut*(*θc*) = arg min*<sup>c</sup>* ∑ *tc*∈*Nc Cc*(*xtc* , *ytc* , *utc* , *des*, *dtc* ) *s*.*t*. *xtc* ≤ *xtc*<sup>+</sup><sup>1</sup> = *Atc xtc* + *Btcutc* + *Ctcdtc* ≤ *xtc ytc* ≤ *ytc* = *Dtc xtc* + *Etcutc* + *Ftcdtc* ≤ *ytc utc* ≤ *utc* ≤ *utc* , *dtc* ≤ *dtc* ≤ *dtc <sup>θ</sup><sup>c</sup>* <sup>≤</sup> *<sup>θ</sup><sup>c</sup>* = [*x<sup>T</sup> tc*<sup>=</sup>0, *<sup>y</sup><sup>T</sup> tc*<sup>=</sup>0, *dtc* , *des*] *<sup>T</sup>* <sup>≤</sup> *<sup>θ</sup><sup>c</sup>* (10)

where *s* and *u* denote the scheduling and control decisions, respectively. Note that the proposed formulation postulates explicit expressions for the scheduling and control strategies as functions of a set of parameters, *θ*, which includes the design of the process. The design dependence of the operational strategies allows for their direct integration in the MIDO formulation. The postulated formulation has two main benefits, (i) due to the explicit form of the follower problems, the multi-level MIDO problem is reduced to a single level, and (ii) only the design variables are left as the degrees of freedom of the problem, since the remaining are determined as a function of the design.

#### **5. Current Challenges and Future Directions**

The PSE community has achieved unequivocally remarkable progress in realizing and advancing the set goals of Professor Sargent on systematic design optimization in five decades. Today, using design optimization tools to at least some extent has long become the standard practice in many industries. Commercial modeling and simulation software tools such as gPROMS (https://www.psenterprise.com/products/gproms) and Aspen Plus Dynamics (https://www.aspentech.com/en/products/pages/aspen-plus-dynamics) have been featuring robust and efficient solvers for dynamic optimization problems for a few years. Despite these milestones in PSE, we still must make significant assumptions and simplifications regarding the operational decisions in the process design phase, even though the impact of their interdependence on process economics and operability has been articulated in numerous studies. Hence, the academia still needs to mature the theoretical foundations and the applicability of unified design optimization approaches before it gains wide industrial recognition. Here, we discuss some of the bottlenecks and potential directions to improve the state-of-the-art for industrial practice.

#### *5.1. The Need for an Industrial Benchmark Problem*

As we have presented in this paper, there is a plethora of proposed modeling techniques and solution approaches for the next generation unified design optimization problems. Therefore, it is clear that we need a generally accepted benchmark problem, preferably in industrial scale, to validate the effectiveness of proposed methodologies. The PSE community has benefited greatly from such standardized problems, such as the famous Tennessee Eastman Process detailed by Downs and Vogel (1993) [113] for process control studies. We believe that a well-defined problem will clarify the objectives in unified design optimization and accelerate the research towards industrial expectations. The problem should describe at least the following.


3. *Product demand and availability of the utility, raw materials, and operating units over a time horizon*. Production allocation and timing is a key aspect of scheduling problem, which are heavily dictated by the product demand and availability of resources. However, it is not a trivial practice to estimate the future of these quantities. Therefore, probability distributions of these components will be beneficial to determine their expected values, while being able to take into account their worst-case scenarios.

#### *5.2. Robust Advanced Control and Scheduling Strategies*

Incorporation of advanced control schemes seamlessly in the design optimization problem requires the controller to capture the dynamics of the process for the entire range of design variables. Burnak et al. (2019) [6] attempted to approximately model the design configuration as a right-hand uncertainty in the constraint set, validated by closed-loop simulations and closed-loop MIDO problems. However, the design variables impose uncertainty in the left-hand side of the constraints, as well as the nonlinear and bilinear terms in the objective function. Therefore, robust control strategies need to be developed for accurate predictions of future states in the control level prior to the realization of the design, and to guarantee the stability of the closed-loop operations in simultaneous approaches.

Analogously, scheduling schemes should be robustified in the design optimization to minimize the rescheduling due to unexpected disruptive events, such as unit failure, drastic changes in product demand rate and raw material availability. Excluding these events in the scheduling scheme may result in steep changes in the target operation, and thus unattainable set points for the controller.

#### *5.3. Considering Flowsheet Optimization, Process Intensification, and Modular Design Opportunities*

Optimization-based plant design techniques have been used and developed for more than four decades [114,115]. These techniques postulate "superstructures" that systematically simulate and compare every combination of flowsheet possibility to determine the optimal process. More recently, superstructures have been formulated at the phenomena level to capture the fundamental relations between the mass and energy, which in turn yields intensified processes [116–122]. Such intensified processes are expected to deliver significantly increased operational efficiency and decreased unit volumes, making them very attractive options both in academia and industry [123]. This rapidly growing interest in intensified processes is one of the most pronounced directions that the PSE community has been taking. Therefore, studying these intensified processes in the context of unified design optimization will attract a wider audience from the industry. Clearly, modeling the spatial (synthesis/intensification) and temporal (scheduling/control) decisions simultaneously in a single problem formulation will capture even more synergistic interactions, which will increase the process profitability.

Furthermore, the researchers studying process intensification can benefit from the tools and methodologies on unification of design, scheduling, and control. Baldea (2015) [124] reported a theoretical justification for the loss of control degrees of freedom due to process intensification, which poses a significant limitation on intensification activities. Tian and Pistikopoulos (2019) [125] and Dias and Ierapetritou (2019) [126] discuss the limitations on the operability of such intensified systems and potential directions to overcome these limitations in their excellent review papers. The researchers on process intensification technologies can adopt the techniques, ranging from steady-state and dynamic flexibility to integration of scheduling and control decisions, in order to address the operability issues.

#### *5.4. Theoretical and Algorithmic Developments in MIDO*

The most limiting bottleneck of the simultaneous approaches is the size of the integrated MIDO problems. The time component of the problem significantly increases the computational complexity, yielding infinitely many NP-hard problems to acquire an optimal solution profile. However,

tailored algorithms can be developed by using the special structure of such integrated problems. For instance, the open loop design optimization problem is relatively simpler than the integrated MIDO, and constitutes a lower bound on the optimal solution of the overall problem. Such properties can be exploited in decomposing the MIDO into subproblems to significantly reduce the search space for faster algorithms.

#### *5.5. Software Development*

Despite the theoretical and practical advances in the unified design problem among the academia, there is no commercially available platform or a software prototype. Such a tool will make the integrated approaches more accessible to the process designers in industry who are not necessarily experts on process control and scheduling, and it will attract more researchers from different disciplines and backgrounds. Pistikopoulos et al. (2015) [127] introduced the PARametric Optimization & Control (PAROC) framework to design explicit controllers based on high-fidelity models, which can be a viable option to address the grand unification challenge [6,84,102,128,129]. However, it is clear that more progress is needed to engage a wider audience.

**Author Contributions:** Conceptualization, B.B. and E.N.P.; writing, B.B.; review, N.A.D. and E.N.P.; resources E.N.P.; supervision, E.N.P.

**Funding:** Financial support from the National Science Foundation (Grant No. 1705423) and Texas A&M Energy Institute, Shell Oil Company, Rapid Advancement in Process Intensification Deployment (RAPID) Institute, and Clean Energy Smart Manufacturing Innovation Institute (CESMII) is greatly acknowledged.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **References**







c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Optimization-Based Scheduling for the Process Industries: From Theory to Real-Life Industrial Applications**

#### **Georgios P. Georgiadis 1,2, Apostolos P. Elekidis 1,2 and Michael C. Georgiadis 1,2,\***


Received: 27 May 2019; Accepted: 4 July 2019; Published: 10 July 2019

**Abstract:** Scheduling is a major component for the efficient operation of the process industries. Especially in the current competitive globalized market, scheduling is of vital importance to most industries, since profit margins are miniscule. Prof. Sargent was one of the first to acknowledge this. His breakthrough contributions paved the way to other researchers to develop optimization-based methods that can address a plethora of process scheduling problems. Despite the plethora of works published by the scientific community, the practical implementation of optimization-based scheduling in industrial real-life applications is limited. In most industries, the optimization of production scheduling is seen as an extremely complex task and most schedulers prefer the use of a simulation-based software or manual decision, which result to suboptimal solutions. This work presents a comprehensive review of the theoretical concepts that emerged in the last 30 years. Moreover, an overview of the contributions that address real-life industrial case studies of process scheduling is illustrated. Finally, the major reasons that impede the application of optimization-based scheduling are critically analyzed and possible remedies are discussed.

**Keywords:** process scheduling; optimization; process system engineering; mixed-integer programming

#### **1. Introduction**

Scheduling is concerned with the allocation of scarce resources among competing activities over time. It is a decision-making process aiming to optimize one or more objectives by taking into account the processes taking place and their interactions with the environment. Scheduling problems exist in many manufacturing and production systems, in transportation and distribution of people and goods, and in other types of industries. The three elements which need to be mapped out are time, tasks and resources: The time at which the tasks have to be performed needs to be optimized considering the availability and restrictions on the required resources. The resources may include processing, material storage and transportation equipment, manpower, utilities (e.g., steam, electricity), any supplementary equipment and so on. The tasks typically include processing operations (e.g., reaction, separation, blending, packaging) as well as other activities like transportation, cleaning in place, changeovers, etc. Both external and internal elements of the production need to be considered. The external element originates from the need to co-ordinate manufacturing and inventory levels based on a given demand, as well as arrival time of raw materials and even maintenance activities. The internal element considers the execution of tasks in an appropriate sequence and time, while taking into account all external considerations and resource availabilities. Overall, the sequencing and timing of tasks over time and the assignment of appropriate resources to the tasks must be performed in an efficient manner, that

will, as far as possible, optimize an objective. Typical objectives include the minimization of cost or maximization of profit, the maximization of throughput, the minimization of tardy jobs, etc.

Flexible multipurpose plants are able to produce a wide range of different products using a variety of production routes. This characteristic makes such plants particularly effective for the manufacture of classes of products that exhibit a large degree of diversity and which are subject to fast-varying demands. Due to their inherent flexibility, the scheduling of such plants is a problem of high complexity. Compared to other parts of the supply chain management (e.g., distribution management and inventory control), the production scheduling is often by far the most computationally demanding part. The most general "multipurpose" plants can be viewed as collections of production resources (e.g., raw materials, processing and storage equipment, utilities, manpower) shared by a number of processing operations manufacturing a number of products over a given time horizon. The process may involve intermediates shared among two or more products, recycles of unconverted material, and multiple routes to the same end product. Single or multiple stage multi-product plants are thus special cases of multipurpose plants.

Roger Sargent was one of the first researchers in chemical engineering who foresaw the value of, and need for optimization in the design, control, and operation of process systems. One of the major steps in his research was in the area of multipurpose batch and continuous process scheduling, where the introduction of the state-task network (STN) concept was a major breakthrough. This generic representation allowed researchers to efficiently address arbitrary configurations of recipe-based batch operations. The major novelty was that equipment was not preassigned, like in previous contributions [1]. The utilization of the STN representations for the production scheduling problem resulted to a discrete time mixed-integer linear programming (MILP) model. Despite its novelty, that paper was rejected twice, since at that time MILP algorithms were considered inefficient and incapable of solving large-scale problems, which of course has changed drastically over the last decades [2]. Furthermore, progress with the STN modelling approach was also due to the improved formulation proposed by Nilay Shah in which big-M constraints were replaced by fewer and tighter sets of constraints [3].

An indicative example of the impact of Sargent's research work in the process scheduling area is his paper, "A General Algorithm for Short-Term Scheduling of Batch-Operations. 1. MILP Formulation" by E. Kondili, C. C. Pantelides, and R. W. H. Sargent, Computers & Chemical Engineering, 17, 211–227 (1993). This is one of the most widely cited contributions in the PSE (Process Systems Engineering) community (over 800 citations in SCOPUS as of April 2019), and has been recognized by researchers all over the world as the major framework for mathematically modeling batch and continuous operations through the state-task-network representation.

The formulation of Kondili et al. [4] relies on binary variables that specify whether a task starts in an equipment at the start of each time period. Other key variables denote the amount of material in each state, and the amount of utility required for processing tasks over each time interval. Equipment and utility usage constraints as well as material balances and capacity constraints are considered in the formulation. A common, discrete time grid is employed to capture all plant resource utilizations in a straightforward manner. This approach was hindered in its ability to handle large problems mainly due to the limitations of discrete-time approaches that require relatively large numbers of grid points, thus resulting to large-sized models.

Inspired by the work of Kondili et al. [1] a number of contributions appeared in the literature utilizing the STN representation and the overall mathematical programming framework to address general classes of batch and continuous process scheduling problems. Shah et al. [3] was able to generate the smallest possible integrality gap for this type of formulation by efficiently modifying the allocation constraints. They additionally proposed a tailored branch-and-bound solution procedure that uses a significantly smaller LP (Linear Programming) relaxation in order to further improve integrality at each node. In the same research, authors addressed the cyclic scheduling problem, where they simultaneously derived optimal schedules as well as the frequency at which they should be repeated [5]. Papageorgiou and Pantelides [6,7] further expanded this work to cover the case of multiple campaigns. Yee and Shah [8,9] also considered variable elimination to improve the performance of general discrete-time scheduling models. More specifically, they recognized that only about 5–15% of the variables are active at the integer solution, and it would be beneficial to identify and eliminate as far as possible inactive variables prior to solving the scheduling problem. To achieve that, they introduce an LP-based heuristic, alongside a flexibility and sequence reduction technique and a formal branch-and-price method.

Pantelides et al. [10] presented an STN-based approach for the scheduling problem of pipeless plants, where material is transferred between processing stations in vessels, thus requiring the simultaneous scheduling of the movement and processing operations. Pantelides [11] criticized the STN, arguing that despite its advantages, it inherently suffers from a number of drawbacks. For example, the fact that each equipment if treated as a distinct entity that results to solution degeneracy in case of multiple equivalent items exist. Therefore, he proposed a differentiated representation, the resource-task network (RTN), which is based on the equable description of all resources [11]. Contrary to the STN representation, where different states are consumed or produced by a task utilizing the equipment and the utilities, in this approach even the items of equipment or the plants' utilities are considered as resources. Production units are assumed to be consumed at the start and produced at the end of a task. Furthermore, different equipment conditions (e.g., "clean" or "dirty") can be treated as separate resources, with different activities (e.g., "processing" or "cleaning") consuming and generating them—this allows for a simple representation of changeovers. Pantelides illustrated that the integrality gap of RTN formulations is never worse than the most efficient form of STN formulation, and the ability to adapt additional problem features in a straightforward way, made it a favorable framework for future research.

The review above has mainly focused on the development of discrete-time models. As pointed out by Schilling [12], while discrete-time models have been capable to handle numerous industrially-relevant problems (see, e.g., [13]), they are characterized by a number of inherent drawbacks:


In order to address these issues a number of researchers have attempted to develop scheduling models that employ a continuous representation of time. As a result, fewer grid points are required leading to fewer variables and smaller model sizes.

Dimitriadis et al. [14] describe two rolling horizon procedures for medium-term planning and scheduling, based on the more general RTN formulation. They take advantage of the unique properties of Wilkinson et al. [15] and aggregation in this context. In the forwards rolling horizon algorithm, the horizon is divided into two-time blocks. The first is relatively short and modelled in detail, while the second is relatively long and modelled using the aggregate scheduling formulation. The solution of this MILP gives rise to a detailed solution for the first period and an aggregate one for the second. Dimitriadis et al. [14] recognized that, rather than fix all the variables in the first period at the next iteration of the procedure, it makes sense only to fix the complicating integer variables and leave the continuous ones free for further optimization. At the next iteration, there are three time blocks, the first one with fixed integer variables, the second one modelled in detail and the third (the remainder of the horizon) modelled at an aggregate level. The algorithm proceeds until a detailed solution is obtained for the entire horizon.

As noted in the excellent review by Shah [16] a common conclusion in most PSE contributions is that one of the most important advances in the area of process scheduling over the past 25 years has been the increasing usage of rigorous mathematical programming approaches. In addition, the importance of the establishment of frameworks for process scheduling which can be used for the description of a wide variety of processes and for the development of general solution algorithms has been emphasized.

The contributions described above, inspired many researchers from the PSE community to further investigate the production scheduling problem. Numerous novel approaches have been proposed by different research teams, providing novel efficient models and solution techniques. Network representations [17], event-based formulations [18] and precedence-based models [19] have been developed. Furthermore, a high interest has been expressed for real-life industrial study cases and problem specific solutions have been generated for real industrial facilities. Moreover, the ever-increasing computational power, allowed the handling of larger problem instances. However, there is still a significant gap between the academic research and the industrial practice, as only a few contributions have been successfully applied in real-life scheduling problems.

The rest of this paper is organized as follows. In Section 2, a detailed analysis of the theoretical concepts of optimization-based process scheduling is presented, including a classification of the different mathematical models, as well as a characterization of the problems they are able to address. Section 3 illustrates a systematic review on the application of optimization methods in real-life industrial scheduling problems in the process industries. The main modelling features and the industrial case study characteristics are summarized. In Section 4, we highlight the major challenges in applying optimization methods in real industrial problems and discuss potential remedies to close the existing gap between theoretical advancements and practical implementations. Finally, Section 5, draws up the main concluding remarks of this work

#### **2. Theoretical Aspects of Optimization-Based Process Scheduling**

As noted by Gabow [20], all scheduling problems are NP-hard, meaning that no known solution algorithms exist that are of polynomial complexity in the problem size. Therefore, the development of efficient optimization-based solution strategies for production scheduling has been a great challenge to the research community. As a result, a significant contribution emerged in the last decades aiming to develop either tailored algorithms for specific problem instances or efficient generic methods.

#### *2.1. Classification of Scheduling Problems*

Main goal of all scheduling problems is to propose a schedule that reaches the production targets, while respecting all operational, logistical and technical constraints, and achieves a certain objective, such as the maximization of profit, the minimization of the total cost, earliness and/or tardiness, and production makespan.

The general scheduling problem seeks to optimally answer the following questions (Figure 1):


Note that depending on the specifics of the problem in hand, some of these decisions are not considered in the scheduling level. When developing a model for the optimal scheduling all characteristics of the production must be considered to ensure the feasibility of the proposed schedules. However, the production needs to be portrayed in an abstract way to reduce the computational complexity of the problem. This is even more crucial when dealing with real-life industrial applications, which are typically characterized by complex structures, ever-expanding product portfolios and a huge number of constraints that must be considered.

**Figure 1.** Decisions of production scheduling in the process industries.

Traditionally, scheduling problems are defined in terms of a triplet α/β/γ [21]. The α field describes the production environment, while the β field denotes the special characteristics and production constraints. Finally, field γ describes the problem's objective, e.g., minimization of cost. The entries of this triplet can be extremely diverse between process industries, since a great variety of aspects needs to be considered when developing optimization models for process scheduling. As a result, many classes of scheduling problems exist. However, the general production scheduling problem can be summarized as follows:


The first term denotes the characteristics of the facility and can be considered static input to the scheduling problem, since it remains the same for all problem instances of a facility, unless any redesign studies are considered. The remaining terms are inputs from other decision-making processes in the manufacturing environment. Scheduling is not a standalone problem; it is part of the manufacturing supply chain and has strong connections to other planning functions. Production targets and materials availability come from the planning level, while resource availability is an output of the control level, thus there is a significant flow of information from other planning functions to scheduling (Figure 2).

**Figure 2.** Information flow towards scheduling level.

Scheduling is a critical decision-making process in all process industries, from the chemical and pharmaceutical to the food and beverage and the petrochemical sector. Besides the aforementioned general description of scheduling, industrial applications display strong differences to each other, due to the facility itself, the production policy or market and business considerations. First step when approaching an industrial scheduling problem is to identify its problem specifics, in order to accurately portray the problem in hand. Moreover, a strong correlation between different classes of scheduling problems and the available mathematical modelling frameworks exist. The scheduling problems found in process industries are classified in terms of: (a) The production facility, (b) the interaction with the rest of the production supply chain, and (c) the specific processing characteristics and constraints. A short description of these terms follows, the interested reader can find details in the excellent reviews of Maravelias [22] and Harjunkoski et al. [23].

#### 2.1.1. The Production Facility

At this point we should note that the following analysis focuses on production scheduling. However, many scheduling problems in the process industries target to the optimization of material transfer operations rather than production operations. Characteristic examples are crude oil and pipeline scheduling. With this in mind, the production facility is classified based on the type of process (batch/continuous) and the production environment (sequential or network).

#### Process Type

The type of production processes found in the process industries can be defined as continuous or batch. In the continuous mode, units are continuously fed and yield constant flow. Continuous processes are appropriate for mass production of similar products, since they can achieve consistency of product quality, while manufacturing costs are reduced, due to economies of scale. The main characteristic of batch processes is that all components are completed at a unit before they continue to the next one. Batch production is advantageous for production of low-volume high-added value products, or for production of seasonal demands which are difficult to forecast. One of the main advantages of batch production is the reduced initial capital investment, therefore, it is especially profitable for small business or trial runs of new facilities. From a scheduling point-of-view, both batch and continuous processes require the same type of decisions. Tasks can be characterized as batches or lots. Assignment (batches/lots to units), sequencing (between batches/lots) and timing (of batches/lots) decisions are identical, while selection and sizing of tasks (batching/lot-sizing) display more degrees of freedom in continuous processes. Capacity restrictions in continuous processes refer to processing rates and processing times and are usually unrestricted, thus a given order can be satisfied in a single lot (campaign) or multiple shorter ones. On the other hand, batch production is capacitated by the amount of processed material that a unit can process, thus affecting the number and size of batches to be scheduled. Another difference lies in the way inventory levels are affected. At this point, it is worth mentioning that many facilities are characterized by more than one type of processes. A characteristic example is the "make-and-pack" type of production, where several batch or continuous processing stages are followed by a packing (continuous) stage. This production flow is very common in the food and beverage and the consumer goods industries and requires the consideration of both the characteristics of batch and continuous production processes [24,25].

#### Production Environment

Production facilities can be classified as sequential or network based on the material handling restrictions. In sequential processing, each batch/lot follows a sequence of stages based on a specific recipe. Throughout its recipe a batch retains its identity, since it cannot be mixed with other batches or split into multiple downstream batches. Network facilities are characterized as more general and complex and have usually an arbitrary topology. Moreover, no restrictions exist for the handling of input and output materials, thus mixing and splitting operations are included. Based on their topological characteristics, sequential facilities can be further categorized into the following:

• Single stage: Production facility that consists of just one processing stage, which may consist of a single unit or multiple parallel units. The product to unit compatibility may be fixed (batch can be processed in a single unit) or flexible (batch can be processed in multiple units), but in all cases each batch must be processed in a single unit.

• Multistage: Each batch must be processed in more than one processing stages, each consisting of a single unit or multiple parallel units. The multistage environment can be further categorized into multiproduct and multipurpose, depending on the imposed routing restrictions. Multiproduct facilities are equivalent to flowshop environments in discrete manufacturing, where all products go through the same sequence of processing stages. In contrast, a facility is characterized as multipurpose when the routings are product-specific, or when a processing unit belongs to different processing stages depending on the product, thus being equivalent to jobshop environments in discrete manufacturing.

Early studies mainly focused on scheduling problems that are characterized as sequential [26,27]. Process industries with a sequential environment are very similar to discrete manufacturing, from a scheduling point-of-view. Sequential facilities can be easily modelled in terms of batches and production stages, like jobs and operations in discrete manufacturing. However, this does not hold true for network facilities, thus they cannot be modelled in a similar straightforward manner. Kondili et al. [1] followed by Pantelides [11] were the first to propose general representations of network facilities (STN, RTN), allowing the development of optimization models for scheduling problems of such complex structures. A classification of the production environments for process industries is illustrated in Figure 3.

**Figure 3.** Categorization of scheduling problems based on the production environment.

#### 2.1.2. Interaction with Other Planning Functions

Scheduling is strongly interconnected to the rest of the planning functions of the manufacturing supply chain; therefore, it cannot be approached as a standalone problem. The interactions between scheduling and the other decision making processes in a manufacturing environment must be accounted for, since they determine significant aspects of the scheduling problem; in particular: (a) the input parameters of the scheduling problem, (b) the decisions to be optimized by the scheduler, (c) the type of scheduling problem to be solved and (d) the problem's objective.

Planning and scheduling are two interdependent, however, distinct decision-making processes. Their differences lie in the level of detail of the used models, the time horizon and the problem's objective. In contrast to production scheduling, aggregate models are usually employed in planning, in order to specify the required produced amounts and storage levels that are able to satisfy a given demand at the minimum cost. Moreover, the planning horizon is much larger as it spans from weeks to months. The solution of planning determines the input of the scheduling problem in terms of production targets like order sizes, due dates and release dates. Additionally, batching/lot-sizing decisions can be made in the planning level, thus affecting the type of decisions that needs to be made in the scheduling level. In that case batching/lot-sizing decisions are pre-fixed and the scheduling decisions are narrowed down to just unit to task assignment, sequencing and timing of tasks. There is

also an important flow of information between scheduling and control; more specifically, the optimized schedule provides the reference points to the control level while resource availability is in turn provided to the scheduling level. Most studies until the early 2000s, approach production scheduling as a standalone problem. However, the scientific community acknowledged the importance of integrating the decision-making process of the various functions (planning, scheduling and control) that comprise the supply chain of a process industry [28]. The integrated planning and scheduling problem has been studied in multiple works in the last decades [29,30] and also implemented in industrial case studies with great success [31]. In contrast the integrated scheduling and control and integrated planning, scheduling and control problems have been only recently examined [32,33].

The demand volume and variability defined by the market environment in which an enterprise operates plays a pivotal role, since it specifies the type of the scheduling problem to be solved. On the one hand, high-volume production with relative constant demand based on forecasting favors a "make-to-stock" production policy, while the low-volume production with irregular demand follows a "make-to-order" policy. In the former the generated schedule is repeated periodically ("cyclic scheduling"), while in the latter a short-term schedule must be frequently generated. The choice of a meaningful objective for any production scheduling problem is a challenging task due to the numerous competing goals. The production characteristic that usually imposes the objective function is the relation between the capacity of the plant and the demand to be satisfied. In particular, when the demand overcomes the capacity of the plant, then objectives such as, the minimization of backlogs or the maximization of throughput are favored. On the contrary, if the capacity is enough to satisfy the demand, then the minimization of total cost is usually preferred as the overarching production goal. However, the definition of the production scheduling objective also strongly depends on market considerations and goals originating from other planning functions. For example, the maximization of throughput cannot be a valid objective in a production that must be fixed to the amounts defined in the planning level. It must be also noted that production scheduling is an inherently dynamic process, so the objective can be adjusted at any time due to market-related reasons, e.g., new or changed contracts or fluctuations in the demand.

#### 2.1.3. Processing Characteristics and Constraints

Scheduling problems may refer to facilities that exhibit various special processing characteristics and constraints. These aspects complicate the problem but must be considered, in order to ensure the feasibility of the generated production schedules. In the next section we will shortly review some of them and further details can be found in [34].

Resource considerations, aside from task-unit assignments and task-task sequences, are of great importance. These may involve auxiliary units (e.g., storage vessels), utilities (e.g., steam and water) and manpower. Resources are mainly classified into renewable (recover their capacity after being used in a task, e.g., labor) and non-renewable (their capacity is not recovered after being consumed by a task, e.g., raw materials). Renewable resources can be further classified into discrete (e.g., manpower) and continuous (e.g., electricity, cooling water). Another important characteristic in process industries is the handling of storage, which is usually referred to as the storage policy. Depending on the duration a material can be stored, the storage policies are described as (i) unlimited intermediate storage (UIS), (ii) non-intermediate storage (NIS), (iii) finite intermediate storage (FIS) and (iv) zero wait (ZW). Setups are a critical factor in most processing facilities as they represent operations like re-tooling of equipment, cleaning or transitions between steady states. They are associated with a specific downtime that can be sequence-independent or sequence-dependent (changeovers) and a cost is induced to the production process. To reduce the complexity associated with the consideration of setups, products are categorized into families. In that case setups exist only between products of different families.

This classification illustrates the complexity of scheduling problems and the tremendous diversity of aspects that must be accounted for when dealing with real industrial applications (Figure 4). The inherent diversification of scheduling problems in the process industries hindered the initial efforts

of the academic community to propose a unified general mathematical framework. Therefore, research turned into the development of less general methods that can address industrial cases that share similar characteristics. As a result, a multitude of efficient specialized methods for the optimization of scheduling in the process industries have been proposed in the last 30 years.

**Figure 4.** Information extracted from problem characteristics.

#### *2.2. Classification of Modelling Approaches*

As mentioned in the previous subsection, scheduling problems in the process industries are defined by extremely diverse features (e.g., production environment, processing characteristics etc.), while different aspects need to be taken into account based on external parameters, like the market environment in which the industry under study operates. Therefore, the initial attempts of proposing a mathematical framework that would constitute a panacea to all scheduling problems, were unsuccessful and soon solutions that take advantage of the problem-specific characteristics emerged. The struggle to overcome the computational complexity associated with scheduling problems, gave rise to numerous scheduling models. It should be noted that in this work we focus on optimization-based approaches, more specifically, the models presented are mixed-integer programming (MIP) models. Nevertheless, we should mention that an abundance of alternative solution approaches, e.g., constraint programming models [35,36], heuristics [37] and metaheuristics [38], exist in the literature. These methods can provide fast and feasible solutions, thus being a very attractive option for industrial case studies. However, their superiority in terms of computational complexity comes with a cost, since optimality of the generated schedules is not ensured. To combine the advantages of both optimization and non-optimization approaches, hybrid methods have emerged that are able to provide near-optimal solutions in low computational time [39].

The three main aspects that describe all optimization models for scheduling are: (i) The optimization decisions to be made, (ii) the modelling elements and (iii) the representation of time (Figure 5).

**Figure 5.** Main aspects of models for optimal production scheduling.

#### 2.2.1. Optimization Decisions

The optimization decisions are affected by the handling of batches/lots. As we underlined in Section 2.1.2, batching decisions may be optimized in the planning level, thus be prefixed and be an input to the scheduling problem. Even if this is not the case, the scheduler has the flexibility to decide whether the batching decisions will be part of the optimization model. For example, the decision-maker can heuristically specify the number and size of batches and then utilize an optimization approach for the unit allocation, sequencing and timing decisions. Usually models for sequential environments favor this two-step approach. In contrast, a monolithic approach, consisting of batching/lot-sizing, unit assignment, sequencing and timing decisions, is used for network environments. Few recent works have proposed a monolithic approach to deal with scheduling problems in sequential environments [40–42]. In some special cases, like in the single machine problems, only sequencing and timing decisions are optimized, thus reducing the scheduling problem to a traditional travelling salesman problem.

#### 2.2.2. Modelling Elements

According to the entity used to ensure the resource constraints on processing units, modelling approaches are classified into two categories: Batch-based and material-based. In sequential environments, where the identity of each batch remains the same throughout the processing stages, batch-based approaches are used. On the contrary a material-based approach is favoured, when dealing with network environments, where batches are mixed or split. It is important to mention that the modelling elements used are tied to the optimization decisions. More specifically, in monolithic approaches the scheduling problems are modelled using a material-based approach, while a batch-based approach is followed, whenever the batching decisions are known a priori.

The modelling elements are strongly tied with the representation of the manufacturing process, which is the core of every scheduling model. The goal of a successful representation is to translate the real problem (orders, units, stages) into mathematical entities (variables, constraints) in an abstract way, that will allow for the fast generation of optimal and feasible schedules. Even a simple manufacturing process may consist of multiple operations, therefore, the use of a simplified representation is essential. The oldest type of manufacturing process representation is utilized to model scheduling problems of sequential production environment and is based on (i) processing stages, (ii) processing units in each stage and (iii) batches or products (depending on whether batching decisions are prefixed or not). The second type of representation emerged in the early 1990s from the novel works of Kondili [1] and Pantelides [11], who introduced the STN and RTN, both based on the modelling of materials, tasks, units and utilities. The STN represents manufacturing processes as a collection of material states (feeds, intermediate final products) that are consumed or produced by tasks. The main difference between STN and RTN is that in the latter states, units and utilities are represented uniformly as resources that are produced and consumed by tasks. While originally introduced for

scheduling problems in network environments, recent works have addressed problems in sequential environments [43,44] using the RTN representation.

#### 2.2.3. Time Representations

The most studied topic and the one that mostly differentiates optimization models for scheduling is the representation of time. Depending on the way sequencing and timing of tasks are considered, modelling approaches are categorized in two broad approaches, in particular precedence-based and time-grid-based. Based on their type, precedence-based models are classified into general, immediate and unit-specific general precedence models and time-grid-based into discrete and continuous. Continuous-time formulation may employ single or multiple-time grids. Figure 6 illustrates the various time representation approaches in optimization models for scheduling.

**Figure 6.** Categorization of modelling approaches based on time representation.

All precedence-based models consist of unit-task allocation and task-task sequencing constraints [45]. The latter are expressed as precedence relationships between tasks processed in the same unit, while the former ensure that each batch/lot is processed by exactly one unit in each stage. Binary sequencing variables are introduced to enforce the precedence relationships and ensure the generation of a feasible schedule (no processing of multiple tasks simultaneously in the same unit). Another main characteristic of any precedence model is that the timing variables are not mapped onto an external time reference, rather their exact values are specified within the scheduling horizon based on the interactions (timing constraints) between pairs of batches/lots or between processing stages of the same batch. Two types of precedence variables exist: (i) General, where precedence relationships are established between all pairs of batches/lots and (ii) immediate, where they are established only between consecutive pairs. General precedence models require fewer variables, so they are more computationally efficient. However, these models do not identify subsequent tasks, making it difficult to consider changeover costs and heuristics, such as pre-fixing or forbidding certain processing sequences. To overcome this limitation, Kopanos et al. [39] proposed the unit-specific general precedence approach that combines both general and immediate sequencing variables. In all cases precedence-based models can provide high quality solutions with low computational cost, thus being an attractive alternative when dealing with real-life industrial problems. One of the main disadvantages of this approach is the quadratic increase of the size of the model with the number of batches/products considered. The use of information such as product families or pre-fixing of sequences mitigates this phenomenon and vastly improves the efficiency of the models [46].

Time-grid-based models enforce timing and sequencing constraints through the utilization of a single or multiple time grids, onto which events (e.g., starting or completion of task) are mapped. A great variety of time-grid-based approaches exist depending on the representation of events (time slots, global periods, time points or events), which are classified into discrete and continuous. In discrete-time models the time-grid is portioned into a pre-fixed number of global time periods of a known duration, both of which need to be specified by the modeler. Most discrete formulations use a common time frame for all shared resources. However, Velez and Maravelias [47] proposed a discrete model that employs multiple time frames. One of the main challenges when setting up discrete models is the proper

selection of the number of time periods that needs to be employed. A fine grid results to solutions of higher quality but in cost of larger less computationally efficient models. An advantage of discrete-time models is their capability of monitoring inventory and backlog levels, material balances, as well as the availability and consumption of utilities without introducing nonlinearities. Moreover, time-dependent utility-pricing and holding and backlog costs can be linearly modelled, while integration with higher planning levels is straightforward [48]. Additionally, discrete-time formulations are superior to their continuous counterparts in terms of solution quality [49]. Nevertheless, discrete formulations result in very large, however tight, models, especially when small discretization of time is mandatory. In continuous models, the horizon is subdivided into a fixed number of periods of variable length, which is defined as part of the optimization procedure. Single, common and multiple, unit-specific time frames have all been successfully employed to continuous-time models. Continuous formulations can alleviate some of the computational issues associated with discrete-time models, since fewer time periods, thus variables, are required for the representation of the same scheduling problem. However, they are not necessarily more computationally efficient compared to their discrete counterparts. Finally, it should be mentioned, that few models that utilize multiple ways of representing time have been proposed, thus combining both the advantages of discrete- and continuous-time formulations [29,50].

#### *2.3. Alternative MILP Models for Process Scheduling*

We already illustrated a classification of the various scheduling problems as well as the main modelling approaches that have been suggested in the last 30 years. A scheduling model is determined by both externally specified (problem class) and user selected (modelling approach) factors. On the one hand, the model should be suitable for the examined problem environment and the processing specifics of the facility under study, and on the other, it should be developed in terms of the chosen modelling approach's characteristics. A given problem can be represented in multiple ways, however there is a significant relationship between these two aspects. In this subsection we will demonstrate the basic aspects of the mathematical models that have been proposed by the scientific community. More specifically, we present an overview of the models based on the problems they are used for and we analyse the basic constraints and variables of representative models. Further details on the different mathematical models for production scheduling can be found in the excellent review of Méndez et al. [34].

#### 2.3.1. Models for Network Production Environments

In network environments batches do not maintain their identity, since mixing and splitting of batches is allowed. Therefore, the problem is presented utilizing either the STN or the RTN process representation (batch-based approaches). Moreover, the complexity of the production arrangement, with tasks consuming or producing multiple materials and materials being processed in different tasks and units, requires the proper monitoring of material balances, status of units and utility and inventory levels. This necessitates the utilization of a time-grid based approach.

A plethora of modelling formulation emerged after the introduction of the discrete STN and RTN models. Reklaitis and Mockus [51] were the first to propose a continuous-time STN formulation. A single common grid is used, in which the timing of the grid points ("event orders") was determined by the optimization. The model is an MINLP, which can be further simplified to a mixed integer bilinear problem that is solved using an outer-approximation algorithm. Zhang and Sargent [52,53] developed an RTN-based continuous time formulation that can address both batch and continuous operations. The ensued MINLP model is solved using a local linearization procedure in combination with a column generation algorithm.

One of the major drawbacks of the first models developed according to the continuous STN and RTN mathematical frameworks was the large integrality gap. This deficiency was addressed by Schilling and Pantelides [12,54]. They modified the formulation of Zhang and Sargent [53], simplifying it and improving its general solution characteristics, while they developed a hybrid branch-and-bound

solution method which branches in the space of the interval durations as well as in the space of the integer variables.

Castro et al. [17] proposed a relaxation of Schilling's formulation [12], allowing tasks to last longer than the actual processing time. Consequently, their model is less degenerate and less CPU time is required. Some of the co-authors further improved this formulation in [55], allowing the optimization of continuous processes. A novel common-grid STN-continuous formulation was introduced by Giannelos and Georgiadis [56]. They utilized a non-uniform time grid, that eliminates any unnecessary time events, thus leading to small MILP models. Maravelias and Grossmann [57] suggested a general continuous STN-model that accounts for various processing characteristics such as different storage policies, shared storage, changeover times and variable batch sizes. The contribution of Sundaramoorthy and Karimi [58] is another well-known continuous MILP model that introduced the idea of several balances (resource, time, masses etc.).

The concept of multiple unit-specific time grids was first proposed by Ierapetritou and Floudas [18]. This approach decouples the task events from the unit events, thus less slots are required. As a result, smaller MILP models are generated, leading to a significant decrease in computational effort. Multiple works have been proposed ever since, improving the computational characteristics and expanding the scope of the initial formulation [59–61].

Velez and Maravelias [47] were the first to introduce the concept of multiple, non-uniform discrete time grids. The multiple grids can be unit-, task- and material-specific. The same authors extended this work in [62] with the consideration of general resources and characteristics like changeovers and intermediate storages. It should be noted that while these formulations were initially proposed for network facilities, they can be also used for the scheduling of sequential environments.

We will now focus our attention on two representative scheduling models for network environments. First, we will consider the continuous common-grid model by Castro et al. [55]. Here an RTN representation is employed, while the model utilizes a common grid to express the timing constraints. More specifically, a set of global time points *T* is predefined throughout the scheduling horizon. The major decisions are expressed through the binary allocation variable *Ni*,*t*,*t* that is enabled whenever a task starts at time point *t* and is completed at or before point *t* - . The rest of the decision variables are continuous and express the exact time that corresponds to each time point *Tt*, the size of a batch/lot of a task ξ*i*,*<sup>t</sup>* and the amount of resource being consumed at each time point *Rr*,*t*. The major constraints of the model can be summarized as follows:

$$T\_{t\prime} - T\_t \ge \sum\_{i \in I\_r} (\alpha\_i \cdot \overline{N}\_{i,t,t\prime} + \beta\_i \cdot \overline{\xi}\_{i,t,t\prime}), \qquad \forall r \in \mathbb{R}^l, t \in T, \ t\prime \in T, \ t \le t\prime \tag{1}$$

$$T\_{t\prime} - T\_t \le H \cdot \left(1 - \sum\_{i \in I\_t} \overline{N}\_{i, t, t\prime} \right) + \sum\_{i \in I\_t} \left(a\_i \cdot \overline{N}\_{i, t, t\prime} + \beta\_i \cdot \overline{\xi}\_{i, t, t\prime} \right)\_{\prime} \qquad \forall r \in R^I, \ t \in T, \ t\prime \in T, \ t \le t\prime \tag{2}$$

$$\mathcal{V}\_{i}^{\min} \cdot \overline{\mathcal{N}}\_{i,t,\mathsf{t}\prime} \leq \overline{\mathcal{I}}\_{i,t,\mathsf{t}\prime} \leq \mathcal{V}\_{i}^{\max} \cdot \overline{\mathcal{N}}\_{i,t,\mathsf{t}\prime} \qquad \forall i \in I, \mathsf{t} \in T, \mathsf{t}\prime \in T, \mathsf{t} \leq \mathsf{t}\prime \tag{3}$$

$$\mathcal{V}\_{i}^{\text{min}} \cdot \overline{\mathcal{N}}\_{i,t+1} \le \sum\_{r \in \mathbb{R}\_{i}^{ST}} \mathcal{R}\_{r,t} \le \mathcal{V}\_{i}^{\text{max}} \cdot \overline{\mathcal{N}}\_{i,t,t+1} \quad \forall i \in I, t \in T \tag{4}$$

$$V\_i^{\min} \cdot \overline{N}\_{i, t - 1, t} \le \sum\_{r \in \mathbb{R}\_i^{ST}} R\_{r, t} \le V\_i^{\max} \cdot \overline{N}\_{i, t - 1, t\_r} \quad \forall i \in I, t \in T \tag{5}$$

$$R\_r^{\text{min}} \le R\_{r,t} \le R\_r^{\text{max}} \quad \forall r \in R, \ t \in T \tag{6}$$

$$R\_{r,t} = R\_{r,t-1} + \sum\_{i \in I\_{r}} \left[ \sum\_{t \le t} \left( \mu\_{r,i}^{p} \cdot \overline{N}\_{i,t,t} + \boldsymbol{\nu}\_{r,i}^{p} \cdot \underline{\boldsymbol{\varsigma}}\_{i,t,t} \right) + \sum\_{t \ge t} \mu\_{r,i}^{c} \cdot \overline{N}\_{i,t,t} + \boldsymbol{\nu}\_{r,i}^{c} \cdot \overline{\boldsymbol{\varsigma}}\_{i,t,t} \right] + \tag{7}$$
 
$$\sum\_{i \in I\_{ST}} \left( \mu\_{r,i}^{p} \cdot \overline{N}\_{i,t-1,t} + \mu\_{r,i}^{c} \cdot \overline{N}\_{i,t,t+1} \right), \quad \forall r \in R, \ t \in T, t > 1$$

Assuming that no more than one task can be executed in each unit at a certain time (unary resource), constraint sets (1) and (2) guarantee that the time difference between any pair of time points *t* and *t*' must be at least equal to the processing time of all tasks starting and finishing at those points. Furthermore, the batch/lot size is bounded by the unit capacity (3), while constraints (4) and (5) impose the storage constraints. They ensure that in case of a resource excess at time t, the corresponding storage task has to take place for both *t* − 1 and *t*. Finally, constraint (7) guarantees that the resource balance considerations are not violated.

Next, we present the model of Janak et al. [60]. In contrast to the previous model, an STN-representation is chosen. Tasks are mapped onto multiple time grids through the concept of event points. These are time instances located along the time axis of each unit that represent the starting of a task. Due to the incorporation of a unit-specific grid, fewer time points are required compared to common-grid formulations, thus the number of binary variables is significantly reduced. The main variables of the model are *Wi*,*t*, *W<sup>s</sup> <sup>i</sup>*,*<sup>t</sup>* and *<sup>W</sup><sup>f</sup> <sup>i</sup>*,*<sup>t</sup>* denoting that a task *i* is active, started or finished at event point *t*, accordingly. This formulation is one of the most general of the ones that employ a unit-specific grid, since it has the ability to account for various storage policies, batch splitting and mixing, changeovers and variable batch sizes. As a result, it involves a huge number of constraints, hence only the major ones will be presented here.

$$\sum\_{i \in I\_j} W\_{i,t} \le 1, \qquad \forall j \in J, t \in T \tag{8}$$

$$\mathcal{W}\_{i,t} = \sum\_{t\nu \le t} \mathcal{W}\_{i,t\nu}^{\text{st}} - \sum\_{t\nu \le t} \mathcal{W}\_{i,t\nu'}^{f} \qquad \forall i \in I, t \in T \tag{9}$$

$$\sum\_{t \in T} \mathcal{W}\_{i,t}^{s} = \sum\_{t \in T} \mathcal{W}\_{i,t'}^{f} \qquad \forall i \in I \tag{10}$$

$$B\_{i,j,t}^{s} \le B\_{i,j,t}, \qquad \forall i \in I, j \in J\_i, t \in T \tag{11}$$

$$B\_{i,j,t}^{s} \le V\_i^{\max} \cdot \mathcal{W}\_{i,t'}^{s} \qquad \forall i \in I, j \in J\_{i\prime}, t \in T \tag{12}$$

$$B\_{i,j,t}^{s} \ge B\_{i,j,t} - V\_i^{\max} \cdot \left(1 - W\_{i,t}^{s}\right)\_{\prime} \qquad \forall i \in I\_{\prime}, j \in J\_{i\prime}, t \in T \tag{13}$$

$$S\_{s,t} = S\_{s,t-1} + \sum\_{i \in I\_s^P} \rho\_{i,s} \cdot B\_{i,j,t-1}^f + \sum\_{i^{ST} \in I\_s^{ST}} B\_{i^{ST},t-1}^{ST} - \sum\_{i \in I\_s^C} \rho\_{i,s} \cdot B\_{i,j,t}^f + \sum\_{i^{ST} \in I\_s^{ST}} B\_{i^{ST},t'}^{ST} \qquad \forall s \in S, t \in T, t > 1 \tag{14}$$

$$T\_{i,j,t}^f \ge T\_{i,j,t}^s \qquad \forall i \in I, j \in J\_i, t \in T \tag{15}$$

$$T\_{i,j,t}^f \le T\_{i,j,t}^c + H \cdot \mathcal{W}\_{i,t}, \qquad \forall i \in I, j \in J\_{i\prime}, t \in T \tag{16}$$

$$T\_{i,j,t}^s \ge T\_{i\nu, j, t-1}^f + H \cdot (1 - W\_{i\nu, t-1}), \qquad \forall j \in \mathcal{J}, i \in \mathcal{I}\_j, i\nu \in \mathcal{I}\_j, i \ne i\nu, t \in T, t > 1 \tag{17}$$

$$\begin{aligned} T\_{i,j,t}^{\mathbb{S}} \ge T\_{i\nu,j\iota,t-1}^{f} + H \cdot \left(1 - W\_{i\nu,t-1}^{f}\right)\_{\prime} \quad & \forall s \in \mathbb{S}, i \in I\_{s\prime}^{c}, i\nu \in I\_{s\prime}^{p}, j \in \mathbb{J}, \\ & j \in I\_{i\nu}; j\nu \in I\_{i\nu}, j \neq j\nu, t \in T, t > 1 \end{aligned} \tag{18}$$

The major assignment constraints (8)–(10) impose that: (i) At most, one task can be executed by unit *j* at time *t* (unary resource), (ii) the assignment variable *Wi*,*<sup>t</sup>* will be active only if the task *i* has started but not finished at or before time *t* and (iii) each task *i* must start and finish within the given scheduling horizon. Batch-size considerations are employed by constraints (11)–(13). In particular, they bound the amount of material starting processing at time *t, Bs <sup>i</sup>*,*j*,*<sup>t</sup>* according to the unit capacity and relate it to the amount undertaking task *i* in unit *j* at time *t, Bi*,*j*,*t*. Constraint set (14) enforces the material balances, stating that the amount of state *s* at time *t*, *Ss*,*t*, is increased by the amount produced and stored at time *t* − 1 and decreased by the amount consumed or stored at time *t*. The timing constraints (15) and (16) relate the starting *T<sup>s</sup> <sup>i</sup>*,*j*,*<sup>t</sup>* and completion *<sup>T</sup><sup>f</sup> <sup>i</sup>*,*j*,*<sup>t</sup>* times of a task *i* in unit *j* and at time *t*. More specifically, they impose that the completion time must be larger or equal to the starting time, and that if the task *i* is not processed in unit *j* at time *t*, the completion time is set equal to the starting time. Constraint (17) ensures that if task i finishes at time *t* − 1 and task *i* starts at time *t* in the same unit, then task *i* must start after the completion of task *i* - . Finally, let us consider a task *i* that produces a state *s* at time point *t* − 1 that is used by task *i* in time *t*. To respect the production recipe task *i* must start after the completion of task *i* - . This sequencing consideration is enforced by constraint (18).

#### 2.3.2. Models for Sequential Production Environments

Scheduling problems of sequential environments do not share the same complexity, in terms of problem representation, with the ones encountered in network environments. Therefore, both precedence-based and time-grid based approaches can be employed. Each of these approaches display specific advantages and drawbacks. On the one hand, precedence-based models generate smaller, more intuitive models that provide high quality solutions, on the other hand time-grid based models are usually tighter and computationally superior. As a result, a great variety of models have been proposed to address sequential production environments.

One of the most impactful time-grid based models is [63] from Pinto and Grossmann. They described an MILP model for the minimization of earliness of orders for a multiproduct plant with multiple equipment items at each stage. The interesting feature of the model is the representation of time, where two types of individual time grids are used: One for units and one for orders. Castro and Grossmann [64] proposed a non-uniform time grid representation for the scheduling problem of multistage multiproduct plants. They tested their formulation for various objectives, e.g., minimization of makespan, total cost and total earliness and compared it with other known formulations, concluding that the efficiency of a model highly depends on the objective and the problem characteristics. The same authors extended their work in [43] with the consideration of sequence-dependent setup times.

Unlike to most of the other contributions, which propose continuous-time models, the work of Maravelias and co-workers thoroughly investigated the employment of discrete-time models in sequential environments. Sundaramoorthy et al. [65] suggested a discrete time model to incorporate utility constraints for the scheduling problem of multistage batch processes. Merchan et al. [66] developed four novel formulations, two of them based on the STN and RTN representation and two more inspired by the resource-constrained project scheduling problem (RCPSP). Moreover, the authors introduced tightening constraints and reformulations that allowed for significant computational enhancements. Recently, Lee and Maravelias [67] presented two new MIP models for scheduling in multipurpose environments using network representations. Interestingly, states and tasks were defined based on batches instead of materials, making possible the consideration of material handling constraints in sequential production environments. The authors displayed the potential of the proposed models by incorporating important process features, such as time-varying data and limited shared resources, and by solving medium-size problem instances to optimality.

The concept of precedence has been extensively studied by the PSE community [68–70]. Numerous unit-specific immediate [71], immediate [72] and general precedence [19,73] models have been proposed for scheduling problems in sequential environments. In initial studies the batches to be scheduled was a problem data, however later contributions suggested models for the simultaneous batching and scheduling problem [74].

Let us consider the general scheduling problem of a multistage multiproduct facility with multiple units operating in parallel in each stage. Moreover, we assume that the batching decisions are fixed and provided to the scheduler from the planning decision level. This problem can be efficiently tackled by numerous precedence-based models. Here we use the formulation proposed by Méndez et al. [19] and present its core constraints. The main decision variable of all precedence-based models is a Boolean indicating the sequential relation between any pair of orders. More specifically, in the presented formulation *Xo*,*o*-,*<sup>l</sup>* defines whether an order *o* is processed prior to order *o* at stage *s*. Other characteristic decision variables are the binary allocation variable *Yo*,*j*, defining whether an order *o* is executed by unit *j* or not, and *Co*,*<sup>l</sup>* that denotes the completion of order *o* in each stage. The main constraints of the model are illustrated below:

$$\sum\_{j \in \mathbb{J}\_{sl}} Y\_{o,j} = 1, \qquad \forall o \in O, l \in L \tag{19}$$

$$\mathbb{C}\_{o,l} \ge \sum\_{j \in I\_{o,l}} \left[ \mathbb{Y}\_{o,j} \cdot \left[ pt\_{o,j} + su\_{o,j} \right] \right], \qquad \forall o \in O, l \in L \tag{20}$$

$$\begin{split} \mathbb{C}\_{o\nu,l} - pt\_{o\nu j} \ge & \mathbb{C}\_{o,l} + su\_{o\nu,j} + \tau\_{o\rho\nu,j} - M \cdot \left(1 - \mathbf{X}\_{o\rho\nu,l}\right) - M \cdot \left(2 - \mathbf{Y}\_{o,j} - \mathbf{Y}\_{o\nu,j}\right)\_{l} \\ & \forall o, o\nu \in \mathcal{O}\_{\boldsymbol{\nu}} o\boldsymbol{\nu} > o, l \in \mathcal{L}\_{o,\boldsymbol{\mu}\boldsymbol{\nu}}, j \in \mathcal{I}\_{o,\boldsymbol{\mu}\boldsymbol{\nu}} \end{split} \tag{21}$$

$$\begin{aligned} \mathsf{C}\_{o,l} - pt\_{o\prime j} \ge \mathsf{C}\_{o\prime,l} + su\_{o,j} + \tau\_{o\prime,j,j} - M \cdot \mathsf{X}\_{o,\nu\prime,l} - M \cdot \left(2 - \mathsf{Y}\_{o,j} - \mathsf{Y}\_{o\prime,j}\right)\_{\mathsf{V}} \\ \forall o, o \in \mathsf{O}, o\prime > o, l \in \mathsf{L}\_{o,\nu\prime}, j \in \mathsf{I}\_{o,\nu\prime,l} \end{aligned} \tag{22}$$

$$C\_{o,l} \le C\_{o,l+1} - \sum\_{j \in I\_{o,l+1}} (pt\_{o,j} \cdot Y\_{o,j}) \tag{23}$$

Constraint (19) ensures that each order *o* is processed by exactly one unit *j* in each stage *l*. The main timing considerations are specified by constraint (20), which enforces the completion time of an order *o* executed by unit *j* to be at least equal to the required processing and setup time. Big-M parameters are employed to express the sequencing constraints (21) and (22), between any pair of orders *o* in each stage *l.* A major difference to immediate precedence models, is that here only one sequencing variable is defined for every pair of orders, as a result the size of the model is significantly reduced. Moreover, both constraints become active only when both orders are processed by the same unit, i.e., *Yo*,*<sup>j</sup>* = *Yo*-,*<sup>j</sup>* = 1, therefore the unit index is omitted from the precedence variables. If order *o* is processed before order *o* in the same unit, constraint (21) becomes active, ensuring that order *o*- will be completed after the completion of order *o* plus the required processing time of *o* and any sequence-dependent or -independent setup times, while constraint (22) becomes redundant. In the opposite case where order *o* is processed earlier than *o*', constraint (22) is activated and (21) becomes redundant. Finally, constraint (23) guarantees the correct sequence between processing stages for the same order.

At this point we should emphasize that no modelling approach exists that is computationally superior to the others in every type of scheduling problem. While discrete-time approaches generate tighter models, their continuous-time counterparts (precedence-based, continuous time-grid-based) require less variables, thus generating smaller-sized models. Extensive comparative studies on scheduling problems in sequential environments conclude that time-grid-based models tend to be generally superior to precedence-based ones [43,64]. However, we must note that the computational efficiency of a model can drastically change even with small alterations in the facility characteristics and the final objective. This will be more evident in our analysis in Section 3, which accentuates the case-specific nature of the problem. Finally, consider that most modelling developments have been tested in small or medium sized study cases, that usually do not represent real-life industrial scheduling problems. Consequently, the computational efficient of any optimization-based model itself is not sufficient enough to address large-sized industrial problems. Thus, as we present in the following section, the introduction of techniques, such as heuristics and decomposition algorithms, is inevitable.

#### **3. Real-Life Process Systems Industrial Applications**

As described in the previous section, a plethora of different mathematical models has been proposed for tackling the production scheduling problem. Except from solving literature problem examples, several researchers, mainly from the PSE community, expressed a high interest for handling real-life industrial case studies. Numerous modelling approaches and methods can be found in the open literature, addressing a great variety of industrial process scheduling problems. A categorization based on the industrial sectors, such as chemical, pharmaceutical, petrochemical, steel, food and consumer goods industries, is presented below, along with the proposed modeling approaches. We focus our attention on MILP-based approaches for the offline scheduling problem, excluding other solution methods (e.g., heuristic rules, metaheuristic algorithms etc.).

#### *3.1. Chemical Industries*

One of the main industrial sectors widely studied, considers chemical plants, where a variety of new products is produced via the chemical transformation of multiple raw materials. The use of mixed batch and continuous processes, the special equipment technologies and the necessity to achieve a specific quality of products are the main challenges in these problems. In chemical plants, various types of products can be manufactured via the same or a similar sequence of operations by sharing the several plant's production units, intermediate materials, and other production resources. Lin and Floudas [75] proposed a continuous time, event-based MILP scheduling model and a decomposition methodology, to solve large-scale industrial cases of multiproduct batch plants. A real-life study case of a chemical plant, including 3 stages, 35 final products and 10 pieces of equipment is considered. To systematically apply the proposed approach, a graphical user interface is developed. Depending on each problem instance, the computational time of the proposed approach ranges from 15 min to 7 h. Janak et al. [76] extended the previous approach, by adapting intermediate due dates and other technical constraints. A unit specific, event-based formulation is applied in parallel with a decomposition-based approach, utilizing the rolling horizon technique. Problem instances with up to 67 product orders have been considered and a termination criterion of 3 h CPU time has been used. Westerlund et al. [77] introduced a mixed discrete-continuous time formulation to tackle short-term and periodic scheduling problems of multi-product plants, including intermediate storage constraints. As the suggested approach is focused on industrial applications, good quality solutions are targeted in reasonable computational times instead of global optimal solutions. The mixed discrete-continuous model provides better solutions in smaller computational times, in comparison with the discrete-time approach. A strategic planning tool was developed based on the proposed model and applied to an industrial plant, importing demand data from the plant's ERP (Enterprise Resource Planning) system. Additionally, four scheduling approaches have been developed by Velez et al. [78]. Here the idea of multiple discrete-time grids is utilized, as each material, task and unit has its own time grid. Upper bounds on the total production of each material are defined using the concept of the effective time window for the executed tasks. Further extensions are adapted in order to solve a variety of different problems. The introduced methods have been applied to benchmark problem instances that can be found in literature [79] and to a real case study from Dow company [80], including five main product lines. Near optimal solutions are achieved in 1 h CPU time on average. A comparison of the proposed approaches and four other continuous time formulations has been also presented. The results indicate that the discrete time models generate better solutions in less computational time.

#### *3.2. Pharmaceutical Industries*

A special subsector of the chemical plants is the pharmaceutical industry. The majority of the operations taking place in these facilities are batch, as there is a high necessity to ensure the quality of the final products. Moniz et al. [81] motivated by a real-world scheduling problem of a chemical-pharmaceutical industry, developed a case-specific discrete-time MILP scheduling model, for batch plants. All the data used by the mathematical formulation, is taken automatically from a decision-making tool and a process representation, developed as a prototype in Microsoft Visio. A representative industrial case, including four products, nine shared processing units and 40 tasks, has been studied. The solutions can be generated in acceptable computational times according to the plant operators and even for larger problem instances, suboptimal but good quality solutions are provided in 1 h CPU time. Stefansson et al. [82] studied a large-scale industrial case study from a pharmaceutical company, including even 73 products and 35 product families. Mathematical frameworks based both on discrete and continuous time representations have been proposed and a comparison of them is illustrated. The initial problem is decomposed into two subproblems and the stage which constitutes the main production bottleneck is scheduled first. The continuous-time formulation can provide better solutions even for larger problem instances. Case studies with up to 400 orders can be solved by utilizing the continuous time formulation and schedules with 9.8% integrality gap are generated in 1408 min. Optimal schedules for smaller case studies, involving up to 150 product orders, can be generated in less than 1 h. On the other hand, only up to 75 products can be scheduled to optimality by utilizing the discrete time formulation, as suboptimal solutions with 10% integrality gap are generated for instances with up to 300 products. Castro et al. [83], presented a decomposition-based algorithm for tackling the high complexity of large-scale problems of multiproduct facilities. The production orders are inserted iteratively into the generated schedule, allowing some flexibility to provide better solutions. A case study comprising of 50 production orders, 17 units and six stages is efficiently solved in less than 1 min. The same pharmaceutical study case has been also considered by Kopanos et al. [39]. They proposed a decomposition-based solution strategy relying on two precedence-based MILP models in order to optimize different objectives, such as makespan, changeover-time and cost minimization. A feasible schedule is rapidly generated, and it is enhanced by applying an improvement algorithm. High quality solutions are provided for industrial cases with up to 60 products allocated to 17 units. Liu et al. [84] focused on the production and maintenance planning problem of biopharmaceutical process, consist of a fermentation and a purification stage. Maintenance activities related to the regeneration of the column resin, taken place in the purification stage, are considered. Two industrial indicative problem instances are illustrated to assess the applicability of the proposed MILP model and global optimal solutions are found, without exceeding the time limitation of 1 h CPU time. An event-based continuous time mathematical framework based on the STN representation has also been proposed for a general multiperiod biopharmaceutical scheduling problem [85]. Optimal solutions can be found in computational time in the range of 1–2 min.

#### *3.3. Petrochemical Industries*

A special interest is expressed for the scheduling problem of oil refineries or petroleum industries. A variety of products are produced by this specific industrial sector, such as gasoline, diesel jet fuel and others. Many different and complex processes are taken place in the oil refineries; therefore, their efficient scheduling constitutes a great challenge. Shah et al. [86] motivated by a study case provided by Honeywell Process Solutions (HPS), considered an MILP based heuristic algorithm. The initial oil refinery problem is spatially decomposed into two subproblems, one considering the production and blending and the other the delivery of the finished products. Feasible solutions are generated by solving the two subproblems iteratively, via a six-step heuristic algorithm as the resolution of the direct proposed MILP model is characterized by a high computational cost. Ten different problem instances were presented, for the production of diesel and jet fuel and nearly optimal solutions were generated within less computational time. In particular, the computational time ranged from 2 s to 1 h depending on the cases' complexity. Zhang and Hua [87] deployed a plant-wide multi-period planning model, aiming to the integration of the plant processes and the utility system, in order to reduce the energy consumption. The plant-wide model is extended by considering the utility system model and constraints referring to the utilities' balances such as steam, fuel oil and gas are adapted. The maximization of the total profit of the whole refinery plant is considered as the objective. The optimal operation modes of units and stream flow are defined by the model. Product blending and maintenance activities are taken into account. As the process system and the utility system are optimized separately by the suggested hierarchical method and the subproblems' complexity decreased, good quality, but not global optimal solutions are generated in acceptable computational times. The applicability of the approach was illustrated in a real study case that considers a refinery industry, located in South China. The refinery, except from importing of electricity to cover its power needs, was also able to export the surplus power back to the network or other power companies. The integrated problem of investment planning and operation scheduling of offshore oil facilities was also addressed, by utilizing a multiperiod MILP model in order to maximize the general profit [88]. Various operational nonlinear constraints related to the reservoir performance and to other resources are efficiently adapted into the proposed model which was solved by utilizing a decomposition algorithm in order to handle the high complexity. A real-life, large-scale illustrative example was considered. Although a feasible solution can't be returned by solving the exact MILP model, a feasible solution within 6600 s was obtained, by utilizing the proposed decomposition algorithm. The operation scheduling of a crude oil terminal has been considered by Assis et al. [89]. A real-life case study, oriented by the national refinery of Uruguay was considered and near optimal solutions were obtained by using an iterative two step MILP-NLP algorithm within a time limit of 3600 s. A domain reduction relaxation was also adapted for handling the emerging bilinear terms.

Other than the processes taking place in the refinery industries, a special interest has been also expressed from the PSE community, in the scheduling of liquid transportation via pipeline systems in petroleum supply chain. The crude oil is gathered and transported to the refineries, as the final refined products are sent to the retail market and distributed to customers. In order to reduce the transportation time, pipelines are preferred instead of using trucks or other means of transport, providing also more safety and lower CO2 emissions. Castro and Mostafaei [90] motivated by the scheduling problem of liquid transportation, proposed an event-point MILP formulation for treelike transportation systems, where a single input node leads to multiple outputs. A continuous time representation was utilized and novel constraints for ensuring the avoidance of forbidden product sequences were adapted. A real-life study case from the Iranian Oil Pipelines and Telecommunication Company network was considered and the optimal schedules could lead to even a 6.2% capacity increase, as the given demand can be efficiently covered fourteen hours earlier. A comparison with previous methods, proposed from one of the co-authors [91], indicates the efficiency of the approach. A time termination criterion of 5 h has been used for the proposed formulation. The number of the event points has been identified as key parameter with high impact on the computational time and solution quality. Nearly optimal solutions can be generated in less than 1 h by reducing the available event points. A similar problem, referring to the scheduling of a transportation system of petroleum products, produced from a single oil refinery industry was tackled by Cafaro and Cerdá [92]. They proposed an MILP continuous time model in order to define the optimal lot size, the batch sequence, as well as the delivery time of batch order. A variety of constraints were taken into account, such as tank availability and quality control operations. A real-life study case consisting of six different oil derivatives produced by a unique oil refinery to a single distribution center was scheduled, and the results indicated that better solutions were produced in comparison with other approaches for the same problem, in less than 60 s CPU time. The same problem was also addressed by Cafaro et al. [93], but now allowing simultaneous product deliveries, thus providing more realistic solutions. The proposed two-level MILP-based solution technique aimed for the minimization of the total number of operations in order to reduce the number of restarts and stoppages of the pipeline. On the upper level, the feasibility of the problem was ensured, as more detailed decisions such as lot sizing, lot sequencing and timing decisions were defined on the lower level. A study case related to REPLAN refinery industry, consisting of five distribution centers at Brazil, is used to illustrate the applicability of the model. Significant savings were noticed in CPU time using the multiple delivery policy, as the illustrative examples under consideration can be solved in less than 125 s CPU time. Rejowski and Pinto [94], inspired also from the REPLAN refinery,

proposed two discrete-time, MILP-based models, to solve a real-life problem, including the distribution of various petroleum products to five depots. An indicative instance of a 75-h time horizon that is discretized in 5-h intervals was presented. A good quality solution with integrality gap of 5.8% was returned within the time limit of 10,000 CPU s. Boschetto et al. [95], proposed an MILP-based solution algorithm for solving a large scale real-life pipeline network problem, by determining the delivery and the pumping times of 14 different oil products and ethanol, to a number of distribution centers. Efficient heuristic rules were utilized in parallel with a continuous time representation, for tackling the daily scheduling problem, in reasonable computational times within 3–5 min. The generated solution for various studied cases, have been also validated by the planners.

#### *3.4. Food Industries*

The PSE community has also shown significant interest for the scheduling of food industries. Common characteristics of food processing industrial facilities, such as intermediate due dates, shelf life considerations and multiple mixed batch and continuous processing stages, substantially complicate the optimization of scheduling decisions. The above combined with market trends that enforce the gradual increase of the product portfolio, the demand profile (high variability-low volumes), and the multiple identical machines and shared resources, make the consideration of real-life industrial cases extremely challenging.

As the food industry focuses mainly on the production of perishable final products a make-to-stock production policy is not efficient, as the generation of high inventory levels should be avoided. A plethora of industrial study cases have been considered from various subsectors of the food industry. Baldo et al. [96] motivated by a real study case from a Portuguese brewery industry, proposed a novel MILP-based relax and fix heuristic algorithm, for the integrated fermentation and packing problem. The time horizon is discretized in two subperiods. The first subperiod is scheduled in detail, as for the second subperiod only the main planning decisions, such as the inventory levels, are optimized. Small and big sized problem instances have been considered, with five filling lines and up to 40 products. Although a direct comparison with the company plan was not possible, good quality schedules were generated, using a termination criterion of runtime limit equal to 3600 s or 7200 s. An immediate precedence-based MILP formulation for the packing stage of a brewery company was developed using a mixed discrete-continuous time representation in [29]. The scheduling decisions were defined in a continuous manner, while material balances were expressed at each discrete time period to ensure the generation of feasible schedules. The idea of grouping the products into product families leads to significant reduction of the computational cost. Changeover times among sequential time periods were also taken into account. The industrial study case under consideration consists of eight processing units and 162 products are produced in total, which are grouped into 22 product families. The generated solutions were better than the ones extracted by commercial tools. An upper bound of 300 s CPU time was utilized, for all cases under consideration. Abakarov and Simpson [97] investigated the scheduling problem of food canneries focusing on the sterilization stage and allowing the possibility of the simultaneous sterilization of different products in the same retort. A graphic user interface, able to identify the nondominated simultaneous sterilization vectors, was connected to the proposed MILP model. Different cases were solved, including 16 products with randomly generated product demand values, depicting a reduction of up to 25% in total plant operation time. The usage of COIN-OR as software tool can decrease the model's computational time to 7.38 s. Georgiadis et al. [98] studied the integrated sterilization and packing stage scheduling problem in a large-scale canned fish Spanish industry. An MILP based decomposition algorithm was utilized to tackle the high computational cost, as the products are inserted in an iterative way until the final schedule was generated. A general precedence model efficiently describes the batch (sterilization) and the continuous (packing) processes of the plant. Nearly optimal schedules of a large-scale problem instance, with 100 final products and 362 product batches, have been generated for both stages, in less than 20 min. A study case of a real-world edible-oil deodorized industry was studied by Liu et al. [99]. The plant was described

as a single-stage multiproduct batch process. The final products were grouped into product families having the same due date. The proposed approaches relied on mixed discrete and continuous MILP mathematical formulations and classic TSP (travelling salesman problem) constraints. A real study case of 128 hours' time horizon of interest was studied. 70 product orders of 30 different final products of seven groups of different release time were scheduled. The new formulations are shown to be more efficient than previously proposed methods found in the literature. Solutions with approximately 2% integrality gap can be generated in 20 CPU s without allowing the backlog generation and 1075 CPU s by allowing the possibility of backlogs. Polon et al. [100] studied a sausage production industry aiming to the profit maximization by solving an MILP scheduling model for batch processes. The packaging stage, which often constitutes the main production bottleneck has not been considered. The plant operates in a single campaign mode and eight products are produced in total.

A special subsector of food industries is dairy manufacturing. Numerous products are produced, such as yoghurt, cheese and butter and distributed to customers worldwide. Doganis and Sarimveis [101] solved the scheduling problem of a single yoghurt production line taking into account inventory, manpower and capacity restrictions. The model was tested using data from a yoghurt production line of a Greek dairy industry, where 18 products were produced and global optimum schedules have been generated in less than 15 s. The integrated planning and scheduling problem of a small size Balkan type semi-continuous yoghurt facility, with eight final product types, produced by three intermediates has also been investigated [102]. The evaluation of the proposed MILP approach has been utilized via a simulation model. Thirty-two different scenarios were considered and a significant decrease in the total waste and makespan was achieved in approximately 1 h of CPU time. Touil et al. [103] deployed an MILP model for a small multiproduct milk industry, located in Morocco, aiming at the minimization of makespan. The stages of homogenization, pasteurization and packaging were scheduled for four final products, seven packing lines, two pasteurization units and one homogenizer. Efficient solutions were illustrated for the cases under consideration, as optimal schedules can be found in 2 min CPU time. A novel mixed discrete-continuous MILP formulation was deployed by Kopanos et al. [104] for the scheduling problem of a Greek yoghurt production facility. The idea of "product families" was adapted similarly to the other aforementioned works from the same authors. The packing stage was scheduled in detail, but mass balance constraints related to the production stage were also adapted, using a discrete time representation. Ninety-three final products (grouped into 23 product families) were allocated in four packing lines. Novel resource constraints can adapt realistic limitations to various types of resources (e.g., manpower) and ensure the generation of feasible optimal solutions in less than 10 min, depending on the case complexity. Based on a similar approach, the scheduling problem of another large-scale Greek dairy industry has been studied [105]. A rolling horizon technique was embedded to reactively adjust the schedule in case of disturbances, like the cancellation or modification of orders, or the sudden arrival of new orders or any digressions from the planned production. One hundred and fifty-eight final products (grouped into 44 product families) were allocated to six parallel packing lines, while the time horizon of interest was five days. A total cost decrease of 20% was achieved in comparison with the schedules generated by the company. An integrated software tool with a user-friendly graphical interface has been developed to connect the proposed MILP model to the input data, located in excel files (parameter values such as changeover times etc.) and the ERP system (providing the demand values). As a result, optimal solutions can be generated automatically in less than 10 min.

#### *3.5. Consumer Goods Industries*

Consumer goods, or final goods, are described as products consumed by the average customer. Depending on the shelf-life duration, they can be further categorized to durable goods (such as detergents) and nondurable (e.g., beverages). One of the main consumer goods group is the fast-moving consumer goods (FMCG), which are characterized by frequent purchases, rapid consumption and low prices. Elzakker et al. [106], presented a problem-specific model for the short-term scheduling

problem, considering a fast-moving consumer goods (FMCG) industry. An algorithm based on a unit-specific, continuous time interval MILP model is proposed. Dedicated time intervals to specific product types are adapted to decrease the computational time. In order to assess the efficiency and the applicability of the proposed formulation ten industrial case studies are considered, as provided by Unilever, related to an ice cream production process. Optimal schedules have been generated for problem instances of up to 73 batches of eight products allocated to six storage tanks and two packing lines within 170 s. The time-horizon under consideration was 120 h. The production scheduling problem of an ice cream facility has also been tackled by Kopanos et al. [107]. A real-life study case of eight final ice cream products, two packing lines and six aging vessels is addressed. The simultaneous optimization of all processing stages is achieved, and 50 problem instances are optimally solved. An MILP-based decomposition strategy is proposed to handle scheduling problems of large-scale food process industries. High quality solutions were generated for larger cases of up to 24 final products utilizing the proposed decomposition technique. Industry related needs imposed the adaptation of a 600 CPU s as a time limit, but global optimal solutions can be found in less than 10 s for smaller problem instances.

Giannelos and Georgiadis [108] developed an MILP model to address the scheduling problem in fast consumer goods manufacturing processes. The proposed MILP model relied on the STN process representation and a continuous time formulation was used to reduce the computational complexity of the problem. The formulation was tested on a medium-sized industrial consumer goods manufacturing process, considering cases with up to 35 final products and 5 packing lines. Feasible schedules were generated within a 5–10% integrality gap in computational times, smaller than 5 min. Méndez and Cerdá [109] proposed a general precedence MILP formulation based on a continuous time representation, while they introduced constraints related to sequence dependent setup times and products' due dates. Furthermore, efficient preordering rules were considered in order to provide solutions for industrial study cases of up to 18 final products, produced from five intermediates over a scheduling period of five days. The CPU time needed was gradually depleted to less than 10 s for small problem instances and even to 10 min to larger ones, by applying the proposed approaches and the suggested preordering rules. Baumann and Trautmann [110] proposed a hybrid method for large-scale, short-term scheduling problems of packed consumer goods products. A subset of the operations were scheduled iteratively by solving a general precedence MILP model [24]. Also, an iterative improvement step is applied to the initial schedule, by following a reinsertion policy identifying some critical operations. Ten large-sized instances provided by The Procter and Gamble Company that consisted of up to 1391 operations have been solved within reasonable CPU times of less than 1 h, as a 5 s time limit has been set for each iteration [111]. Smaller-scale problem instances with known optimal solutions, have also been considered, and optimal or near-optimal schedules were generated by applying the aforementioned hybrid method. Elekidis et al. [112] investigated the short-term scheduling problem of a large-scale consumer's goods industry. An immediate-general precedence-based model was illustrated, focusing mainly on the packing stage. Constraints related to the previous stages were also taken into account. The production orders are inserted iteratively, utilizing a decomposition algorithm. Various real-life study cases have been considered that include up to six packing lines and 130 final products. Near optimal schedules are generated and significant savings in the changeover time are noticed within a CPU time of 10 min.

Georgiadis et al. [113] presented two different scheduling approaches, based on the RTN and the STN representations respectively. The work was focused on the scheduling problem of large-scale manufacturing industries of electrical appliances. A case study provided from a large manufacturing company located in Greece was used to assess the applicability of the proposed approaches. The generated schedules can be visualized via the Microsoft Excel application. A significant decrease in the operational cost was reported in a variety of problem instances. Although, the necessary computational time was in the range of some seconds, this could differ as the considered problem instances are described as data-driven.

#### *3.6. Steel Plants*

Another important field of interest is the steel-making process industries. Various challenges arise, due to the large variety of final products, the complex process that take place and the volatile electricity prices. The steel production is often divided into three stages: Molten steel is produced first (melt shop) and then the produced slabs are transformed (the hot rolling) into intermediate or final-products, (e.g., coils, billets etc.). In the last stage (cold casting), the dimensions and the desired mechanical properties are achieved. Biondi et al. [114], studied the scheduling problem of a hot rolling mill in a steel plant. Strict production constraints related to metallurgic production were taken into account. Decisions regarding the production planning on a hot rolling mill were taken, by applying intelligent heuristics, resulting to efficient production programs. A slot-based MILP formulation was proposed and the sequence of the aforementioned programs was defined. The technique described above, has also been implemented as a web-service, that gathers information from both the ERP system and the DCS (Distributed Control System). Rolling programs including up to 3000–5000 coil orders can be generated within some seconds. Yang et al. [115], proposed an MILP mathematical formulation to tackle the scheduling problem by optimizing the byproduct gas systems in steel plants. Optimal solutions can be found in approximately 10 s. A representative case study from a steel plant in China has been considered to illustrate the proposed approach. A significant reduction of even 7.8% in the operation cost was noticed. Li et al. [116] considered the scheduling problem of steel making industries, focusing mainly on the steelmaking continuous casting process, as it constitutes the main production bottleneck. A novel unit-specific event-based continuous-time MILP model was proposed, relied on material continuity and other technological requirements constraints in order to ensure the generation of feasible schedules. An extension of previous rolling horizon approaches [75,76] is also applied due to the high complexity of the large-scale problems under consideration. Four representative industrial problems have been considered, to assess the efficiency of the proposed approach. Although, not even a feasible solution could be returned by solving the MILP model for the two larger problem instances within the time limit of 80,000 s, good quality solutions were generated by using the new proposed approach in 3 s and 12,287 s respectively.

Gajic et al. [117] studied the integrated scheduling and electricity optimization problem of a hot rolling mill, taking also into account electricity costs and prices. An MILP-based model was proposed in parallel with intelligent heuristics, aiming to group the individual heats into casting sequences and decompose the large problem into several sub-problems with lower complexity. The proposed approach has been successfully implemented via offline tests and it has been deployed in the melt shop at Acciai Speciali Terni S.p.A., a member of ThyssenKrupp AG and one of the world's leading producers of stainless steel based in Italy. The scheduling solutions could be generated within a few minutes and it has been shown that the electricity costs can be reduced by 3% as the coordination among the different production stages was significantly improved. Hadera et al. [118] proposed a new general precedence MILP scheduling model adapting energy awareness. Optimal production schedules were generated, simultaneously optimizing the electricity purchase and solving the load commitment problem. The case of selling electricity back to the grid was also taken into account. In order to handle large-scale industrial problem instances of a melt shop section of a stainless-steel plant, a bi-level heuristic algorithm was used. Solutions in the range of 9% and 25% integrality gap were obtained within the time limit of 600 s.

The scheduling problem of multiproduct plants with parallel units, implementing energy intensive tasks has also been considered. Continuous and discrete time RTN-based mathematical formulations have been proposed and tested to a few study cases from a real industrial problem [119]. Efficient solutions were generated and optimality gaps of 1% were achieved within 5 min of computational time. Significant savings of electricity costs, even up to 20%, were reported [120]. The same problem was also considered by Kong et al. [121]. The authors proposed an MILP model, targeting the minimization of the operational cost by optimizing the by-product gas distribution. The proposed MILP model has been successfully tested in a real-life case study, provided by an iron and steel plant located in

China. The results show that the operational cost was reduced by up to 6.2% by applying the proposed approach. Wang et al. [122] investigated the bi-objective single machine batch scheduling problem of a real-world scheduling problem in a glass company located in Shanghai, China. An exact ε-constraint method was adapted to the MILP model in order to minimize the makespan and the total energy costs. Two heuristic methods were proposed to tackle the high complexity of larger scale problems. A representative real-life case study, including 13 batches has been studied and a pareto curve was generated, illustrating the tradeoffs between the two objectives. Not even a feasible solution can be obtained after 31,823.64 CPU s by utilizing the direct ε-constraint method. However, approximate Pareto fronts, including 11 different solution points, can be generated in 5648.05 and 6359.25 CPU s, by using the two proposed heuristic methods respectively.

#### *3.7. Paper Industries*

A special interest has been expressed for the problem of trim loss minimization, mainly in the paper industry. In case the final products are to be divided in sizes of specific dimensions, significant trim losses are unavoidably generated, leading to important increase in the operational cost. Westerlund et al. [123] studied the trim-loss problem of a Finnish paper-converting mill. A two-step optimization procedure based on an MILP model was solved by CPLEX solver in fractions of a second, resulting to waste savings of 2% of the turnover. A sequential updating procedure [124] has been also presented. The researchers proposed that by resolving an MILP model iteratively, feasible schedules can be promptly extracted. A case with up to 10 products and three machines was tackled and the provided solutions are better than the ones generated manually. According to Roslöf et al. [125] various sophisticated heuristics can be utilized in large scale industrial problems to provide feasible suboptimal solutions in reasonable computational times. Extending the previous approach, an extra improvement reordering step was introduced that can lead to nearly optimal solutions. Various problem instances have been considered and the proposed approach has been compared with a heuristic based policy, and manually generated schedules. Better solutions can be provided by targeting at tardiness and makespan minimization utilizing the proposed method. A real-life case study provided by a Finnish paper mill included 61 scheduling jobs and a single processing unit was solved in 3755.1 s of CPU time. Giannelos and Georgiadis [126] proposed a slot-based MILP scheduling model, which relied on a continuous time representation, to examine the problem of cutting operations on parallel slitting machines. A CPU time limit of 1500 s has been adapted, and good quality solutions have been obtained within 6% to 9% integrality gap. The proposed approach has been applied to an industrial case study, provided by a paper mill company (Macedonian Paper Mills, S.A., Greece), including eight final products and consisting of three parallel cutting machines. Castro et al. [127] proposed an MILP and an MINLP mathematical model, which were based on a continuous and a discrete time RTN representation. The non-linearity of the second formulation can be eliminated assuming a constant throughput. The aforementioned frameworks were applied to an industrial case study from a pulp mill plant, located in Portugal. According to the detailed comparisons, the discrete time formulation seems to be more efficient, as optimal schedules can be generated in less than 600 s of CPU time, depending on the different level of discretization and the differentiation of the problem instances. On the other hand, the computational cost of the continuous time formulation seems prohibitively high, as by increasing the number of the event points by one, an increase of one order of magnitude appeared also to the computational time. Castro et al. [128] proposed an RTN-based formulation and showcased a detailed comparison between continuous and discrete time models by applying them in an industrial problem consisting of three raw materials, five intermediates and five product qualities. Novel recycling policies were also adapted and as a result, significant reduction of raw materials can be achieved, providing higher profits and lower waste. On the contrary, with the previous research, the proposed continuous time formulation led to a better quality and faster solution than the discrete one. A similar pattern also appeared in the computational time, as an increase of one order of magnitude appeared by increasing the number of the event points by one.

Table 1 presents the contributions reviewed in this section, illustrating the industrial sector for which the optimization method has been developed, as well as their main features.


**Table 1.** Summary of the industrial applications using optimization-based scheduling approaches.


#### **Table 1.** *Cont*.

#### **4. Industrial Applications of Optimization-Based Scheduling—Challenges**

In the previous section, a wide range of real-life applications using various scheduling frameworks have been described. It can be noticed that most of the methods have efficiently handled small or medium sized problem instances, with just a few of them addressing large scale industrial problems.

A major issue, referred to the applicability of the scheduling approaches, is the accuracy of the input data. In many cases the generated schedules are manually modified by the planners in a time-consuming process, as some important information is missing [111]. If the data given by the plant managers are not accurate, the assessment of the model's efficiency is extremely difficult or even impossible. We can conclude that the connection of the proposed solution method with the plants' ERP and the other plant systems, plays a key role for a successful model implementation. Integrated software should be developed, to provide an easy way to transport the necessary input data to the model solver in an automated way. In addition, the direct communication between the

mathematical model and the ERP systems, could also make possible the consideration of the several planning decisions of the plant (e.g., planned cleaning or maintenance activities).

The aforementioned suggestion could also give the chance to solve and analyze a plethora of different cases, in order to test the efficiency of the model. All possible operational scenarios should be provided from the plant managers to the model developers. Further important production parameters could also be identified by checking the feasibility of the schedules generated by the model during the test runs.

The proposed models could efficiently solve all problem cases. The company should identify the largest or the most complex cases of the plant. Different models could also be addressed for different problem instances. Also, as the scheduling and planning decisions are strongly connected, there is a high necessity to ensure that the planning decisions are feasible. For example, the customer's demand, given usually by the ERP systems, must be covered by the plant's capacity during the available time horizon. Otherwise, not even a feasible solution could be returned from the model. However, even in cases when the optimization fails, an output should be exported [130].

The model's output should also be visualized in an interactive way. The planners should be able to modify the schedules by hand, or adding late order deliveries. We have to keep in mind that the planners are the ones responsible for the final schedule and an optimization model can provide a suggested solution to them. Flexible Gannt charts, allowing the right or left swift movement of the production orders, could also be an efficient tool for the planners. Also, in that way, a number of modifications, based on the planners' experience can be easily adapted to the final schedule [23].

Moreover, as the industrial environment is highly dynamic and it is characterized by a high level of uncertainty, a lot of daily modifications should be applied. Therefore, efficient rescheduling methods should also be provided, utilizing fast execution models. The utilization of decomposition algorithms or heuristic rules in parallel with the MILP models could also be a good quality option, as nearly optimal solutions can be generated in reasonable computational times [96,107]. It should also be noticed, that the updated schedule has to consider the previous schedule as an initial condition. In this way more realistic solutions can be provided.

As in the majority of the plants, production scheduling is done on a daily basis, the computational time of the utilized models should not be extremely high. Even good quality suboptimal solutions could be acceptable in time limitations defined by the company. Moreover, in practice due to unexpected events, the initial schedule rarely is completely applied.

An analysis of the expected savings could also be accomplished. The consideration of case studies of past weeks could provide a good estimation of the optimization savings range. The company should analyze the investment that has to be made to install a new integrated software and the possible time or cost savings of it. The necessary investment has to lead to higher profits. The company personnel should also be trained to use and modify efficiently the new optimization tools.

The proposed model approach should also be able to adapt new information or different types of modifications. New products could be inserted, by identifying their basic production features, as well as new pieces of equipment. Parameter values, such as, units' production rates, or new product allocation policies should be easily changed by the planners. The new integrated software should be flexible enough to adapt to the plant's changes, as otherwise the solution approach could be considered obsolete. The optimization tools should be easily used, and modified by the planners in a daily basis, without the supervision of experts in model development [131].

#### **5. Conclusions**

Prof. Sargent is one of the greatest personalities in the area of chemical engineering. He is the pioneer and father of the field of process system engineering. He was one of the first researchers who foresaw the value of the need for the optimization of process scheduling. His contributions inspired numerous researches in the last three decades, resulting in a plethora of mathematical programming formulations for general classes of process scheduling problems. This work presents a review of the theoretical modelling aspects and solution approaches for the process scheduling problem. Although there is no single comprehensive approach for handling all possible industrial problems, and the majority of the proposed models are characterized as problem-specific, a categorization is possible by identifying basic common features. Various problem instances, as well as different optimization models are classified, presenting their most important characteristics and the main contributions of them. Furthermore, an overview on the applications of the modeling formulations in real-life industrial case studies is presented. The industrial studies under consideration are categorized according to the different process industries subsectors, focusing on chemical, pharmaceutical, food, consumer goods, steel and paper industries. It is concluded that only a limited number of the contributions are able to solve large-scale industrial problems. In most of the large-scale problem instances, a variety of decomposition techniques and heuristic rules are applied in parallel with the mathematical programming models, and as a result; good but suboptimal solutions are obtained. A very small number of the generated solutions were directly compared with the plant's schedules. The development of integrated software tools, aiming to the direct data transfer between the mathematical models and the plant's ERP systems, according also and to the general trend of digitalization, was identified as a crucial step for the successful industrial applications of scheduling methods. Another scheduling challenge to be confronted is the efficient and flexible visualization of the generated solutions, allowing also for manual modifications.

**Author Contributions:** All authors contributed to this work. Individual contributions of the authors can be specified as: "Conceptualization M.C.G.; Analysis and Investigation G.P.G. and A.P.E.; Writing – Original draft G.P.G. and A.P.E.; Writing – Review & Editing M.C.G., G.P.G. and A.P.E.; Supervision M.C.G.".

**Funding:** This research was funded by the European Union's Horizon 2020 research and innovation program (SPIRE PPP framework) grant number 723575 (Project CoPro).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**


*Parameters*



*Variables*


#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Adjustable Robust Optimization for Planning Logistics Operations in Downstream Oil Networks**

#### **Camilo Lima 1, Susana Relvas 1,\*, Ana Barbosa-Póvoa <sup>1</sup> and Juan M. Morales <sup>2</sup>**


Received: 27 April 2019; Accepted: 19 July 2019; Published: 2 August 2019

**Abstract:** The oil industry operates in a very uncertain marketplace, where uncertain conditions can engender oil production fluctuations, order cancellation, transportation delays, etc. Uncertainty may arise from several sources and inexorably affect its management by interfering in the associated decision-making, increasing costs and decreasing margins. In this context, companies often must make fast and precise decisions based on inaccurate information about their operations. The development of mathematical programming techniques in order to manage oil networks under uncertainty is thus a very relevant and timely issue. This paper proposes an adjustable robust optimization approach for the optimization of the refined products distribution in a downstream oil network under uncertainty in market demands. Alternative optimization techniques are studied and employed to tackle this planning problem under uncertainty, which is also cast as a non-adjustable robust optimization problem and a stochastic programing problem. The proposed models are then employed to solve a real case study based on the Portuguese oil industry. The results show minor discrepancies in terms of network profitability and material flows between the three approaches, while the major differences are related to problem sizes and computational effort. Also, the adjustable model shows to be the most adequate one to handle the uncertain distribution problem, because it balances more satisfactorily solution quality, feasibility and computational performance.

**Keywords:** distribution; planning; oil supply chain; robust optimization; uncertainty

#### **1. Introduction**

In the oil industry, companies develop a set of activities and processes so that crude oil can be duly transformed into oil products demanded by final consumers such as petrochemical industries, airports and individual users [1]. The associated decision-making process is framed within a deeply complex context and addresses the oil exploration, production and transportation at upstream, oil refining at midstream, and oil product distribution and marketing at the downstream segments [2]. Most of these activities are difficult tasks and need to be properly accomplished. Decision-support tools based on mathematical programming techniques are often applied to assist such activities and processes across the oil supply chain [3]. These optimization tools are tailored to cover the underlying reality under study, and hence they are often stochastic large-scale programs, which normally increases the problem complexity [4].

In management sciences and engineering, the most widely-used modeling framework to cope with uncertainty has been two-stage stochastic programming [2], where the first-stage variables must be decided before the uncertainty realizations and the second-stage variables are decided in accordance with the uncertainty outcomes [5]. When either the information is missing or even the probability distribution of the uncertainty is not available, fuzzy programming could be an alternative framework to handle uncertainty in optimization problems by considering random parameters as fuzzy numbers

and treating constraints as fuzzy sets [6]. Robust optimization (RO) can be used, which also requires moderate information about uncertainty described by means of uncertainty sets, i.e., the key building blocks in this modeling framework [7]. Other frameworks have also been proposed to support decision making under uncertainty such as parametric programming, which relies on sensitivity analysis theory, and also chance-constrained programming, which uses probabilistic theory [8]. Currently, data-driven optimization methods are a new and timely research line to investigate the decision-making process under uncertainty, supported not only by the crescent access to uncertainty data, but also by the development in machine learning approaches [9]. Data-driven methods aims at building uncertainty sets directly from uncertainty [10]. In general, the choice of the mathematical modeling framework to solve a particular problem under uncertainty relies on the specific characteristics of this problem, as well as on the level of information about uncertainty.

The use of decision-support tools specifically for the DOSC problem under uncertainty has been a research challenge [11], and different approaches have been proposed in this direction. These include stochastic programming [1,2,4,8,12–21]; fuzzy programming [22,23]; and robust optimization [16,24]. In these optimization models, demand, supply, price, cost, process capacity, product yield, conversion rate and technology development are regarded as uncertain parameters [3].

Nowadays, RO has become attractive to solve optimization problems under uncertainty as it only requires the mean and range of the uncertain data, and other probabilistic information can be gradually incorporated when it becomes available [25]. In addition, it balances adequately solution quality, feasibility and computational performance [26], being formulated to hedge the optimization problem against any disturbances in the input data series [27]. As single-stage non-adjustable robust optimization (NARO) models may lead to over-conservative solutions [28–30], some authors have been emphasizing the need of exploring other methods such as the two-stage or multistage adjustable robust optimization (ARO) models, since less conservative solutions can be [28]. The first-stage and second-stage decisions are robust against all uncertain realizations of the random parameters, as well as the second-stage solution has full adaptability to the uncertainty [25].

Against this background, the purpose of this paper is to address the problem of refined products distribution under demand uncertainty in the downstream oil supply chain (DOSC). Extending previous work developed by the authors [1], where a stochastic programming approach was used considering uncertainty on price and demand for crude oil. Scenarios were built by using time series analysis. In the present work a robust modeling framework is developed that considers the problem as a two-stage adjustable robust optimization (ARO), aiming to conclude on the model conservativeness as well as on the associated solution performance. The framework includes a polyhedral uncertainty set, where the first-stage decision refers to purchase of oil volumes and the second-stage decisions include oil refining planning, oil product distribution planning, inventory management, international trade, and customer fulfillment. To evaluate the proposed model, a real-case study in the Portuguese oil industry is used. Results are compared with those obtained using the equivalent non-adjustable, stochastic and deterministic models. The robust model proposed in this paper is, to the best of the authors' knowledge, the first model to consider an adjustable robust framework to deal with uncertainty in DOSC and, therefore, an academic contribution for this area.

The rest of the paper is structured as follows: Section 2 describes the problem under study. In Section 3, the background on two-stage ARO modeling framework is explained. Section 4 proposes the ARO mathematical formulation. In Section 5, the case study based on the Portuguese oil industry is explored. In Section 6, the obtained numerical results are discussed. Finally, Section 7 draws some conclusions and presents some directions for future research.

#### **2. Problem Description**

The problem under study consists of determining the robust tactical and operational planning of a downstream oil network under uncertainty over a planning time horizon. It must be modeled satisfactorily in a robust optimization framework so as to maximize the worst-case DOSC profit. Profit is computed as the difference between the total cost of oil purchase associated with the first-stage decisions, and the worst-case operational profit associated with the second-stage recourse decisions, in order to satisfy the uncertain product demand, subject to material balance, refining capacity, product yield, supply availability, stock limitations and transportation capacity constraints.

The network under study is composed by three echelons: oil refineries, storage depots and local markets. Crude oil is acquired by oil refineries, which process it into oil products in order to satisfy the internal and external demands. The oil products are adequately stored across the network facilities, and thus conveyed by pipelines, tanker ships, tank wagons and tank trucks from oil refineries through the primary distribution, and by tank trucks from storage depots in the secondary distribution. Oil refineries are allowed to import and export oil products, while storage depots can only import. In summary, the problem can be generally defined as follows:

*Given* a downstream oil system, in which an oil company manages the refining processes, storage activities, logistics, distribution, and marketing;

*Determine* the robust tactical and operational planning under uncertainty so as to define the production levels, yield fractions, material flows, inventory levels, demand fulfillment, transport assignment and international trade throughout the network;

*Subject to* material balance, product yield, equipment capacities, supply limitations, and transportation capacity constraints;

*So as to* maximize the worst-case network profit that comprises the worst-case recourse profit minus the total cost of crude oil procurement.

To solve this problem, a two-stage adjustable robust modeling framework, considering uncertainty in demand for refined products, is developed.

#### **3. Two-Stage ARO Formulation**

In the two-stage ARO modeling framework, the first-stage variables comprise the set of decisions that shall be made before the uncertainty is revealed and the second-stage variables include the set of decisions that must be made after the uncertainty is disclosed. The second-stage decision variables have fully adaptability to any realization of the uncertainty [25]. Considering a linear formulation and let *<sup>x</sup>* and *<sup>y</sup>* be the corresponding first-stage and second-stage decision variables, while <sup>U</sup> is either a polyhedral or a discrete uncertainty set, according to Zeng and Zhao [28], the general form of a two-stage ARO formulation can be defined as follows:

$$\max\_{\mathbf{x}} \mathbf{c}^T \mathbf{x} + \min\_{\mathbf{u} \in \mathcal{U}} \max\_{\mathbf{y} \in F(\mathbf{x}, \mathbf{u})} \mathbf{b}^T \mathbf{y} \tag{1}$$

$$\text{s.t.}\\\text{A.x} \le d, \mathbf{x} \in \mathbb{S}\_{\mathbf{x}} \tag{2}$$

where *<sup>F</sup>*(*x*, *<sup>u</sup>*) <sup>=</sup> *<sup>y</sup>* <sup>∈</sup> *Sy* : *Gy* <sup>≤</sup> *<sup>h</sup>* <sup>−</sup> *Ex* <sup>−</sup> *Mu* with *Sx* <sup>∈</sup> <sup>R</sup>*<sup>n</sup>* <sup>+</sup> and *Sy* <sup>∈</sup> <sup>R</sup>*<sup>m</sup>* <sup>+</sup>. If the uncertainty set U is continuous, the optimization problem (Equations (1) and (2)) will have an infinite number of constraints in the second-stage problem [27]. However, when U is a polyhedral uncertainty set, only its extreme points, i.e., its vertices, might belong to an optimal solution of that second-stage problem [25]. In accordance with Zeng and Zhao [28], when considering only a finite and discrete set of extreme points <sup>U</sup> = {*u*1, ... , *ur*}, as well as the set of the respective recourse decision variables *y*1, ... , *yr* , the prior optimization problem can be recast as presented below:

$$\max\_{\mathbf{x}} \mathbf{c}^T \mathbf{x} + \eta \tag{3}$$

$$\text{s.t.}\\A\mathfrak{x} \le d$$

$$\eta \le \mathbf{b}^T \mathbf{y}\_{v^\nu} v = 1, \dots, \tag{5}$$

$$\mathbf{E}\mathbf{x} + \mathbf{G}y\_{\upsilon} \le \mathbf{h} - \mathbf{M}u\_{\upsilon}, \upsilon = 1, \dots, r \tag{6}$$

$$\mathbf{x} \in \mathbb{S}\_{\mathbf{x}}, \boldsymbol{y}\_{\upsilon} \in \mathbb{S}\_{\mathbf{y}}, \upsilon = 1, \ldots, r. \tag{7}$$

consequently, we can find a solution for the two-stage ARO model (Equations (1) and (2)) through solving its mixed-integer linear program counterpart (Equations (3)–(7)). However, defining this equivalent formulation and determining its solution can often be difficult and dependent on the size of the finite and discrete uncertainty set U = {*u*1, ... , *ur*} [28].

#### **4. Adjustable Robust Mathematical Programming Model**

Having defined in a generic form the two-stage ARO approach, the tactical and operational planning for the DOSC is now considered, which aims at maximizing the worst-case network profit under refined product uncertainty. In the established two-stage ARO problem, the first-stage decisions, i.e., the *here-and-now* decisions, comprise the crude oil purchase, while the second-stage decisions, i.e., the recourse (*wait-and-see*) decisions, include production levels, product yields, material flows, inventory levels, demand fulfillment, transport allocation and international trade, which are made in accordance with the realization of the product demand uncertainty. The two-stage ARO model for the general distribution problem under uncertain demand for refined products is formulated hereinafter.

#### *4.1. Robust Objective Function*

The notation used to formulate the two-stage ARO model is outlined in Table A1 of Appendix A. As shown in Equation (8), the objective function aims at maximizing the worst-case network profit (WCP), which can be stated as follows:

$$\max\_{XCO} \mathcal{WCP} = -\sum\_{i \in I, t \in T} \mathcal{XCO}\_{i,t} \mathcal{PO}\_{t} + \mathcal{Q}\_{\mathcal{U}}(\mathcal{XCO}) \tag{8}$$

$$\text{is.t.}\\XCO\_{i,t} \le DO\_{i,t} \\ \forall i \in I, t \in T \tag{9}$$

where *XCOi*,*<sup>t</sup>* represents the crude oil purchase at oil refinery *i* in time point *t* and refers to the first-stage decisions. QU(·) is defined as a function of the first-stage (*here-and-now*) decisions and denotes the worst-case recourse profit, which must be determined for the worst-case scenario of the uncertain product demands in the second-stage problem. Here, uncertainty is modeled through a polyhedral uncertainty set U, where the uncertain parameters can take values in, because such fact guides to define a robust counterpart that corresponds to a linear programming problem [27]. Equation (9) limits the amount of crude oil purchased by each oil refinery *i* at time point *t* by the corresponding oil demand.

The worst-case recourse profit QU(·) is equal to the objective function value of the next robust formulation, Equation (10), which can be formulated as follows:

$$\begin{split} \mathcal{Q}\_{\mathcal{U}}(\cdot) &= \min\_{\substack{\text{ADD}\_{k,p},\text{eff}\,\text{Rpt}\in\mathcal{T}}\,\text{I}} \max\_{\substack{\text{F}\in\text{I},\text{Lpt}\,\text{T} \\ j=1,p\in\text{P},t\in\text{T}}} \left( \sum\_{i\in\text{I},t\in\text{T}} \text{RR}\_{i,t} + \sum\_{j\in\text{I},p\in\text{P},t\in\text{T}} \text{MD}\_{j,p,t} + \sum\_{k\in\text{K},p\in\text{P},t\in\text{T}} \text{MC}\_{k,p,t} \right) \\ &- \sum\_{i\in\text{I},p\in\text{P},t\in\text{T}} \text{CE}\_{i,p,t} - \sum\_{h\in\text{H},p\in\text{P},t\in\text{T}} \text{CI}\_{h,p,t} - \sum\_{i\in\text{I},p\in\text{P},t\in\text{T}} \text{CP}\_{i,p,t} \\ &- \sum\_{j\in\text{I},p\in\text{P},t\in\text{T}} \text{CID}\_{j,p,t} - \sum\_{i\in\text{I},t\in\text{T}} \text{CID}\_{i,t} - \sum\_{j\in\text{I},p\in\text{P},t\in\text{T}} \text{CR}\_{i,p,t} \\ &- \sum\_{j\in\text{I},p\in\text{P},t\in\text{T}} \text{CID}\_{j,p,t} - \sum\_{k\in\text{K},p\in\text{P},t\in\text{T}} \text{CIM}\_{k,p,t} - \sum\_{k\in\text{K},p\in\text{P},t\in\text{T}} \text{CID}\_{k,p,t} \Big) \end{split} \tag{10}$$

where U represents the polyhedral uncertainty set that is used to describe the stochastic parameters, R ∈ <sup>R</sup><sup>+</sup> denotes the set of recourse decisions, i.e., <sup>R</sup> <sup>=</sup> {*XR,XOP,XCO,XP,XS,XRP,XU,XE,XI,IO,IR,ID,IM*}, H ∈ <sup>R</sup><sup>+</sup> depicts the set of *here-and-now* decisions, <sup>H</sup> <sup>=</sup> {*XCO*}, and <sup>F</sup> is a set defined as a function of the set of *here-and-now* decisions H and the worst-case realization of the deviations Δ*DPk*,*p*,*<sup>t</sup>* from the nominal values *DPk*,*p*,*t*. Notice that the set of recourse decisions R is made after the outer minimization

problem selects the worst-case realization of the deviations of the product demands Δ*DPk*,*p*,*t*. In other words, the *min*-*max* problem picks the worst realization of the demand that minimizes the network profit. Therefore, the robust modeling framework concerns with the maximization of the recourse profit in the worst-case scenario of the random parameter within the uncertainty set U. In summary, the two-stage RO model addresses the volumes of exported and imported products, refined oil, product yield fractions, primary and secondary product distributions, inventories of oil and products, and met and unmet demand as the second-stage decisions, which are adjusted to the volumes of purchased oil, i.e., the first-stage decisions, after uncertainty is disclosed.

The worst-case recourse profit QU(·), i.e., Equation (10), includes the revenues from the sales of refined products by the oil refineries, margins from buying and selling refined products by the storage depots and local markets, and costs for exporting, importing, transporting, and storing refined products throughout the network. Before presenting the equations of the recourse optimization problem QU(·), the associated revenues, margins and costs are now discussed, being outlined by the set of Equations (11)–(22) below:

$$RR\_{i,t} = \sum\_{p \in P} \left[ X R\_{i,p,t} \times \left( P P\_{a,p,t} - T N\_{r,p} \right) \right] \forall i \in I, t \in T, a = r == r \text{-}primary \tag{11}$$

$$\begin{split} MD\_{j,p,t} &= \sum\_{(k,m)\in Rautc\_{i,j,m}} XS\_{j,k,m,p,t} \times \left( PP\_{a,p,t} - TN\_{r,p} \right) \\ &- \sum\_{(i,m)\in Rautc\_{i,j,m}} XP\_{i,j,m,p,t} \times PP\_{a\_2,p,t} \\ &\dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \end{split}$$

$$\forall j \in I, p \in P, t \in T, a = r = \text{depot}, a\_2 = r \text{fiver} \\ y \tag{12}$$

$$\begin{aligned} \mathsf{MC}\_{k,p,t} &= \left[ \begin{pmatrix} \mathsf{X} \mathsf{R} P\_{k,p,t} - \mathsf{X} \mathsf{U}\_{k,p,t} \end{pmatrix} \times \left( \mathsf{P} P\_{a\_1, p,t} - \mathsf{T} \mathsf{N}\_{r,p} \right) \right] \\ -\left[ \begin{pmatrix} \sum\_{(i,m) \in \mathsf{Raut}\_{i,k,m}} \mathsf{X} P\_{i,k,m,p,t} \end{pmatrix} \times \left( \begin{matrix} \mathsf{P} P\_{a\_1, p,t} - \mathsf{T} \mathsf{N}\_{r,p} \end{matrix} \right) \right] \\ \left( \begin{matrix} \sum\_{(i,m) \in \mathsf{Raut}\_{i,k,m}} \mathsf{X} S\_{j,k,m,p} \end{matrix} \right) \times \left( \begin{matrix} \mathsf{P} P\_{a\_1, p,t} \end{matrix} \right) \\ \left( \begin{matrix} \sum\_{(i,m) \in \mathsf{Raut}\_{i,k,m}} \mathsf{T} P\_{i,k,m,p,t} \end{matrix} \right) \right] \end{aligned} \tag{13}$$

$$\forall k \in K, p \in P, t \in T, a\_1 = r = r \\ \text{return} \\ \text{l}, a\_2 = r \\ \text{e} \\ \text{f} \\ \text{r} \\ \text{r} \\ y, a\_3 = \text{depot} \\ \tag{13}$$

$$\text{CE}\_{i,p,t} = \text{XE}\_{i,p,t} \times \text{TN}\_{l,p} \quad \forall i \in I, p \in P, t \in T, r = \text{export} \tag{14}$$

$$\text{CLI}\_{h,p,t} = \text{XI}\_{h,p,t} \times \left( \text{PP}\_{a,p} + \text{TN}\_{r,p} \right) \quad \forall h \in \text{H}, p \in \text{P}, t \in T, a = r = import \tag{15}$$

$$\text{CPD}\_{i,p,t} = \sum\_{(l,m)\in\text{Raut}\mathfrak{c}\_{l,l,m}} \left( \text{XP}\_{i,l,m,p,t} \times \text{CT}\_{m,p} \times \text{Dist}\_{i,l,m} \right) \quad \forall i \in I, p \in \text{P}, t \in T \tag{16}$$

$$\text{CSD}\_{j,p,t} = \sum\_{(k,m)\in\text{Ront}\_{j,k,m}} \left( \text{XS}\_{j,k,m,p,t} \times \text{CT}\_{m,p} \times \text{Dist}\_{j,k,m} \right) \quad \forall j \in I, p \in P, t \in T \tag{17}$$

$$\text{CIO}\_{i,t} = \text{CKI} \times \text{IO}\_{i,t} \times \text{PO}\_{t} \quad \forall i \in I, t \in T, p = \text{oil} \tag{18}$$

$$\text{CIR}\_{i,p,t} = \text{CII} \times \text{IR}\_{i,p,t} \times \text{PP}\_{a,p,t} \quad \forall i \in I, p \in \mathcal{P}, t \in T, a = r efficiency \tag{19}$$

$$\text{CID}\_{j,p,t} = \text{CKI} \times \text{ID}\_{j,p,t} \times \text{PP}\_{a,p,t} \quad \forall j \in I, p \in P, t \in T, a = \text{depot} \tag{20}$$

$$\text{CIM}\_{k,p,t} = \text{CKI} \times \text{IM}\_{k,p,t} \times \text{PP}\_{a,p,t} \quad \forall k \in \mathcal{K}, p \in \mathcal{P}, t \in T, a = \text{retail} \tag{21}$$

$$\text{CVID}\_{k,p,t} = \text{XLI}\_{k,p,t} \times \text{TN}\_{\tau p} \quad \forall k \in \text{K}, p \in P, t \in T, r = \text{num} \, t \tag{22}$$

Equation (11) defines the refinery revenues obtained from the sales of refined products. Equations (12) and (13) determine the storage depot margins and the local market margins, respectively, being the difference between the revenues from selling refined products and the costs for buying these ones. Equations (14) and (15) specify, respectively, the costs of exporting and importing refined products. Equations (16) and (17) determine the costs for dispatching refined products through the primary and secondary distributions, in that order. Equation (18) defines the crude oil inventory costs at oil refineries, whereas Equations (19)–(21) define, respectively, the refined product inventory costs at oil

refineries, storage depots and local markets. Equation (22) represents the costs for not supplying the required product demand.

#### *4.2. Equations of the Recourse Problem* QU(·)

Equations (23)–(45) of the second-stage problem QU(·) depend on the worst-case realization of the stochastic parameters and are used to model the refinery operation, network flow allocation, inventory policy, and so on:

$$XOP\_{i,t} \le PC\_i \quad \forall i \in I, t \in T \tag{23}$$

$$XR\_{i,p,t} = XOP\_{i,t} \times YF\_{i,p} \quad \forall i \in I, p \in P, t \in T \tag{24}$$

$$IO\_{i,t} = ISO\_i + XCO\_{i,t} - XOP\_{i,t} \quad \forall i \in I, t \in T, t = t\_1 \tag{25}$$

$$IO\_{i,t} = IO\_{i,t-1} + XCO\_{i,t} - XOP\_{i,t} \quad \forall i \in I, t \in T, t > t\_1 \tag{26}$$

$$IR\_{i,p,t} = ISP\_{i,p} + XR\_{i,p,t} + XI\_{i,p,t} - \sum\_{(l,m)\in Rautc\_{i,l,m}} XP\_{i,l,m,p,t} - XE\_{i,p,t}$$

$$\forall i \in I, p \in P, t \in T, t = t\_1 \tag{27}$$

$$IR\_{i,p,t} = IR\_{i,p,t-1} + XR\_{i,p,t} + XI\_{i,p,t} - \sum\_{(l,m)\in Raut\varepsilon\_{i,l,m}} XP\_{i,l,m,p,t} - XE\_{i,p,t}$$

$$\forall i \in I, p \in P, t \in T, t > t\_1 \tag{28}$$

$$I\text{ID}\_{j,p,t} = \text{ISP}\_{j,p} + \text{XI}\_{j,p,t} + \sum\_{\substack{(i,m)\in\text{Raut}\ t:i,m}} \text{XP}\_{i,j,m,p,t} - \sum\_{\substack{(k,m)\in\text{Raut}\ t\_{j,k,m}}} \text{XS}\_{j,k,m,p,t} \tag{20}$$
 
$$\forall i\in I\text{ }n\in\text{P}\ t\in T\ t-t. \tag{21}$$

$$\forall j \in \mathcal{J}, p \in \mathcal{P}, t \in T, t = t\_1 \tag{29}$$

$$ID\_{j, p, t} = ID\_{j, p, t-1} + XI\_{j, p, t} + \sum\_{\substack{i \in \mathcal{I} \ i \in \mathcal{I} \ s}} XP\_{i, j, m, p, t} - \sum\_{\substack{i \in \mathcal{I} \ s \in T, t \in T}} XS\_{j, k, m, p, t}$$

$$(i, \mathfrak{m}) \in \mathbb{R} \\ \operatorname{aut}\_{i, j, \mathfrak{m}} \quad \quad \quad \quad (k, \mathfrak{m}) \in \operatorname{Raut}\_{j, k, \mathfrak{m}}$$

$$\forall j \in J, p \in P, t \in T, t > t\_1 \tag{30}$$

$$IM\_{k,p,l} = ISP\_{k,p} + \sum\_{\substack{(\downarrow m) \text{ or } Raut\_{\langle k,m \rangle}}} XP\_{\langle k,m,p,l \rangle} + \sum\_{\substack{(\downarrow m) \text{ or } Raut\_{\langle k,m \rangle}}} XS\_{\langle k,m,p,l \rangle} + XII\_{k,p,l} - XRP\_{k,p,l} + \Delta DP\_{k,p,l}$$

$$\forall k \in K, p \in P, t \in T, t = t\_1 \tag{31}$$

$$IM\_{k,p,l} = IM\_{k,p,l-1} + \sum\_{\substack{(i,m)\in Raut \ (k,m)}}XP\_{i,k,m,p,l} + \sum\_{\substack{(j,m)\in Raut \ (j,m)}}XS\_{j,k,m,p,l} + XI\_{k,p,l} - XRP\_{k,p,l} + \Delta DP\_{k,p,l}$$

$$\forall k \in \mathbb{K}, p \in P, t \in T, t > t\_1 \tag{32}$$

$$XRP\_{k,p,t} \le \left(\widehat{DP}\_{k,p,t} + \Lambda DP\_{k,p,t}\right) \quad \forall k \in K, p \in P, t \in T \tag{33}$$

$$\text{U}\,\,\text{SSO}\_{i}\times\text{SCO}\_{i}\le\text{IO}\_{i,t}\le\text{SCO}\_{i} \quad \forall i\in I, t\in T\tag{34}$$

$$\text{SSP}\_{i} \times \text{SC}\_{i,p} \le \text{IR}\_{i,p,t} \le \text{SC}\_{i,p} \quad \forall i \in I, p \in P, t \in T \tag{35}$$

$$\text{SSP}\_{j} \times \text{SC}\_{j,p} \le \text{ID}\_{j,p}, t \le \text{SC}\_{j,p} \quad \forall j \in \text{J}, p \in \text{P}, t \in T \tag{36}$$

$$\text{SSP}\_k \times \text{SC}\_{k,p} \le \text{IM}\_{k,p,t} \le \text{SC}\_{k,p} \quad \forall k \in K, p \in \text{P}, t \in T \tag{37}$$

$$\sum\_{l,l \in Rautc\_{i,l,m}} XP\_{i,l,m,p,t} \le ASPD\_{i,m,p} \quad \forall i \in I, m \in M, p \in P, t \in T \tag{38}$$

$$\sum\_{k \in \text{Raut}\_{j,k,m}} \text{XS}\_{j,k,m,p,t} \le ASSD\_{j,m,p} \quad \forall j \in \text{J, } m \in M, p \in P, t \in T \tag{39}$$

$$\sum\_{p \in P} XP\_{i,l,m,p,t} \le A \text{CPD}\_{i,l,m} \quad \forall (i,l,m) \in Rout\_{i,l,m}, t \in T \tag{40}$$

$$\sum\_{t:p\in\mathcal{P}}XS\_{j,k,m,p,t} \le ACSD\_{j,k,m} \quad \forall (j,k,m) \in Rout\_{j,k,m}, t \in T \tag{41}$$

$$\sum\_{t \in (j,m) \in Raut \cdot \varepsilon\_{i,j,m}} XP\_{i,j,m,p,t} + XI\_{j,p,t} \le SC\_{j,p} \times TCM\_{j,p} \quad \forall j \in J, p \in P, t \in T \tag{42}$$

$$\sum\_{(\cdot,m)\in\text{Raut}\_{\cdot,k,m}} \text{XP}\_{i,k,m,p,t} + \sum\_{(\cdot,m)\in\text{Raut}\_{\cdot,k,m}} \text{XS}\_{j,k,m,p,t} \le \text{SC}\_{k,p} \times \text{TCM}\_{k,p} \quad \forall k \in \text{K}, p \in \text{P}, t \in T \tag{43}$$

$$XP\_{i,l,m,p,t} = 0 \quad \forall i \in I, l \in L, p \in P, t \in T, m = rand, \mu\_{i,j,m} > MTD \tag{44}$$

$$\forall X S\_{j,k,m,p,t} = 0 \quad \forall j \in f, k \in K, p \in P, t \in T, m = r \\ \text{and} \\ \mu\_{j,k,m} > MTD \tag{45}$$

Equation (23) limits the oil processed at oil refineries. Equation (24) defines the fractions obtained through the crude oil distillation process. Equations (25) and (26) control the crude oil balance, as well as Equations (27) and (28) handle the refined product balance, at oil refineries. Likewise, Equations (29) and (30) control the refined product balance at storage depots, as well as Equations (31) and (32) do it at local markets. Equation (33) limits the product supply at local markets by the nominal product demand plus the associated deviation. Equations (34)–(37) determine the storage capacities for the facilities in the network. Equation (38) limits the capacity of supplying refined products in the primary distribution, as well as Equation (39) does so in the secondary distribution. Equation (40) defines the arc capacities for the primary distribution, whereas Equation (41) for the secondary distribution. Equations (42) and (43) determine the receiving capacities for the storage depots and local markets, respectively. Equation (44) and (45) suppress the primary and secondary distributions for long distances, respectively.

#### *4.3. Definition of Uncertainty Set* U

The uncertainty set is the key building block of the robust modeling framework [7]. It is assumed that the deviation of the product demand Δ*DPk*,*p*,*<sup>t</sup>* takes values within a convex and budget uncertainty set (Equations (46)–(48)) as follows—see [31]:

$$
\Delta DP\_{k,p,t} = \overline{DP}\_{k,p,t} - \overline{DP}\_{k,p,t} \quad \forall k \in K, p \in P, t \in T \tag{46}
$$

$$\left| \Delta DP\_{k,p,t} \right| \le \Delta DP\_{k,p,t}^{\max} \quad \forall k \in K, p \in P, t \in T \tag{47}$$

$$\sum\_{t} \frac{\left| \Delta DP\_{k,p,t} \right|}{\Delta DP\_{k,p,t}^{\max}} \le \Gamma\_{k,p} \quad \forall k \in \mathbb{K}, p \in P \tag{48}$$

Equation (46) defines the demand deviation Δ*DPk*,*p*,*<sup>t</sup>* as the difference between the uncertain value of the product demand *DPk*,*p*,*<sup>t</sup>* and the nominal value of the demand *DPk*,*p*,*t*. Equation (47) models the symmetrical and bounded range for the deviations of the product demand Δ*DPk*,*p*,*t*. Equation (48) ensures that the sum of the normalized absolute values of the demand deviations for each pair of location and product, across all time points *t*, must not exceed the user-defined budget of uncertainty Γ*k*,*p*. Notice that the normalized absolute value corresponds to the quotient between the absolute and maximum values of the demand deviation, while the uncertain value of the product demand refers to the unknown realization of the product demand. Appendix E includes an example showing how to define the vertices of a polyhedral budget uncertainty set.

#### *4.4. The Adaptive Robust Formulation*

In this section, the above robust formulation, that is, the set of Equations (8)–(48), is reformulated to be fully adaptive to the realizations of the uncertain deviations of the product demands Δ*DPk*,*p*,*t*. In order to transform the prior robust formulation into a single-level approach, the worst-case recourse profit QU(·) may be denoted by an auxiliary variable α, which must be lower or equal to the difference among the financial items—see Equation (49). Moreover, the set of recourse decisions of the second-stage problem QU(·) is reformulated as a function of the uncertain deviations of the product demands Δ*DPk*,*p*,*t*, which are defined by the uncertainty set U (Equations (46)–(48)). Because this uncertainty set U defines a polyhedron with infinite number of points, this reformulation would result in a recourse problem with an infinite number of variables and equations. For this reason, the uncertainty set U must be partitioned into *K* parts in order to make the resulting formulation tractable. According to Bertsimas et al. [25], only an extreme point (i.e., a vertex) of the defined polyhedron can be part of the optimal solution of the recourse problem QU(·). Hence, only the set of vertices *v* ∈ *V* of the polyhedron defined through the uncertainty set U is considered in this adjustable robust formulation. Therefore, the adjustable robust formulation is defined in the following set of Equations (49)–(86):

$$\max\_{\mathcal{R}, \mathcal{H}, \alpha} \mathsf{WCP} = -\sum\_{i \in I, t \in T} \mathsf{XCO}\_{i,t} \times \mathsf{PO}\_{t} + \alpha \tag{49}$$

$$\text{dist.}XCO\_{i,t} \le DO\_{i,t} \forall i \in I, t \in T \tag{50}$$

<sup>α</sup> <sup>≤</sup> min <sup>Δ</sup>*DPk*,*p*,*t*∈U max R∈F (H,Δ*DP*max *k*,*p*,*t* ) ⎛ ⎜⎜⎜⎜⎝ *i*∈*I*,*t*∈*T RRi*,*t*,*<sup>v</sup>* + *j*∈*J*,*p*∈*P*,*t*∈*T MDj*,*p*,*t*,*<sup>v</sup>* + *k*∈*K*,*p*∈*P*,*t*∈*T MCk*,*p*,*t*,*<sup>v</sup>* <sup>−</sup> *i*∈*I*,*p*∈*P*,*t*∈*T CEi*,*p*,*t*,*<sup>v</sup>* <sup>−</sup> *h*∈*H*,*p*∈*P*,*t*∈*T CIh*,*p*,*t*,*<sup>v</sup>* <sup>−</sup> *i*∈*I*,*p*∈*P*,*t*∈*T CPDi*,*p*,*t*,*<sup>v</sup>* <sup>−</sup> *j*∈*J*,*p*∈*P*,*t*∈*T CSDj*,*p*,*t*,*<sup>v</sup>* <sup>−</sup> *i*∈*I*,*t*∈*T CIOi*,*t*,*<sup>v</sup>* <sup>−</sup> *i*∈*I*,*p*∈*P*,*t*∈*T CIRi*,*p*,*t*,*<sup>v</sup>* <sup>−</sup> *j*∈*J*,*p*∈*P*,*t*∈*T CIDj*,*p*,*t*,*<sup>v</sup>* <sup>−</sup> *k*∈*K*,*p*∈*P*,*t*∈*T CIMk*,*p*,*t*,*<sup>v</sup>* <sup>−</sup> *k*∈*K*,*p*∈*P*,*t*∈*T CUDk*,*p*,*t*,*<sup>v</sup>* ⎞ ⎟⎟⎟⎟⎠ ∀*v* ∈ *V* (51)

$$RR\_{i,t,\upsilon} = \sum\_{p \in P} \left[ XR\_{i,p,t,\upsilon} \times \left( PP\_{a,p,t} - T\aleph\_{r,p} \right) \right] \quad \forall i \in I, t \in T, \upsilon \in V, a = r = r\epsilon\\primary \tag{52}$$

$$MD\_{j,p,t,\nu} = \sum\_{(k,m)\in Raut\_{i,j,m}} \text{XS}\_{j,k,m,p,t,\nu} \times \left( PP\_{a,p,t} - TN\_{r,p} \right) - \sum\_{(i,m)\in Raut\_{i,j,m}} \text{XP}\_{i,j,m,p,t,\nu} \times \text{PP}\_{a2,p,t} \tag{53}$$

$$\forall j \in \mathcal{J}, p \in \mathcal{P}, t \in T, \nu \in V, a = r = \text{depot}, a\_2 = r \text{f}\\ \text{nary} \tag{53}$$

$$\begin{aligned} ML\_{k,p,t,\nu} &= \left[ \left( XRP\_{k,p,t,\nu} - XIL\_{k,p,t,\nu} \right) \times \left( PP\_{a\_1,p,t} - TN\_{r,p} \right) \right] \\ -\left[ \left( \sum\_{(j,m)\in Raut\varepsilon\_{j,k,m}} XP\_{i,k,m,p,t,\nu} \right) \times PP\_{a\_2,p,t} + \left( \sum\_{(j,m)\in Raut\varepsilon\_{j,k,m}} XS\_{j,k,m,p,t,\nu} \right) \times PP\_{a\_3,p,t} \right] \end{aligned}$$

$$\forall k \in \mathbb{K}, p \in P, t \in T, v \in V, a\_1 = r = \text{retail}, a\_2 = r \text{e} \\ \text{inner}, a\_3 = \text{depot} \tag{54}$$

$$CE\_{i,p,t,\nu} = XE\_{i,p,t,\nu} \times TN\_{r,p} \quad \forall i \in I, p \in P, t \in T, \nu \in V, r = \text{export} \tag{55}$$

$$\text{C}\newline \text{CI}\_{h,p,t,\upsilon} = \text{XI}\_{h,p,t,\upsilon} \times \left( \text{PP}\_{a,p} + \text{TN}\_{r,p} \right) \quad \forall h \in H, p \in P, t \in T, \upsilon \in V, a = r = \text{import} \tag{56}$$

$$\text{CPD}\_{i,p,t,\upsilon} = \sum\_{(l,m)\in\text{Raut}\_{i,l,m}} \left( \text{XP}\_{i,l,m,p,t,\upsilon} \times \text{CT}\_{m,p} \times \text{Dist}\_{i,l,m} \right) \quad \forall i \in I, p \in P, t \in T, \upsilon \in V \tag{57}$$

$$\text{CSD}\_{j,p,t,\nu} = \sum\_{(k,m)\in \text{Ront}\_{j,k,m}} \left( \text{XS}\_{j,k,m,p,t,\nu} \times \text{CT}\_{m,p} \times \text{Dist}\_{j,k,m} \right) \quad \forall j \in \text{J}, p \in P, t \in T, \upsilon \in V \tag{58}$$

$$\text{CIO}\_{i,t,\upsilon} = \text{CKI} \times \text{IO}\_{i,t,\upsilon} \times \text{PO}\_{t} \quad \forall i \in I, t \in T, \upsilon \in V, p = oil \tag{59}$$

$$\text{CIR}\_{i,p,t,\upsilon} = \text{CKI} \times \text{IR}\_{i,p,t,\upsilon} \times \text{PP}\_{a,p,t} \quad \forall i \in I, p \in P, t \in T, \upsilon \in V, a = r \\ \text{primary} \tag{60}$$

$$\text{CID}\_{j,p,t,\upsilon} = \text{CKI} \times \text{ID}\_{j,p,t,\upsilon} \times \text{PP}\_{a,p,t} \quad \forall j \in \mathcal{J}, p \in \mathcal{P}, t \in T, \upsilon \in V, a = \text{depot} \tag{61}$$

$$\text{CIM}\_{k,p,t,\upsilon} = \text{CKI} \times \text{IM}\_{k,p,t,\upsilon} \times \text{PP}\_{a,p,t} \quad \forall k \in K, p \in P, t \in T, \upsilon \in V, a = \text{retail} \tag{62}$$

$$\text{CII}D\_{k,p,t,\upsilon} = \text{XLI}\_{k,p,t,\upsilon} \times \text{TN}\_{r,p} \quad \forall k \in \mathcal{K}, p \in \mathcal{P}, t \in T, \upsilon \in V, r = \text{num} \,\text{et} \tag{63}$$

$$XOP\_{i,t,\upsilon} \le PC\_i \quad \forall i \in I, t \in T, \upsilon \in V \tag{64}$$

$$XR\_{i,p,t,\nu} = XOP\_{i,t,\nu} \times YF\_{i,p} \quad \forall i \in I, p \in P, t \in T, \upsilon \in V \tag{65}$$

$$IO\_{i,t,\upsilon} = ISO\_i + XCO\_{i,t} - XOP\_{i,t,\upsilon} \quad \forall i \in I, \upsilon \in V, t \in T, t = t\_1 \tag{66}$$

$$IO\_{i,t,v} = IO\_{i,t-1,v} + XCO\_{i,t} - XOP\_{i,t,v} \quad \forall i \in I, v \in V, t \in T, t > t\_1 \tag{67}$$

$$IR\_{i,p,t,\nu} = ISP\_{i,p} + XR\_{i,p,t,\nu} + XI\_{i,p,t,\nu} - \sum\_{(l,m)\in Rautc\_{i,l,m}} XP\_{i,l,m,p,t,\nu} - XE\_{i,p,t,\nu}$$

$$\forall i \in I, p \in P, v \in V, t \in T, t = t\_1 \tag{68}$$

$$IR\_{i,p,t,\nu} = IR\_{i,p,t-1,\nu} + XR\_{i,p,t,\nu} + XI\_{i,p,t,\nu} - \sum\_{(l,m)\in Raut\varepsilon\_{i,l,m}}XP\_{i,l,m,p,t,\nu} - XE\_{i,p,t,\nu}$$

$$\forall i \in I, p \in P, v \in V, t \in T, t > t\_1 \tag{69}$$

$$\text{I}\,\text{ID}\_{j,p,t,\text{pr}} = \text{I}\,\text{SP}\_{j,p} + \text{XI}\_{j,p,t,\text{pr}} + \sum\_{(i,m)\in\text{Rautte}\_{i,j,m}} \text{XP}\_{i,j,m,p,t,\text{pr}} - \sum\_{(k,m)\in\text{Rautte}\_{j,k,m}} \text{XS}\_{j,k,m,p,t,\text{pr}}\tag{70.1}$$
 
$$\text{V}\,\text{V}=\text{I}\,\text{V}\,\text{V}=\text{I}\,\text{V}\,\text{V}=\text{I}\,\text{V}\,\text{V}=\text{I}$$

$$\forall j \in I, p \in P, v \in V, t \in T, t = t\_1 \tag{70}$$

$$ID\_{j,p,t,\nu} = ID\_{j,p,t-1,\nu} + XI\_{j,p,t,\nu} + \sum\_{(i,m)\in Raut\varepsilon\_{i,j,m}} XP\_{i,j,m,p,t,\nu} - \sum\_{(k,m)\in Raut\varepsilon\_{j,k,m}} XS\_{j,k,m,p,t,\nu}$$

$$V\_{\varepsilon} = I\_{\varepsilon,m} - D\_{\varepsilon} = V\_{\varepsilon}t - T\_{\varepsilon}t - \varepsilon\_{\varepsilon}$$

$$\forall j \in \{ , p \in P, \upsilon \in V, t \in T, t > t\_1 \tag{71}$$

$$\begin{split} IM\_{k,p,t,\nu} &= IM\_{k,p,t-1,\nu} + \sum\_{(i,m)\in Raut\varepsilon\_{i,k,m}} XP\_{i,k,m,p,t,\nu} + \sum\_{(j,m)\in Raut\varepsilon\_{j,k,m}} XS\_{j,k,m,p,t,\nu} \\ &+ XL\_{k,p,t,\nu} - XRP\_{k,p,t,\nu} + \Delta DP\_{k,p,t,\nu} \end{split} \tag{72}$$

$$\forall k \in \mathcal{K}, p \in P, v \in V, t \in T, t > t\_1 \tag{73}$$

$$XRP\_{k,p,t,\upsilon} \le \left(\widehat{DP}\_{k,p,t,\upsilon} + \Delta DP\_{k,p,t,\upsilon}\right) \quad \forall k \in K, p \in P, t \in T, \upsilon \in V \tag{74}$$

$$\text{SSC}\_{i} \times \text{SCO}\_{i} \le IO\_{i,t,\upsilon} \le \text{SCO}\_{i} \quad \forall i \in I, t \in T, \upsilon \in V \tag{75}$$

$$\text{SSP}\_i \times \text{SC}\_{i,p} \le IR\_{i,p,t,\upsilon} \le \text{SC}\_{i,p} \quad \forall i \in I, p \in P, t \in T, \upsilon \in V \tag{76}$$

$$\text{LSP}\_{j} \times \text{SC}\_{j,p} \le \text{ID}\_{j,p,t,v} \le \text{SC}\_{j,p} \quad \forall j \in I, p \in P, t \in T, v \in V \tag{77}$$

$$\text{SSP}\_k \times \text{SC}\_{k,p} \le \text{IM}\_{k,p,t,\upsilon} \le \text{SC}\_{k,p} \quad \forall k \in K, p \in P, t \in T, \upsilon \in V \tag{78}$$

$$\sum\_{l \in \text{Ront}\_{i,l,m}} \text{XP}\_{i,l,m,p,t,\upsilon} \le \text{ASPD}\_{i,m,p} \quad \forall i \in I, m \in \text{M}, p \in P, t \in T, \upsilon \in V \tag{79}$$

$$\sum\_{\substack{j,k \in Raut\_{j,k,m}}} X S\_{j,k,m,p,t,\nu} \le ASSD\_{j,m,p} \quad \forall j \in I, m \in M, p \in P, t \in T, \nu \in V \tag{80}$$

$$\sum\_{p \in P} XP\_{i,l,m,p,t,\upsilon} \le ACDD\_{i,l,m} \quad \forall (i,l,m) \in Raute\_{i,l,m}, t \in T, \upsilon \in V \tag{81}$$

$$\sum\_{p \in P} X S\_{j,k,m,p,t,\nu} \le A C S D\_{j,k,m} \quad \forall (j,k,m) \in \text{Raut}\_{j,k,m}, t \in T, \upsilon \in V \tag{82}$$

$$\sum\_{(i,m)\in Raut\_{i,j,m}} X P\_{i,j,m,p,t,\nu} + X I\_{j,p,t,\nu} \le S C\_{j,p} \times T C M\_{j,p} \quad \forall j \in I, p \in P, t \in T, \upsilon \in V \tag{83}$$

$$\sum\_{(j,m)\in Rmt\varepsilon\_{j,k,m}} X P\_{i,k,m,p,t,\nu} + \sum\_{(j,m)\in Rmt\varepsilon\_{j,k,m}} X S\_{j,k,m,p,t,\nu} \le S C\_{k,p} \times T\mathcal{C}M\_{k,p} \quad \forall k \in \mathcal{K}, p \in P, t \in T, \nu \in V \tag{84}$$

$$XP\_{i,l,m,p,t,\upsilon} = 0 \quad \forall i \in I, l \in L, p \in P, t \in T, \upsilon \in V, m = r \text{ord}, \mu\_{i,j,m} > MTD \tag{85}$$

$$XS\_{j,k,m,p,t,\upsilon} = 0 \quad \forall j \in \mathcal{J}, k \in \mathcal{K}, p \in P, t \in T, \upsilon \in V, m = \text{round}\_r \mu\_{j,k,m} > MTD \tag{86}$$

#### **5. Case Study**

The ARO model is applied to a real case study on a refined products network distribution in the Portuguese oil industry, which was originally characterized by Lima et al. [1]. Figure 1 represents the network under study.

**Figure 1.** Portuguese downstream oil supply chain (from Lima et al., 2018).

The network comprises two oil refineries (Sines and Matosinhos), three distribution centers (Boa nova, CLC and Mitrena), 278 local markets (Portuguese cities) and four transportation modes (pipeline transport, railway, roadway and waterway). In this network, the oil refineries import the same type of crude oil and produce eight refined products, namely, jet fuel, diesel, propane, butane, fuel oil, gas oil, gasoline 95 and gasoline 98. These refined products are then transported using pipelines, tanker ships, tank wagons and tank trucks from the oil refineries, and using only tank trucks from the distribution centers. Roadway is only considered for distances less than or equal to 250 km. The planning horizon considered is 6 months, which is discretized into monthly cycles. For the sake of comparison, the same database employed by Lima et al. [1] is used in the case study.

#### **6. Results and Discussion**

The proposed two-stage ARO model is implemented in GAMS 24.5.6 and solved using CPLEX 12.6 on an Intel(R) Xeon(R) processor CPU E5-2660 v3 @ with 2.60 GHz (two processors) and 64 GB RAM memory.

#### *6.1. Setup of the ARO Model*

In the practical application, the ARO model was run under a specific configuration of the budget uncertainty set (Equations (46)–(48)). In other words, the ARO model was solved after specifying the parameters of the polyhedral uncertainty set, namely, the nominal value *DPk*,*p*,*<sup>t</sup>* and the maximum deviation in absolute value Δ*DPmax <sup>k</sup>*,*p*,*<sup>t</sup>* of the uncertain product demand *DPk*,*p*,*t*, besides the budget of uncertainty Γ*k*,*p*. The values of these parameters were estimated from the dataset described by Lima et al. [1]. Therefore, the maximum deviation Δ*DPmax <sup>k</sup>*,*p*,*<sup>t</sup>* was estimated as 10% of the nominal value *DPk*,*p*,*t*, while the budget of uncertainty Γ*k*,*<sup>p</sup>* was set to the average of the total relative deviations over the time horizon for the pairs of local market and product. To be more precise, for each pair of location and product (*k*, *p*), the absolute value of the product demand deviation Δ*DPk*,*p*,*<sup>t</sup>* is normalized to the associated maximum deviation Δ*DPmax <sup>k</sup>*,*p*,*<sup>t</sup>* at each time point. We then sum of all the so normalized deviations over the time horizon. Finally, we average these normalized values over all pairs of products and locations, and the result (let's denote it by Γ*k*,*p*) is considered as the value of the budget of uncertainty, that is,Γ*k*,*<sup>p</sup>* for all (*k*, *p*). This results inΓ*k*,*<sup>p</sup>* = 0.895, for all (*k*, *p*), which leads to uncertainty sets (i.e., Equations (46)–(48)) consisting in polyhedrons of 12 vertices.

#### *6.2. Comparison with the Developed Approaches*

In this section, we compare the results given by the proposed adjustable robust optimization (ARO) model with those obtained using the equivalent non-adjustable robust optimization (NARO), stochastic programming (SP) and the deterministic (DM) models. The mathematical formulations, the model notations and the specific settings for the DM, NARO and SP models can be found in Appendix B, Appendix C, Appendix D, respectively.

#### 6.2.1. Computational Performance

Table 1 presents the problem size, economic performance and solving time for all the optimization models.


**Table 1.** Statistics.

From Table 1, we can conclude that the economic performance (i.e., the profit) declines when uncertainty in demand for refined products is taken into consideration by the optimization model. On the other hand, a model in which all the parameters are considered as deterministic (DM) is not realistic, and it should consider random parameters to provide robustness and consistency to the decision maker.

In this way, the SP model aims at maximizing the expected profit for a set of scenarios, while the robust approaches maximize the total profit of the worst-case scenario. Compared to the ARO approach, the profit is higher about 0.16% for the SP approach. When comparing both robust approaches, the ARO model shows a better economic performance, providing a profit that is almost 0.80% higher than the one returned by the NARO model, which is the most conservative approach.

The problem size is also different for all the considered formulations. That is, the size of the deterministic model is smaller than the other methods, and thus it requires less computational time (1.484 CPU seconds). The NARO model does not increase too much the problem size and the computational time, because its size raises as more dual variables and robust constraints are incorporated into the problem, and such an increase is at the same scale as the number of the added random parameters, which in this case is only the demand for refined products. Nevertheless, much more computational burden is demanded to solve the much bigger ARO and SP models (95.3 and 236.8 CPU seconds, respectively). The size of the ARO model relies on the size of the polyhedral uncertainty set (Equations (46)–(48)), which in turn is determined by the budget of uncertainty Γ*k*,*<sup>p</sup>* in Equation (48). The size of the scenario-based SP model depends on the number of the considered scenarios, which in this case amounts to 243, in order to model the uncertainty in the stochastic parameter. All in all, the ARO approach does not increase the problem size as much as the SP approach when compared with the nominal problem (i.e., DM), having better solution performance, while it approximates better the economic result of the SP when comparing to the NARO approach.

#### 6.2.2. Insights about the Network Profitability

The contribution of the different activities to the network profitability is illustrated in Figure 2 for all the optimization models. The oil procurement costs and the refinery revenues are all included into the margins of the refinery sector, so that the comparisons can be possible with the margins of the distribution and market sectors. In general, the profit breakdown is not too affected when the case study is solved by the different models. As the only difference in these approaches is how uncertainty is modeled, the reason comes from that. Firstly, the same uncertainty set U is considered by both robust optimization approaches. Secondly, the stochastic approach provides an average result of the considered scenarios, which are properly defined in accordance with the level of uncertainty considered by the robust approaches. Lastly, the product demand uncertainty does not actually interfere too much in the profit breakdown.

**Figure 2.** Contribution margins of the network echelons to the profitability.

As it can be seen in Figure 2, the market sector delivers the largest contribution to the network profitability, followed by the refining and distribution sectors, respectively. When comparing both robust approaches, the profit breakdown differs. The NARO model increases the shares of the refining and market sectors and decreases the contribution of the distribution sector when compared with the ARO model. The SP and DM models result in the same profit breakdown, which is slightly different from the one determined by the ARO model, where the market contribution rises 1% and the distribution sector diminishes 1%, whereas the share of the refining sector remains unaltered.

Figure 3 compares the effective profitability contribution of the network sectors among the different optimization models. All the optimization approaches provide quite similar results for the refining sector. The same fact is not observed in the distribution and market sectors, where the NARO model returns the lowest economic performance. These facts corroborate that the NARO model is the most conservative approach to handle uncertainty in this case study.

**Figure 3.** Comparison of profitability contributions among the optimization models.

In order to get more insights about the network profitability, the cost breakdown information is also investigated. As abovementioned, refinery margins comprise the revenues for selling oil products and the cost for acquiring crude oil. The oil procurement cost is the major cost item, and accounts for at least 78.53% of the overall cost in all the developed approaches, while the other costs together account for at most 21.47% of the same total. The oil procurement cost amounts to 78.53%, 80.92%, 78.94% and 78.54% of the total network cost in the ARO, NARO, SP and DM approaches, respectively. The cost breakdown information of the studied models is illustrated in Figure 4. Even though the exportation and lost demand costs are included in the estimation of the cost breakdown, their slices of the overall cost are quite small, i.e., almost 0%, and hence they have been omitted to enhance the readability of the charts without loss of generality.

As shown in Figure 4, the total network cost is divided into the same proportions for both ARO and DM approaches, while the proportions are much different for the other approaches. For instance, the SP and NARO approaches have increased the oil procurement cost by 1% and 3%, respectively, and have decreased the total of the other costs by 1% and 3%, respectively, when compared with the previous approaches. In most cases, the importation cost is the second most relevant cost item, followed by the secondary and primary transportation costs, respectively. The inventory cost aggregates the costs for keeping inventories of oil and refined products at the network facilities and amounts to 1% of the overall cost.

Figure 5 compares the absolute values of the cost items among the optimization methods. Notice that the oil procurement cost is not displayed, because it is much larger than the other cost items, and assumes quite similar values in all the methods, i.e., €2,408,811,163 on average. As it can be seen in Figure 5, the ARO and DM methods have much similar network costs, which are usually higher than the cost items for the NARO and SP models. Notice that the export and unmet demand costs are much smaller than the other cost items, and hence they are displayed in a different scale. Even though the NARO presents the best performance with regards to the network costs, its poorer performance in the network margins makes it as the most conservative approach. In turn, the SP model is the less conservative model, but it is not so efficient as the others in order to fulfill the required demand, and thus the lost demand is larger.

**Figure 4.** Cost breakdown for the optimization models.

**Figure 5.** Cost items for the optimization models.

#### 6.2.3. Network Planning for the Refined Products Distribution

Figure 6 shows how the total volume, which is transported over the network, is divided into different proportions by using the optimization models. As it can be observed, the ARO and DM models return the same volume portions, what is fully consistent with the cost breakdown information in Figure 4. However, the SP and NARO models provide different results when compared with the previous models, but some similarities are also observed. For instance, the amount of oil processed also corresponds to the largest slice of the total network volume, and the associated percentages are 26% and 28% in the SP and DM models, respectively. The oil delivery refers to the second biggest piece of the total volume, and accounts for 23% in the ARO, SP and DM models and 24% in the NARO model. Notice that, at all the volume charts, the oil processed is higher than the oil delivery due to the initial oil inventory at oil refineries in the first-time point, which is consumed throughout the planning horizon.

**Figure 6.** Proportions of the total network volumes among the optimization models.

As it can be seen in Figure 6, the portion of volumes conveyed through primary distribution is considerably bigger than the piece of volumes distributed via secondary distribution. The secondary distribution costs are much bigger than the primary distributions cost as displayed in Figures 4 and 5, because the former is only performed via roadway, i.e., the most expensive way to convey refined products in the network, while the latter can be undertaken by any transportation mode. Similarly, the importation volumes are lower than the exportation volumes, but they generate much bigger costs due to the importation tariffs that are paid to bring the refined product into the Portuguese network. The inventory volumes also present a certain relevance in the network and account for 3% or 4% of the total network volume.

Figure 7 displays the actual volumes determined by solving the case study using the optimization approaches. The NARO model presents the most different network flows among the optimization models, as well as the worst overall performance. Although the NARO model defines to purchase and process the same volumes of crude oil seen in the other models, the network flows through the primary and secondary distributions are significantly shorter, besides the inventories are sensibly higher. Hence, the NARO model does not provide a network that is so profitable as the other models do — see Table 1. Notice that these other models present much similar performances. However, the SP model is the only one that has a so evident lost demand.

**Figure 7.** Network flows by solving the case study with different optimization models.

#### 6.2.4. General Aspects about the Developed Modeling Frameworks

All the proposed approaches are useful to handle uncertainty in DOSC problems, and some conclusions can be withdrawn from the previous analyses. Under a second-stage stochastic programming approach, the decision variables are separated into two different groups, i.e., the firstand second-stage variables. This is a troublesome task, which depends on the decision maker's knowledge about the problem under study [8]. Nonetheless, once this separation is successfully performed, it may improve the model robustness against infeasibilities caused by the realization of the random parameters, because the second-stage variables might be properly adjusted to any particular realization of uncertainty [6]. In general, stochastic programming with recourse might be a good option when the probability distribution of the random parameters can be obtained from the historical data, such that a set of scenarios, i.e., a scenario tree, can be generated to represent the underlying uncertainty [32]. The decision maker can precisely model uncertainty by eliminating the undesirable scenarios and specifying the most adequate scenarios. However, the assignment of probabilities to scenarios, as well as the definition of scenario tree frameworks, could not be easy. Additionally, a wide range of scenarios should be considered to model uncertainty adequately, which could result in large-scale or even intractable mathematical programs [31]. Hence, the use of decomposition methods and approximation schemes for their solution are usually employed to solve this class of optimization problems [33].

In contrast to stochastic programming, robust optimization does not assume that uncertainty has a probability distribution [27], but alternatively it assumes that uncertainty is represented through uncertainty sets [29]. In this way, the decision maker can represent uncertainty in the random parameters by defining their nominal values and variation amplitudes from the historical data. Another advantage of employing robust optimization is the computational tractability for solving numerous classes of uncertainty sets and problem types [25]. In this methodology, the decision maker aims at establishing a feasible solution for any realization of the random parameters in a given uncertainty set [34], while the decision maker can control the trade-off between robustness and performance by using a budget of uncertainty that is introduced in the prescribed set [35]. However, a single-stage (non-adjustable) robust optimization approach tends to be very conservative because all decisions are made before uncertainty is revealed [27] while a two-stage adjustable robust optimization has a higher modeling capability, in which the second-stage problem models decision making after the first-stage decisions are made and uncertainty is disclosed [28].

Generally, the proposed modeling frameworks showed to be efficient and effective to cope with DOSC problems under uncertainty. These methodologies present different goals and consider specific assumptions and simplifications to represent uncertainty. Therefore, the most adequate methodology depends on the considered problem, the dataset and the decision maker preferences.

#### **7. Conclusions**

In this paper, a two-stage adjustable robust optimization model is developed to deal with demand uncertainty in the tactical planning of refined products distribution in a downstream oil supply chain. The adjustable robust model is then compared with its non-adjustable, stochastic and deterministic counterparts, whose objectives are different, that is, the robust approaches concern to maximize the profit at the worst-case scenario, the stochastic approach aims at maximizing the expected profit for a set of scenarios, and the deterministic approach at maximizing the total profit of the nominal problem. However, all the optimization approaches provide comparable results in terms of economic performance and material flows, and the major discrepancies occur with regards to problem sizes and computational properties.

Specifically, the obtained results show that the non-adjustable model is the most conservative, while the stochastic model is in turn the least conservative. However, the main drawback of the stochastic approach is the limitation of problem size due to the computational burden, while the major advantage of the non-adjustable robust approach is that the problem size is not overly enlarged with regards to the nominal deterministic approach. In comparison, the adjustable robust model establishes a problem that is not as short as its non-adjustable counterpart, as well as is not so big as its stochastic counterpart, so that the model tractability issues are reasonable. All the optimization approaches provide different network flows, but too comparable. The adjustable approach presents the best performance in this respect among the developed approaches to cope with uncertainty, because it provides the highest service level in order to fulfill the required demand for refined products. In contrast, the non-adjustable approach has the most inferior performance over the supply chain, whose network flows are majorly lower in the primary and secondary distributions, for example.

In summary, all the developed optimization approaches are valuable to deal adequately with refined products distribution under uncertainty. Each approach has a specific objective and assumes distinct assumptions and simplifications in order to model and represent adequately uncertainty. Hence, the most appropriate method depends on the problem under study, as well as on the available dataset and on the decision maker preferences.

As future work, and as a direct extension of the present work, the two-stage ARO model could be further explored to include other uncertainty sets, as well as modeling more than one uncertainty type simultaneously, ex. crude oil and product prices. Also, other approaches to deal with uncertainty could be explored such as fuzzy programming and chance-constrained programming. In addition, Markov chain and game theory might be investigated and employed to cope with the stochastic parameters. Finally, the studied approaches to deal with uncertainty could be also applied to the strategic and tactical problem of the downstream oil network allowing for the design and planning of such system.

**Author Contributions:** Writing—original draft, C.L.; Major revision—review and edition, S.R. and A.B.-P.; Mathematical model analysis—significant comments and suggestions, J.M.M.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **Appendix A. Nomenclature**


**Table A1.** Model notation.

**Table A1.** *Cont.*


#### **Appendix B. Deterministic Mathematical Formulation**

The deterministic formulation is briefly reported below through its objective function (Equations (A1)–(A13)), network equations (Equations (A14)–(A37)) and model notation (Table A2). For the sake of brevity, the descriptions of the equations are omitted, but the adjustable robust optimization (ARO) counterpart can be consulted within the full paper for further details.

Objective function

$$\begin{array}{c} \text{maximize} \pi = \sum\_{i \in I, t \in T} MR\_{i, t} + \sum\_{j \in J, p \in P, t \in T} MD\_{j, p, t} + \sum\_{k \in K, p \in P, t \in T} MC\_{k, p, t} - \sum\_{i \in I, p \in P, t \in T} CE\_{i, p, t} \\ -\sum\_{\substack{h \in H, p \in P, t \in T}} CI\_{h, p, t} - \sum\_{\substack{i \in I, p \in P, t \in T}} PD\_{i, p, t} - \sum\_{j \in J, p \in P, t \in T} SD\_{j, p, t} \\ -\sum\_{i \in I, t \in T} CI\_{i, t} - \sum\_{\substack{i \in I, p \in P, t \in T}} CR\_{i, p, t} - \sum\_{j \in J, p \in P, t \in T} CD\_{j, p, t} \\ -\sum\_{k \in K, p \in P, t \in T} CM\_{k, p, t} - \sum\_{k \in K, p \in P, t \in T} CID\_{k, p, t} \\ \end{array} \tag{A1}$$

$$MR\_{i,t} = \sum\_{p \in P} \left[ XR\_{i,p,t} \times \left( PP\_{a,p,t} - TN\_{r,p} \right) \right] - \left[ XCO\_{i,t} \times PO\_t \right] \quad \forall i \in I, t \in T, a = r = r\_{\text{refner}} \tag{A2}$$

$$\begin{aligned} \textit{MD}\_{j,p,t} &= \sum\_{(k,m)\in\textit{Raut}\_{i,j,m}} \textit{XS}\_{j,k,m,p,t} \times \left( \textit{PP}\_{a,p,t} - \textit{TN}\_{r,p} \right) - \sum\_{(i,m)\in\textit{Raut}\_{i,j,m}} \textit{XP}\_{i,j,m,p,t} \times \textit{PP}\_{a2,p,t} \\ &\forall j \in \textit{J}, p \in \textit{P}, t \in \textit{T}, a = r = \textit{depot}, a\_2 = r \textit{refner} \, y \end{aligned} \tag{A3}$$

$$\begin{array}{c} \text{MC}\_{k,p,t} = \left[ \left( \text{XRP}\_{k,p,t} - \text{XLI}\_{k,p,t} \right) \times \left( \text{PP}\_{a\_1,p,t} - \text{TN}\_{r,p} \right) \right] - \left[ \left( \sum\_{(i,m) \in \text{Rautte}\_{i,k,m}} \text{XP}\_{i,k,m,p,t} \right) \times \text{PP}\_{i,k,m,p,t} \right] \\ \quad \times \text{PP}\_{a\_2,p,t} + \left( \sum\_{(j,m) \in \text{Rautte}\_{j,k,m}} \text{XS}\_{j,k,m,p,t} \right) \times \text{PP}\_{a3,p,t} \end{array} \tag{A4}$$

$$\forall k \in \mathbb{K}, p \in \mathcal{P}, t \in T, a\_1 = r = \text{retail}, a\_2 = \text{refn} \\ \text{rary}, a\_3 = \text{depot}$$

$$CE\_{i,p,t} = XE\_{i,p,t} \times TN\_{r,p} \quad \forall i \in I, p \in P, t \in T, r = \text{export} \tag{A5}$$

$$\text{CI}\_{\text{h},p,t} = \text{XI}\_{\text{h},p,t} \times \left( \text{PP}\_{a,p} + \text{TN}\_{r,p} \right) \quad \forall h \in H, p \in P, t \in T, a = r = \text{import} \tag{A6}$$

$$\text{CPD}\_{i,p,t} = \sum\_{(l,m)\in \text{Raut}\_{i,l,m}} \left( \text{XP}\_{i,l,m,p,t} \times \text{CT}\_{m,p} \times \text{Dist}\_{i,l,m} \right) \quad \forall i \in I, p \in P, t \in T \tag{A7}$$

$$\text{CSD}\_{j,p,t} = \sum\_{(k,m)\in \text{Ront}\_{j,k,m}} \left( \text{XS}\_{j,k,m,p,t} \times \text{CT}\_{m,p} \times \text{Dist}\_{j,k,m} \right) \quad \forall j \in I, p \in P, t \in T \tag{A8}$$

*CIOi*,*<sup>t</sup>* = *CKI* × *IOi*,*<sup>t</sup>* × *POt* ∀*i* ∈ *I*, *t* ∈ *T*, *p* = *oil* (A9)

$$\text{CIR}\_{i,p,t} = \text{CII} \times \text{IR}\_{i,p,t} \times \text{PP}\_{a,p,t} \quad \forall i \in I, p \in \mathcal{P}, t \in T, a = r efficiency \tag{A10}$$

$$\text{CID}\_{j,p,t} = \text{CKI} \times \text{ID}\_{j,p,t} \times \text{PP}\_{a,p,t} \quad \forall j \in \mathcal{J}, p \in \mathcal{P}, t \in T, a = \text{depot} \tag{A11}$$

$$\text{CIM}\_{k,p,t} = \text{CKI} \times \text{IM}\_{k,p,t} \times \text{PP}\_{a,p,t} \quad \forall k \in \mathcal{K}, p \in \mathcal{P}, t \in T, a = \text{retail} \tag{A12}$$

$$\text{Cull}\_{k,p,t} = \text{X} \text{LI}\_{k,p,t} \times \text{TN}\_{r,p} \quad \forall k \in \text{K}, p \in \text{P}, t \in T, r = \text{num} \text{et} \tag{A13}$$

Network equations

$$XOP\_{i,t} \le PC\_i \quad \forall i \in I\_\prime\\t \in T \tag{A14}$$

$$X R\_{i, p, t} = X O P\_{i, t} \times Y F\_{i, p} \quad \forall i \in I, p \in P, t \in T \tag{A15}$$

$$XCO\_{i,t} \le DO\_{i,t} \quad \forall i \in I, t \in T \tag{A16}$$

$$IO\_{i,t} = ISO\_i + XCO\_{i,t} - XOP\_{i,t} \quad \forall i \in I, t \in T, t = t\_1 \tag{A17}$$

$$IO\_{i,t} = IO\_{i,t-1} + XCO\_{i,t} - XOP\_{i,t} \quad \forall i \in I, t \in T, t > t\_1 \tag{A18}$$

$$IR\_{i,p,t} = ISP\_{i,p} + XR\_{i,p,t} + XI\_{i,p,t} - \sum\_{(l,m)\in Raut\_{i,l,m}} XP\_{i,l,m,p,t} - XE\_{i,p,t} \quad \forall i \in I, p \in P, t \in T, t = t\_1 \tag{A19}$$

$$IR\_{i,p,t} = IR\_{i,p,t-1} + XR\_{i,p,t} + XI\_{i,p,t} - \sum\_{(l,m)\in Rautc\_{i,lm}} XP\_{i,l,m,p,t} - XE\_{i,p,t} \quad \forall i \in I, p \in P, t \in T, t > t\_1 \tag{A20}$$

$$\begin{aligned} lID\_{j,p,t} &= ISP\_{j,p} + XI\_{j,p,t} + \sum\_{(i,m)\in Rautc\_{i,j,m}} XP\_{i,j,m,p,t} - \sum\_{(k,m)\in Rautc\_{j,k,m}} XS\_{j,k,m,p,t} \\ &\quad \forall j \in I, p \in P, t \in T, t = t\_1 \end{aligned} \tag{A21}$$

$$\begin{cases} ID\_{j,p,t} = ID\_{j,p,t-1} + XI\_{j,p,t} + \sum\_{(i,m)\in Raut\varepsilon\_{i,j,m}} XP\_{i,j,m,p,t} - \sum\_{(k,m)\in Raut\varepsilon\_{j,k,m}} XS\_{j,k,m,p,t} \\ \forall j \in \mathbb{J}, p \in P, t \in T, t > t\_1 \end{cases} \tag{A22}$$

$$\begin{aligned} lIM\_{k,p,t} &= lSP\_{k,p} + \sum\_{(i,m)\in Raut\_{i,k,m}} \text{XP}\_{i,k,m,p,t} + \sum\_{(j,m)\in Raut\_{j,k,m}} \text{XS}\_{j,k,m,p,t} + \text{Xl}l\_{k,p,t} - \text{X}\text{RP}\_{k,p,t} \\ & \forall k \in \mathcal{K}, p \in \mathcal{P}, t \in T, t = t\_1 \end{aligned} \tag{A23}$$

$$IM\_{k,p,t} = IM\_{k,p,t-1} + \sum\_{(i,m)\in Raut\_{i,k,m}} XP\_{i,k,m,p,t} + \sum\_{(j,m)\in Raut\_{j,k,m}} XS\_{j,k,m,p,t} + Xl\_{k,p,t} - XRP\_{k,p,t} \tag{A24}$$
 
$$\forall k \in \mathcal{K}, p \in P, t \in T, t > t\_1$$

$$XRP\_{k,p,t} \le DP\_{k,p,t} \quad \forall k \in K, p \in P, t \in T \tag{A25}$$

$$
\forall SSO\_i \times SOC\_i \le IO\_{i,t} \le SOC\_i \quad \forall i \in I, t \in T \tag{A26}
$$

$$\text{SSP}\_{i} \times \text{SC}\_{i,p} \le \text{IR}\_{i,p,t} \le \text{SC}\_{i,p} \quad \forall i \in I, p \in P, t \in T \tag{A27}$$

$$\text{SSP}\_{j} \times \text{SC}\_{j,p} \le \text{ID}\_{j,p}, \mathbf{t} \le \text{SC}\_{j,p} \quad \forall j \in \mathbf{J}, p \in \mathbf{P}, \mathbf{t} \in T \tag{A28}$$

$$\text{SSP}\_k \times \text{SC}\_{k,p} \le \text{IM}\_{k,p,t} \le \text{SC}\_{k,p} \quad \forall k \in \mathbb{K}, p \in P, t \in T \tag{A29}$$

$$\sum\_{l \in Raut \cdot c\_{l,l,m}} XP\_{i,l,m,p,t} \le ASPD\_{i,m,p} \quad \forall i \in I, m \in M, p \in P, t \in T \tag{A30}$$

$$\sum\_{k \in \text{Raut}: j, k, m} X S\_{j, k, m, p, t} \le A S S D\_{j, m, p} \quad \forall j \in \mathcal{J}, m \in M, p \in P, t \in T \tag{A31}$$

$$\sum\_{\mathbf{p}\in\mathbf{P}} XP\_{i,l,m,\mathbf{p},t} \le A C \mathbf{P} D\_{i,l,m} \quad \forall (i,l,m) \in \text{Raut}\_{i,l,m}, t \in T \tag{A32}$$

$$\sum\_{p \in P} X S\_{\vec{j},k,m,p,t} \le A \text{CSD}\_{\vec{j},k,m} \quad \forall \left(j,k,m\right) \in \text{Rute}\_{\vec{j},k,m,t} \tag{A33}$$

$$\sum\_{t \in (j,m) \in Raut \cdot \varepsilon\_{i,j,m}} XP\_{i,j,m,p,t} + XI\_{j,p,t} \le SC\_{j,p} \times TCM\_{j,p} \quad \forall j \in J, p \in P, t \in T \tag{A34}$$

$$\sum\_{(\mathbf{j},m)\in\text{Raut}\_{\mathbf{j},\mathbf{k},\mathbf{m}}} \text{XP}\_{i,\mathbf{k},\mathbf{m},\mathbf{p},\mathbf{t}} + \sum\_{(\mathbf{j},m)\in\text{Raut}\_{\mathbf{j},\mathbf{k},\mathbf{m}}} \text{XS}\_{j,\mathbf{k},\mathbf{m},\mathbf{p},\mathbf{t}} \le \text{SC}\_{\mathbf{k},\mathbf{p}} \times \text{TCM}\_{\mathbf{k},\mathbf{p}} \quad \forall \mathbf{k}\in\text{K}, \mathbf{p}\in P, \mathbf{t}\in T \quad \text{(A.35)}$$

$$XP\_{i,l,m,p,t} = 0 \quad \forall i \in I, l \in L, p \in P, t \in T, m = rand\_r \mu\_{i,j,m} > MTD \tag{A36}$$

$$XS\_{j,k,m,p,t} = 0 \quad \forall j \in I, k \in K, p \in P, t \in T, m = r \text{and}, \mu\_{j,k,m} > MTD \tag{A37}$$

#### **Table A2.** Model notation.


**Table A2.** *Cont.*


#### **Appendix C. Non-Adjustable Robust Optimization (NARO) Mathematical Formulation**

The non-adjustable robust formulation considers the same polyhedral budget uncertainty set of its equivalent ARO model, presented in the full paper. Hereinafter, the NARO model is introduced through its objective function (Equations (A38)–(A50)), network equations (Equations (A51)–(A77)) and model notation (Table A3). The presentation of the ARO model in the full paper must be consulted for details on the descriptions of the equations.

Objective function

$$\begin{array}{rl}\textbf{maximize} & profit\\hline\\ \int\_{k,p,t}^{+\text{demand}} \int\_{k,p,t}^{-\text{demand}} J\_{k,p,t}^{\text{demand}} & \int\_{k,p,t}^{\text{demand}} \\ \text{s.t.} & profit-\sum\_{i \in I, t \in T} MR\_{i,t} + \sum\_{j \in I, p \in P, t \in T} MD\_{j,p,t} + \sum\_{k \in K, p \in P, t \in T} MC\_{k,p,t} - \sum\_{i \in I, p \in P, t \in T} CE\_{i,p,t} \\ & -\sum\_{h \in H, p \in P, t \in T} CI\_{h,p,t} - \sum\_{i \in I, p \in P, t \in T} CD\_{i,p,t} - \sum\_{j \in I, p \in P, t \in T} SD\_{j,p,t} \\ & -\sum\_{i \in I, t \in T} CI\_{i,t} - \sum\_{i \in I, p \in P, t \in T} CI\_{k,p,t} - \sum\_{j \in I, p \in P, t \in T} CD\_{j,p,t} \\ & -\sum\_{k \in K, p \in P, t \in T} CI\_{k,p,t} - \sum\_{k \in K, p \in P, t \in T} CD\_{k,p,t} \le 0 \\ \end{array} \tag{A38}$$

$$MR\_{i,t} = \sum\_{p \in P} \left[ \mathbf{X} R\_{i,p,t} \times \left( \mathbf{P} P\_{a,p,t} - T \mathbf{N}\_{r,p} \right) \right] - \left[ \mathbf{X} \mathbf{CO}\_{i,t} \times \mathbf{PO}\_t \right] \forall i \in I, t \in T, a = r = r \text{-refner} \, y \tag{A39}$$

$$\begin{aligned} MD\_{j,p,t} &= \sum\_{(k,m)\in Raut\varepsilon\_{i,j,m}} \text{XS}\_{j,k,m,p,t} \times \left( PP\_{a,p,t} - TN\_{r,p} \right) - \sum\_{(i,m)\in Raut\varepsilon\_{i,j,m}} \text{XP}\_{i,j,m,p,t} \times PP\_{a\_2,p,t} \\ &\forall j \in \mathcal{J}, p \in \mathcal{P}, t \in T, a = r = \text{depot}, a\_2 = r \text{f}merry \end{aligned} \tag{A40}$$

$$\begin{array}{l} \text{MC}\_{k,p,t} \\ = \left[ \left( \text{X} \text{PP}\_{k,p,t} - \text{X} \text{LI}\_{k,p,t} \right) \times \left( \text{PP}\_{a\_1, p, t} - \text{TN}\_{r, p} \right) \right] \\ - \left[ \left( \sum\_{(j,m) \in \text{Raut}\_{i,k,m}} \text{XP}\_{i,k,m,p,t} \right) \times \text{PP}\_{a\_2, p,t} + \left( \sum\_{(j,m) \in \text{Raut}\_{j,k,m}} \text{XS}\_{j,k,m,p,t} \right) \right] \end{array} \tag{A41}$$

$$\forall k \in K, p \in P, t \in T, a\_1 = r = r \\ \text{tail}, a\_2 = r \\ \text{if} \\ \text{ir} \\ y, a\_3 = \\ \text{depth}$$
 
$$\text{CE}\_{i, p, t} = \text{XE}\_{i, p, t} \times \text{TN}\_{r, p} \quad \forall i \in I, p \in P, t \in T, r = \text{exp} \\ \text{rot} \tag{A42}$$

$$\text{CI}\_{\text{h},p,t} = \text{XI}\_{\text{h},p,t} \times \left( \text{PP}\_{a,p} + \text{TN}\_{r,p} \right) \quad \forall h \in H, p \in P, t \in T, a = r = \text{import} \tag{A43}$$

$$\text{CPD}\_{i,p,t} = \sum\_{(l,m)\in \text{Raut}\_{i,l,m}} \left( \text{XP}\_{i,l,m,p,t} \times \text{CT}\_{m,p} \times \text{Dist}\_{i,l,m} \right) \quad \forall i \in I, p \in P, t \in T \tag{A44}$$

$$\text{CSD}\_{j,p,t} = \sum\_{(k,m)\in\text{Raut}\_{j,k,m}} \left( \text{XS}\_{j,k,m,p,t} \times \text{CT}\_{m,p} \times \text{Dist}\_{j,k,m} \right) \quad \forall j \in I, p \in P, t \in T \tag{A45}$$

$$\text{CIO}\_{i,t} = \text{CKI} \times \text{IO}\_{i,t} \times \text{PO}\_{t} \quad \forall i \in I, t \in T, p = \text{oil} \tag{A46}$$

$$\text{CIR}\_{i,p,t} = \text{CII} \times \text{IR}\_{i,p,t} \times \text{PP}\_{a,p,t} \quad \forall i \in I, p \in \mathcal{P}, t \in T, a = r efficiency \tag{A47}$$

$$\text{CID}\_{j,p,t} = \text{CKI} \times \text{ID}\_{j,p,t} \times \text{PP}\_{a,p,t} \quad \forall j \in I, p \in P, t \in T, a = \text{depot} \tag{A48}$$

$$\text{CIM}\_{k,p,t} = \text{CKI} \times \text{IM}\_{k,p,t} \times \text{PP}\_{a,p,t} \quad \forall k \in K, p \in \mathcal{P}, t \in T, a = \text{retail} \tag{A49}$$

$$\text{CLID}\_{k,p,t} = \text{XLI}\_{k,p,t} \times \text{TN}\_{r,p} \quad \forall k \in \text{K}, p \in P, t \in T, r = \text{num} \, t \tag{A50}$$

Network equations

$$XOP\_{i,t} \le PC\_i \quad \forall i \in I, t \in T \tag{A51}$$

$$XR\_{i,p,t} = XOP\_{i,t} \times YF\_{i,p} \quad \forall i \in I, p \in P, t \in T \tag{A52}$$

$$XCO\_{i,t} \le DO\_{i,t} \quad \forall i \in I, t \in T \tag{A53}$$

$$IO\_{i,t} = ISO\_i + XCO\_{i,t} - XOP\_{i,t} \quad \forall i \in I, t \in T, t = t\_1 \tag{A54}$$

$$IO\_{i,t} = IO\_{i,t-1} + XCO\_{i,t} - XOP\_{i,t} \quad \forall i \in I, t \in T, t > t\_1 \tag{A55}$$

$$IR\_{i,p,t} = ISP\_{i,p} + XR\_{i,p,t} + XI\_{i,p,t} - \sum\_{(l,m)\in Raut\varepsilon\_{i,l,m}} XP\_{i,l,m,p,t} - XE\_{i,p,t} \quad \forall i \in I, p \in P, t \in T, t = t\_1 \tag{A56}$$

$$IR\_{i,p,t} = IR\_{i,p,t-1} + XR\_{i,p,t} + XI\_{i,p,t} - \sum\_{(l,m)\in Rautc\_{i,l,m}} XP\_{i,l,m,p,t} - XE\_{i,p,t} \quad \forall i \in I, p \in P\_r \\ t \in T\_r \\ t > t\_1 \tag{A57}$$

$$\begin{aligned} IID\_{j,p,t} &= ISP\_{j,p} + XI\_{j,p,t} + \sum\_{(i,m)\in Raut\_{i,j,m}} XP\_{i,j,m,p,t} - \sum\_{(k,m)\in Raut\_{j,k,m}} XS\_{j,k,m,p,t} \\ &\forall j \in I, p \in P, t \in T, t = t\_1 \end{aligned} \tag{A.58}$$

$$\begin{aligned} \text{ID}\_{j,p,t} = \text{ID}\_{j,p,t-1} + \text{XI}\_{j,p,t} + \sum\_{(i,m)\in\text{Raut}\_{i,j,m}} \text{XP}\_{i,j,m,p,t} - \sum\_{(k,m)\in\text{Raut}\_{j,k,m}} \text{XS}\_{j,k,m,p,t} \\ \forall j \in \text{J}\_{\prime}, p \in \text{P}, t \in T, t > t\_1 \end{aligned} \tag{A59}$$

$$\begin{aligned} I\text{-}IM\_{k,p,t} &= ISP\_{k,p} + \sum\_{(i,m)\in\text{Raut};k,m} \text{XP}\_{i,k,m,p,t} + \sum\_{(j,m)\in\text{Raut};j,m} \text{XS}\_{j,k,m,p,t} + \text{XL}\_{k,p,t} - \text{XRP}\_{k,p,t} \\ &\forall k \in \text{K}, p \in P, t \in T, t = t\_1 \end{aligned} \tag{A60}$$

$$IM\_{k,p,t} = IM\_{k,p,t-1} + \sum\_{(i,m)\in Raut\_{i,k,m}} XP\_{i,k,m,p,t} + \sum\_{(j,m)\in Raut\_{j,k,m}} XS\_{j,k,m,p,t} + Xl\_{k,p,t} - XRP\_{k,p,t} \tag{A61}$$
 
$$\forall k \in \mathbb{K}, p \in P, t \in T, t > t\_1$$

$$\begin{split} \text{XRP}\_{k,p,t} - \text{DP}\_{k,p,t} \text{Y}\_{\mathcal{g}} + \text{DP}\_{k,p,t}^{\text{max}} & \left( \boldsymbol{\xi}\_{k,p,t}^{+drmal} + \boldsymbol{\xi}\_{k,p,t}^{-drmal} \right) + \boldsymbol{\Gamma}\_{k,p}^{drmand} \boldsymbol{\eta}\_{k,p}^{drmand} \le 0 \\ \forall \boldsymbol{k} \in \mathcal{K}, p \in \mathcal{P}, t \in \mathcal{T}, \mathcal{g} &= |l| + 1 \end{split} \tag{A62}$$

$$\xi\_{k,pt}^{+,dmaxl} + \frac{\eta\_{k,p}^{drand}}{DP\_{k,p,t}^{max}} \ge +\mathcal{Y}\_{\mathcal{S}} \quad \forall i \in I, \mathcal{g} = |l| + 1, t \in T \tag{A63}$$

$$\xi\_{k,p,t}^{\text{demand}} + \frac{\eta\_{k,p}^{\text{demand}}}{DP\_{k,p,t}^{\text{max}}} \ge -Y\_{\mathcal{S}} \quad \forall i \in I, \mathcal{g} = |I| + 1, t \in T \tag{A64}$$

$$Y\_{\mathcal{S}} = \mathbf{1}\mathcal{g} = |\mathcal{l}| + \mathbf{1} \tag{A65}$$

$$
\forall SSO\_i \times \text{SCO}\_i \le IO\_{i,t} \le \text{SCO}\_i \quad \forall i \in I, t \in T \tag{A66}
$$

$$\text{SSP}\_{i} \times \text{SC}\_{i,p} \le \text{IR}\_{i,p,t} \le \text{SC}\_{i,p} \quad \forall i \in I, p \in P, t \in T \tag{A67}$$

$$\text{LSP}\_{j} \times \text{SC}\_{j,p} \le \text{ID}\_{j,p,t} \le \text{SC}\_{j,p} \quad \forall j \in \text{J}, p \in P, t \in T \tag{A68}$$

$$\text{L.S.P}\_k \times \text{SC}\_{k,p} \le \text{IM}\_{k,p,t} \le \text{SC}\_{k,p} \quad \forall k \in \text{K}, p \in P, t \in T \tag{A69}$$

$$\sum\_{l \in Rautc\_{i,l,m}} XP\_{i,l,m,p,t} \le ASPD\_{i,m,p} \quad \forall i \in I, m \in M, p \in P, t \in T \tag{A70}$$

$$\sum\_{k \in Raut \colon c\_{j,k,m}} X S\_{j,k,m,p,t} \le A S SD\_{j,m,p} \quad \forall j \in \mathcal{J}, m \in M, p \in P, t \in T \tag{A71}$$

$$\sum\_{p \in \mathcal{P}} XP\_{i,l,m,p,t} \le A \text{CPD}\_{i,l,m} \quad \forall (i,l,m) \in Rout \varepsilon\_{i,l,m}, t \in T \tag{A72}$$

$$\sum\_{p \in P} X S\_{\langle j,k,m,p,t \rangle} \le A \text{CSD}\_{\langle j,k,m \rangle} \quad \forall (j,k,m) \in \text{Rout}\_{\langle j,k,m \rangle} t \in T \tag{A73}$$

$$\sum\_{(i,m)\in Raut\_{\mathbb{G}\_{i,j,m}}} XP\_{i,j,m,p,t} + XI\_{j,p,t} \le SC\_{j,p} \times TCM\_{j,p} \quad \forall j \in \mathbb{J}, p \in P, t \in T \tag{A74}$$

$$\sum\_{i,j,m \in \text{Raut}\_{i,k,m}} X P\_{i,k,m,p,t} + \sum\_{(j,m) \in \text{Raut}\_{j,k,m}} X S\_{j,k,m,p,t} \le \text{SC}\_{k,p} \times T \text{CM}\_{k,p} \quad \forall k \in \text{K}, p \in P, t \in T \tag{A75}$$

$$XP\_{iJ,m,p,t} = 0 \quad \forall i \in I, l \in L, p \in P, t \in T, m = r \\ \\ \\ \\ a d\_t \\ \mu\_{i,j,m} > MTD \tag{A76}$$

$$XS\_{j,k,m,p,t} = 0 \quad \forall j \in \mathcal{I}, k \in \mathcal{K}, p \in P, t \in T, m = r \text{and}, \mu\_{j,k,m} > MTD \tag{A77}$$

#### **Table A3.** Model notation.


**Table A3.** *Cont.*


Robust parameters Γ*demand <sup>k</sup>*,*<sup>p</sup>* Budget parameter to adjust the robustness of product demand *DPmax <sup>k</sup>*,*p*,*<sup>t</sup>* Maximum variation in product demand for market *k* and product *p* at time point *t* Robust dual variables η*demand k*,*p* Dual variable associated with the establishment of the budget parameter of product demand ξ+*demand k*,*p*,*t* Quantify the sensitivity to positive deviation in product demand for market *k* and product *p* at time point *t* ξ−*demand k*,*p*,*t* Quantify the sensitivity to negative deviation in product demand for market *k* and product *p* at time point *t*

**Table A3.** *Cont.*

#### **Appendix D. Stochastic Mathematical Programming Formulation**

The two-stage stochastic programming (SP) model is formulated using node-variable formulation, where the decision variables of the optimization problem are associated with the nodes of the scenario tree. The objective function (Equations (A78)–(A92)), network equations (Equations (A93)–(A117)) and model notation (Table A4) of the SP model are presented below. The scenario tree to represent the evolution of the product demands is established considering an optimistic growth of 5% with 0.35 probability, a realistic and unchangeable occurrence of 0% with 0.50 probability and a pessimistic decrease of 10% with 0.15 probability, in accordance with the reasoning developed by Fernandes et al. [36]. For more details about the SP formulation, the interested reader is referred to Lima et al. [1].

Objective function

$$\underset{XCO}{\text{maximize}} - \sum\_{i \in I, t \in T} XCO\_{i,t}PO\_t + \mathcal{Q}(XCO) \tag{A78}$$

$$\text{s.t.}\,XCO\_{i,t} \le DO\_{i,t} \forall i \in I\_\prime\\ \text{t.t.}\tag{A79}$$

$$\begin{split} \mathcal{Q}(\text{XCO}) &= \underset{\text{XCO},\text{O}}{\text{maximize}} \sum\_{s \in \mathcal{S}} P\_{s} \Big( \sum\_{i \in I, t \in T} RR\_{i,t} + \sum\_{j \in I, p \in \mathcal{P}, t \in T} MD\_{j, p, t} + \sum\_{k \in K, p \in \mathcal{P}, t \in T} MC\_{k, p, t} \\ &- \sum\_{i \in I, p \in \mathcal{P}, t \in T} CE\_{i, p, t} - \sum\_{h \in H, p \in \mathcal{P}, t \in T} CI\_{h, p, t} - \sum\_{i \in I, p \in \mathcal{P}, t \in T} PD\_{i, p, t} \\ &- \sum\_{j \in I, p \in \mathcal{P}, t \in T} CSD\_{j, p, t} - \sum\_{i \in I, p \in \mathcal{P}} CIO\_{i, t} - \sum\_{i \in I, p \in \mathcal{P}, t \in T} CR\_{i, p, t} \end{split} \tag{A80}$$

$$\begin{aligned} & -\sum\_{\substack{j \in I, p \in P, t \in T}} \text{CID}\_{j, p, t} - \sum\_{\substack{k \in K, p \in P, t \in T}} \text{CIM}\_{k, p, t} - \sum\_{\substack{k \in K, p \in P, t \in T}} \text{CID}\_{k, p, t} \right) \\ & \qquad \text{RR}\_{i, t} = \sum\_{p \in P} \left[ \text{XR}\_{i, p, t, s} \times \left( \text{PP}\_{a, p, t} - \text{TN}\_{l, p} \right) \right] \quad \forall i \in I\_{\star} \ (t, s) \in \text{TS}, a = r = r \neq \text{refner} \, y \end{aligned} \tag{A81}$$

$$\begin{aligned} \text{MD}\_{j,p,t,s} &= \left[\sum\_{(\mathbf{i},\mathbf{m})\in\text{Raut}\_{j,\mathbf{m}}} \text{XS}\_{j,\mathbf{k},\mathbf{m},p,t,s} \times \left(\text{PP}\_{\mathbf{a},p,t} - \text{TN}\_{r,p}\right)\right] - \left[\sum\_{(\mathbf{i},\mathbf{m})\in\text{Raut}\_{i,j,\mathbf{m}}} \text{XP}\_{i,j,\mathbf{m},p,t,s} \times \text{PP}\_{\mathbf{a},p,t}\right] \\ &\forall j \in \text{J}, p \in \text{P}\_{\text{\textquotedblleft}t}(t,s) \in \text{TS}, a = r = \text{depot}, a\_2 = r \text{f} \text{iner}ry \end{aligned} \tag{A82}$$

$$\begin{array}{c} \text{MC}\_{k,p,t,s} = \left[ \left( \text{XRP}\_{k,p,t,s} - \text{XLI}\_{k,p,t,s} \right) \times \left( \text{PP}\_{a\_1,p,t} - \text{TN}\_{r,p} \right) \right] - \left[ \left( \sum\_{(i,m) \in \text{Baut}\_{i,k,m}} \text{XP}\_{i,k,m,p,t,s} \right) \times \text{PP}\_{a\_1,p,t,s} \right] \\ \quad \times \text{PP}\_{a\_2,p,t} + \left( \sum\_{(j,m) \in \text{Baut}\_{j,k,m}} \text{XS}\_{j,k,m,p,t,s} \right) \times \text{PP}\_{a\_3,p,t} \end{array} \tag{A83}$$

$$\forall k \in K, p \in P, \stackrel{\frown}{(t,s)} \in TS, a\_1 = r = \stackrel{\frown \smile}{\text{retial}}, a\_2 = \stackrel{\frown \smile}{\text{refinary}}, a\_3 = \stackrel{\frown}{\text{depot}}$$

$$\text{CE}\_{i,p,t,s} = \text{XE}\_{i,p,t,s} \times \text{TN}\_{r,p} \quad \forall i \in I, p \in P, (t,s) \in \text{TS}, r = \text{export} \tag{A84}$$

$$\text{CI}\_{h,p,t,s} = \text{XI}\_{h,p,t,s} \times \left( \text{PP}\_{a,p,s} + \text{TN}\_{r,p} \right) \quad \forall h \in H, p \in \mathcal{P}, (t,s) \in \text{TS}, a = r = \text{import} \tag{A85}$$

$$\text{CPD}\_{i,p,t,s} = \sum\_{(l,m)\in Raut\_{i,l,m}} \left( \text{XP}\_{i,l,m,p,t,s} \times \text{CT}\_{m,p} \times \text{Dist}\_{i,l,m} \right) \quad \forall i \in I, p \in P, (t,s) \in \text{TS} \tag{A86}$$

$$\text{CSD}\_{j,p,t,s} = \sum\_{(k,m)\in \text{Raut}\_{j,k,m}} \left( \text{XS}\_{j,k,m,p,t,s} \times \text{CT}\_{m,p} \times \text{Dist}\_{j,k,m} \right) \quad \forall j \in l\_r \\ p \in P, (t,s) \in \text{TS} \tag{A87}$$

$$\text{CIO}\_{i,t,s} = \text{CKI} \times \text{IO}\_{i,t,s} \times \text{PO}\_{a,s} \quad \forall i \in I, (t,s) \in \text{TS}, a = \text{procurement} \tag{A88}$$

$$\text{CIR}\_{i,p,t,s} = \text{CKI} \times \text{IR}\_{i,p,t,s} \times \text{PP}\_{a,p,t} \quad \forall i \in I, p \in \mathcal{P}, (t,s) \in \text{TS}, a = referary \tag{A89}$$

$$\text{CID}\_{j,p,t,s} = \text{CKI} \times \text{ID}\_{j,p,t,s} \times \text{PP}\_{a,p,t} \quad \forall j \in \text{J}, p \in P, (t,s) \in \text{TS}, a = \text{depot} \tag{A90}$$

$$\text{CIM}\_{k,p,t,s} = \text{CKI} \times \text{IM}\_{k,p,t,s} \times \text{PP}\_{a,p,t} \quad \forall k \in \mathcal{K}, p \in \mathcal{P}, (t,s) \in \text{TS}, a = \text{retail} \tag{A91}$$

$$\text{C}\,\text{C}\,\text{U}\,\text{D}\_{k,p,t,\text{s}} = \text{X}\,\text{U}\_{k,p,t,\text{s}} \times \text{TN}\_{r,p} \quad \forall k \in \text{K}, p \in P, (t, \text{s}) \in \text{TS}, r = \text{unment} \tag{A92}$$

Network equations

$$XOP\_{i,t,s} \le PC\_i \quad \forall i \in I, (t,s) \in TS \tag{A93}$$

$$XR\_{i,p,t,s} = XOP\_{i,t,s} \times YF\_{i,p} \quad \forall i \in I, p \in P, (t,s) \in TS \tag{A94}$$

$$IO\_{i,t,s} = ISO\_i + DO\_{i,t} - XOP\_{i,t,s} \quad \forall i \in I\_\prime(t,s) \in TS, t = t\_\mathbf{1} \tag{A95}$$

$$IIO\_{i,t,s} = IO\_{i,t-1,s} + DO\_{i,t} - XOP\_{i,t,s} \quad \forall i \in I\_\star(t,s) \in TS, t > t\_1, s \in \overline{S} \tag{A96}$$

$$\begin{aligned} IR\_{i,p,t,s} &= ISP\_{i,p} + XR\_{i,p,t,s} + XI\_{i,p,t,s} - \sum\_{(l,m)\in Raut\varepsilon\_{i,l,m}} XP\_{i,l,m,p,t,s} - XE\_{i,p,t,s} \\ &\quad \forall i \in I, p \in P, (t,s) \in TS, t = t\_1 \end{aligned} \tag{A97}$$

$$\begin{aligned} IR\_{i,p,t,s} &= IR\_{i,p,t-1,\overline{s}} + XR\_{i,p,t,s} + XI\_{i,p,t,s} - \sum\_{(l,m)\in Raut\varepsilon\_{i,l,m}} XP\_{i,l,m,p,t,s} - XE\_{i,p,t,s} \\ &\quad \forall i \in I, p \in P\_{\prime} \ (t,s) \in TS, t > t\_1, s \in \overline{S\overline{S}} \end{aligned} \tag{A98}$$

$$\begin{aligned} \text{i } ID\_{j,p,t,s} = \text{I}SP\_{j,p} + \text{XI}\_{j,p,t,s} + \sum\_{(i,m)\in Rautte\_{i,j,m}} \text{XP}\_{i,j,m,p,t,s} - \sum\_{(k,m)\in Rautte\_{j,k,m}} \text{XS}\_{j,k,m,p,t,s} \\ \forall j \in \text{J}, p \in \text{P}, (t,s) \in \text{TS}, t = t\_1 \end{aligned} \tag{A99}$$

$$\begin{aligned} \text{IDP}\_{j,p,t,s} &= \text{ID}\_{j,p,t-1,\overline{s}} + \text{XI}\_{j,p,t,s} + \sum\_{(i,m)\in\text{Raut}\_{t,j,m}} \text{XP}\_{i,j,m,p,t,s} - \sum\_{(k,m)\in\text{Raut}\_{j,k,m}} \text{XS}\_{j,k,m,p,t,s} \\ &\forall j \in \text{J}, p \in \text{P}\_{\text{\textquotedblleft}t}(t,s) \in \text{TS}, t > t\_{1\prime}, s \in \overline{S\text{S}} \end{aligned} \tag{A100}$$

$$\begin{aligned} I\text{IM}\_{k,p,l,s} &= I\text{SP}\_{k,p} + \sum\_{(i,m)\in\text{Rautle}\_{i,k,m}} \text{XP}\_{(k,m,p,l,s} + \sum\_{(i,m)\in\text{Rautle}\_{j,k,m}} \text{XS}\_{[k,m,p,l,s} + \text{XL}\_{k,p,l,s} - \text{XRP}\_{k,p,l,s} \\ & \forall k \in \text{K}, p \in \text{P}, (t,s) \in \text{TS}, t = t\_1 \end{aligned} \tag{A101}$$

$$\begin{aligned} lM\_{k,p,l,s} = l\mathbb{C}\_{k,p,l=1,7} + \sum\_{(i,m)\in\text{Raut}\_{i,k,m}} \text{XP}\_{i,k,m,p,l,s} + \sum\_{(i,m)\in\text{Raut}\_{i,k,m}} \text{XS}\_{j,k,m,p,l,s} + \text{XL}\_{k,p,l,s} - \text{XR}\_{k,p,l,s} \\ \forall k \in \text{K}, p \in \text{P}\_{\text{\textdegree}}(t,s) \in \text{TS}, t > t\_1, s \in \overline{\text{SS}} \end{aligned} \tag{A102}$$

$$DPR\_{k,p,t,s} = DP\_{k,p} \forall k \in \mathbb{K}, p \in P, (t,s) \in TS, t = t\_1 \tag{A103}$$

$$DPR\_{k,p,t,s} = DPR\_{k,p,t-1,5} \chi^p \psi^s \quad \forall k \in K, p \in P, (t,s) \in TS, t > t\_1, s \in S\overline{S} \tag{A104}$$

$$\text{VAR}\_{k,p,t,s} \le DPR\_{k,p,t,s} \quad \forall k \in K, p \in P\_\prime\left(t, s\right) \in TS \tag{A105}$$

$$\text{LSSO}\_{\text{i}} \times \text{SCO}\_{\text{i}} \le IO\_{\text{i},\text{s}} \le \text{SCO}\_{\text{i}} \quad \forall \text{i} \in I\_{\text{\textdegree}} \text{(t,s)} \in TS \tag{A106}$$

$$\text{LSSP}\_{i} \times \text{SC}\_{i,p} \le \text{IR}\_{i,p,t,s} \le \text{SC}\_{i,p} \quad \forall i \in I, p \in P, (t,s) \in \text{TS} \tag{A107}$$

$$\text{SSP}\_{j} \times \text{SC}\_{j,p} \le \text{ID}\_{j,p,t,s} \le \text{SC}\_{j,p} \quad \forall j \in \text{J}, p \in P, (t,s) \in \text{TS} \tag{A108}$$

$$\text{SSP}\_k \times \text{SC}\_{k,p} \le \text{IM}\_{k,p,t,s} \le \text{SC}\_{k,p} \quad \forall k \in \mathcal{K}, p \in P\_\prime\left(t, s\right) \in \text{TS} \tag{A109}$$

$$\sum\_{l \in \text{Ront}\_{i,l,m}} \text{XP}\_{i,l,m,p,t,s} \le A \text{SPD}\_{i,m,p} \quad \forall i \in I, m \in M, p \in P\_{\text{\textquotedblleft}}(t,s) \in TS \tag{A110}$$

$$\sum\_{k \in Raut \atop j,k,m} X S\_{j,k,m,p,t,s} \le ASSD\_{j,m,p} \quad \forall j \in \mathcal{J}, m \in M, p \in P\_\prime\left(t,s\right) \in TS \tag{A111}$$

$$\sum\_{p \in P} XP\_{i,l,m,p,t,s} \le ACD\_{i,l,m} \quad \forall (i,l,m) \in Rautte\_{i,l,m}(t,s) \in TS \tag{A112}$$

$$\sum\_{p \in P} X S\_{j,k,m,p,t,s} \le A \text{CSD}\_{j,k,m} \quad \forall (j,k,m) \in \text{Raut}\_{j,k,m}(t,s) \in T \text{S} \tag{A113}$$

$$\sum\_{(i,m)\in Raut\epsilon\_{i,m}} X P\_{i,j,m,p,t,s} + X I\_{j,p,t,s} \le S C\_{j,p} \times T C M\_{j,p} \quad \forall j \in I, p \in P, (t,s) \in T S \tag{A114}$$

$$\sum\_{(j,m)\in Raut\varepsilon\_{j,k,m}} X P\_{i,k,m,p,t,s} + \sum\_{(j,m)\in Raut\varepsilon\_{j,k,m}} X S\_{j,k,m,p,t,s} \le S C\_{k,p} \times T C M\_{k,p} \quad \forall k \in K, p \in P\_{\prime} \ (t,s) \in T S \quad (A115)$$

$$XP\_{i,l,m,p,t,s} = 0 \quad \forall i \in I, l \in L, p \in P, (t,s) \in TS, m = rand, Dist\_{i,j,m} > MTD \tag{A116}$$

$$XS\_{j,k,m,p,t,s} = 0 \quad \forall j \in f\_r \newline k \in K\_r \newline p \in P\_r(t,s) \in TS\_r \newline m = r \newline a \newline d, Dist\_{j,k,m} > MTD \newline \tag{A117}$$

#### **Table A4.** Model notation.


88

**Table A4.** *Cont.*


#### **Appendix E. Considerations about a Typical Polyhedral Budget Uncertainty Set**

In this part, we demonstrate how to enumerate all the possible vertices *v* ∈ *V* in the budget uncertainty set U for a generic pair of location *k* and product *p*, and a specific budget of uncertainty Γ*k*,*<sup>p</sup>* over a time horizon *T*. Consider the Equations (A118) and (A119) below:

$$\left| \Delta DP\_{k,p,t} \right| \le \Delta DP\_{k,p,t'}^{\max} \,\forall k \in K, p \in P, t \in T \tag{A118}$$

$$\sum\_{t} \frac{|\Delta DP\_{k,p,t}|}{\Delta DP\_{k,p,t}^{\max}} \le \Gamma\_{k,p\prime} \forall k \in \mathcal{K}, p \in \mathcal{P} \tag{A119}$$

Equation (A118) determines symmetrical intervals for the deviation of the product demand from the nominal value, while the total deviation across all time points are limited by the budget of uncertainty in Equation (A119). We can omit the indices *k* and *p* in Equation (A120) once they could refer to any pair of location and product. Consider a budget of uncertainty Γ*k*,*<sup>p</sup>* = 1 and a time horizon covering two time points as follows:

$$\frac{|\Delta DP\_1|}{\Delta DP\_1^{\text{max}}} + \frac{|\Delta DP\_2|}{\Delta DP\_2^{\text{max}}} \le 1\tag{A120}$$

Equation (A120) ensures that if the deviation of product demand at time point *t* = 1 is at the lower or upper bound of the range defined by Equation (A118), i.e., <sup>|</sup>Δ*DP*1<sup>|</sup> = <sup>Δ</sup>*DPmax* <sup>1</sup> , the deviation of the product demand at time point *t* = 2 will be necessarily zero, |Δ*DP*2| = 0. Inversely, when the deviation of product demand at time point *t* = 2 is at the lower or upper bound of the range defined by Equation (A118), i.e., <sup>|</sup>Δ*DP*2<sup>|</sup> = <sup>Δ</sup>*DPmax* <sup>2</sup> , the deviation of the product demand at time point *t* = 1 will be certainly zero, |Δ*DP*1| = 0. In such a way, we have enumerated all the possible scenarios *v* ∈ *V* in <sup>U</sup>, i.e., four vertices, where *<sup>V</sup>* = (Δ*DP*1, 0);(−Δ*DP*1, 0);(0,−Δ*DP*2);(0, Δ*DP*2) . Note that the size

of the set of vertices *V* depends on the value of the budget of uncertainty Γ*k*,*p*, which takes values in the range [0; 2]. When Γ*k*,*<sup>p</sup>* = 1.5, there are eight vertices within the set *V*, as shown below:

$$V = \left\{ \begin{array}{c} \left( \Delta DP\_{1\prime} - \frac{\Delta DP\_{2}}{2} \right); \left( \Delta DP\_{1\prime} \frac{\Delta DP\_{2}}{2} \right); \left( \frac{\Delta DP\_{1}}{2}, -\Delta DP\_{2} \right); \left( \frac{\Delta DP\_{1}}{2}, \Delta DP\_{2} \right); \\\left( -\Delta DP\_{1\prime} - \frac{\Delta DP\_{2}}{2} \right); \left( -\Delta DP\_{1\prime} \frac{\Delta DP\_{2}}{2} \right); \left( -\frac{\Delta DP\_{1}}{2}, -\Delta DP\_{2} \right); \left( -\frac{\Delta DP\_{1}}{2}, \Delta DP\_{2} \right) \end{array} \right\}$$

We can generalize that when Γ*k*,*<sup>p</sup>* = 0, the uncertainty set U has a just one vertex, corresponding to the nominal deterministic case. As Γ*k*,*<sup>p</sup>* increases, the size of the uncertainty set U enlarges. As shown before, when Γ*k*,*<sup>p</sup>* takes any value in the interval [0.01; 1], the polyhedron will have four vertices, while if Γ*k*,*<sup>p</sup>* assumes any value within the interval [1.01; 1, 99], the polyhedron will have eight vertices. On the other hand, when Γ*k*,*<sup>p</sup>* = 2, the polyhedron will have four vertices again.

It is important to highlight that this is just a generic illustration to show how to enumerate the vertices of a budget uncertainty set, and it was not used in the case study shown in the full paper. On the other hand, such instance can easily be extended to include a longer time horizon, such that the vertices of more complex polyhedral uncertainty sets can be determined.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Modern Modeling Paradigms Using Generalized Disjunctive Programming †**

#### **Qi Chen and Ignacio Grossmann \***

Center for Advanced Process Decision-Making, Carnegie Mellon University, Pittsburgh, PA 15213, USA; qichen@andrew.cmu.edu


Received: 14 October 2019; Accepted: 6 November 2019; Published: 10 November 2019

**Abstract:** Models involving decision variables in both discrete and continuous domain spaces are prevalent in process design. Generalized Disjunctive Programming (GDP) has emerged as a modeling framework to explicitly represent the relationship between algebraic descriptions and the logical structure of a design problem. However, fewer formulation examples exist for GDP compared to the traditional Mixed-Integer Nonlinear Programming (MINLP) modeling approach. In this paper, we propose the use of GDP as a modeling tool to organize model variants that arise due to characterization of different sections of an end-to-end process at different detail levels. We present an illustrative case study to demonstrate GDP usage for the generation of model variants catered to process synthesis integrated with purchasing and sales decisions in a techno-economic analysis. We also show how this GDP model can be used as part of a hierarchical decomposition scheme. These examples demonstrate how GDP can serve as a useful model abstraction layer for simplifying model development and upkeep, in addition to its traditional usage as a platform for advanced solution strategies.

**Keywords:** process design; process modeling; mathematical programming; MINLP; generalized disjunctive programming

#### **1. Introduction**

Mathematical programming is a powerful tool for process design and optimization, allowing the modeler to consider both continuous and discrete decisions. In process design, discrete degrees of freedom often determine topological structure (selection/activation/ordering of nodes and edges) while continuous variables determine system states such as flow rates or qualities. In the general case, these process design problems can involve nonlinear variable relationships and are addressed as Mixed-Integer Nonlinear Programming (MINLP) problems.

$$\begin{aligned} \min z &= f(x, y) \\ \text{s.t.} &\quad g(x, y) \le 0 \\ &\quad h(x, y) = 0 \\ &\quad \mathbf{x} \in X \subseteq \mathbb{R}^n \\ &\quad y \in Y \subseteq \mathbb{Z}^m \end{aligned} \tag{\text{MINL.P}}$$

The general form for these optimization models is given in (MINLP). An objective function *f*(*x*, *y*) is minimized by selecting values for continuous variables *x* and integer variables *y*, subject to satisfying inequality constraints *g*(*x*, *y*) ≤ 0 and equality constraints *h*(*x*, *y*) = 0. In processes, the continuous

variables usually represent flows, pressures, and temperatures. The integer variables are commonly 0–1 variables for the selection of units, but can also represent the number of units. The inequalities usually describe process and equipment limitations while equality constraints describe physical property relationships. An abundant literature exists for the formulation and solution of MINLP models [1–4]. Particularly in chemical engineering, many models are now formulated using algebraic relationships. Such equation-oriented modeling is becoming progressively more common [5], with differential equations used to describe temporal and spatial dynamics. For process design problems, postulation of alternatives is also an important consideration, with several approaches described in literature [6–15]. However, even with the growth of more complex models, there has been limited emphasis on the link between algebraic relationships and model logic.

Generalized Disjunctive Programming (GDP) represents one effort to systematize the relationship between algebraic relations and logical clauses [2,16,17], in pursuit of a framework to simplify both model formulation and solution of the eventual mathematical programming problem. GDP can be seen as the extension of theoretical work in disjunctive programming from the operations research community [18,19] to formulations involving nonlinear algebraic relationships. GDP gives the modeler a mathematical framework to express high-level logical statements without needing to immediately translate them into algebraic form. The general form for GDP optimization models is given in (GDP).

$$\begin{aligned} \min obj &= f(x, z) \\ s.t. &\quad g(x, z) \le 0 \\ &\quad \bigvee\_{i \in D\_k} \left[ \begin{array}{c} \mathbf{Y}\_{ik} \\ r\_{ik}(x, z) \le 0 \end{array} \right] \quad \forall k \in K \\ &\quad \bigvee\_{i \in D\_k} Y\_{ik} \quad \forall k \in K \\ &\quad \Omega(Y) = True \\ &\quad x \in X \subseteq \mathbb{R}^n \\ &\quad Y\_{ik} \in \{True, False \} \quad \forall i \in D\_k, \forall k \in K \\ &\quad z \in Z \subseteq \mathbb{Z}^n \end{aligned} \tag{GDPs}$$

As with the MINLP formulation, an objective function *f*(*x*, *z*) is minimized. Continuous decisions variables are still represented by *x*, but Boolean variables *Y* now describe selection among discrete alternatives. Remaining integer variables are denoted by *z*. This is preferable, as the conditional constraints *rik*(*x*, *z*) ≤ 0 corresponding to selection of alternative *Yik* can be grouped together and separated from the globally valid constraints *g*(*x*, *z*) ≤ 0 that must hold true for any selection of alternatives. Note that equality constraints are implicitly captured in (GDP) through the use of two inequality constraints. We term these groupings of a Boolean indicator variable with relevant conditional constraints a "disjunct", as they each constitute one term of a disjunction ∨ (logical "OR" relationship). Next, we state that for each disjunction *k* ∈ *K*, exactly one of the disjuncts *i* ∈ *Dk* will be selected, a generalization of the logical XOR ∨. Finally, GDP also allows the explicit specification of logical propositions Ω(*Y*) = *True* to describe logical relationships between selection of the discrete alternatives. These logical propositions are key to the modeling strategies addressed later in this work.

GDP offers two major advantages over the traditional MINLP modeling approach. First, it facilitates more intuitive modeling of process decision-making by allowing explicit specification of logical relationships [20]. The grouping of related constraints in disjuncts also helps to keep GDP models more organized. Second, by exploiting explicit logical structure provided by GDP models [21], advanced solution algorithms can reap benefits in convergence speed and robustness [22–24]. In this work, we focus on the modeling implications of GDP use.

While long-time practitioners of MINLP modeling approaches may sometimes find logical propositions too verbose, we contend that explicit logic is more readable and better preserves a modeler's original intent. Take for example the logical statement in Equation (1), which may correspond to the following process specification: if we purchase the cheap feed (*Y*<sup>1</sup> = *True*), then we will need to install a pretreater (*Y*<sup>2</sup> = *True*) or install an additional separation unit (*Y*<sup>3</sup> = *True*). That is, *Y*<sup>1</sup> implies *Y*<sup>2</sup> or *Y*3.

$$\mathcal{Y}\_1 \implies (\mathcal{Y}\_2 \lor \mathcal{Y}\_3) \tag{1}$$

This resolves to the equivalent algebraic constraint in Equation (2), using binary variables *y* [20,25].

$$-y\_1 + y\_2 + y\_3 \ge 0\tag{2}$$

Both descriptions are valid, but the logical statement is self-documenting and is much clearer to a new modeler. Expert MINLP modelers are accustomed to automatically preprocessing their formulation to encode relevant problem logic in the algebraic constraints. With GDP, this mental overhead and the associated potential for human error is eliminated. As a result, GDP models promise to be easier to develop and maintain.

Compounding this effect is the fact that mathematical programming is frequently an analysis tool in a larger process design procedure [26]. Business needs and customer expectations are not always initially stated in a form amenable to formulation as an algebraic objective or constraint [27]. A process design problem therefore involves several iterations of reassessing assumptions and adjusting constraints to match new business needs or revised customer expectations. This means that neither the data nor the structure for a model formulation can be regarded as static throughout the design workflow. As a result, the optimization model may be rewritten or revised several times in the course of tackling a single process design problem [28]. Considerations such as environmental/social impacts, safety, and operability may also be added to the model as constraints or secondary objectives. Moreover, multiple versions of the same model are often necessary to trade-off process detail versus model tractability for different process sections, increasing the number of model formulations that must be developed and maintained. By separating model logical structure from the underlying algebraic descriptions, GDP reduces the work necessary to revise a model. In doing so, it aims to advance the state-of-the-art in process design [26,29–32]. Later, we also show how GDP can help manage model variants.

GDP can be viewed in a broader context as a logic-based model abstraction layer, facilitating intuitive expression of the discrete decision spaces. These abstractions are a necessary response to complexity [33], to make mathematical programming capabilities accessible to a broader range of process modelers. The ubiquity of commercial chemical process simulators is attributable to their ability to abstract large-scale mathematical computations from chemical engineering decision-making. They provide a drag-and-drop interface to assemble a process structure and a drop-down menu to pick among standard physical property packages. Similar efforts to provide high-level modeling capabilities underpinned by mathematical programming—such as Egret [34] for power systems design, ICAS [35] for process and product design, MIPSYN [36] for process synthesis, and IDAES [37] for advanced energy systems design—can benefit from the logical abstraction provided by GDP. Note that modeling in GDP does not preclude the use of MINLP solution methods. Instead, it can offer a systematic yet flexible approach to generate the appropriate MINLP formulation via reformulation. For instance, Castro and Grossmann [38] derive several traditional scheduling formulations using standardized reformulations of GDP models.

The two most popular ways to reformulate a GDP model as an MINLP model are the Big-M (BM) [16] and Hull Reformulation (HR) [39] methods. BM and HR trade-off problem size (number of variables and constraints) and the tightness (quality) of the continuous relaxation. Other reformulations are also possible [40–43] with different tradeoffs in problem size, relaxation tightness, and computational cost to generate the reformulation. In general, multiple valid MINLP formulations exist to describe the same problem logic, and there exists no general way to determine a priori the most tractable formulation [17]. Direct formulation as an MINLP requires the modeler to

commit to a single formulation approach, while GDP allows multiple algebraic formulations to be systematically generated from a single logical description, so that the most advantageous variant may be utilized.

Despite theoretical progress, lack of computational tools to support GDP modeling has hindered its adoption in both academic and industrial settings. The GAMS algebraic modeling language [44] provides support for GDP models through its Extended Mathematical Programming syntax, with the ability to generate the BM and HR reformulations. Prior to version 23.7, GAMS also supported solution of GDP models using the LOGMIP 1.0 solver [45]. However, as a closed-source commercial platform, academic interest has been limited. More recently, Pyomo.GDP [46] has emerged as an open-source ecosystem for GDP modeling and development, built on top of the Pyomo algebraic modeling language [47] in Python. As an open-source platform, it has been able to incorporate recent innovations in reformulation strategies [48] and logic-based solution algorithms [22]. Powerful options now exist for formulating and solving GDP models. However, compared to the MINLP literature, relatively few formulation examples exist for GDP. This paper aims to address that gap.

In this paper, we focus on GDP as a modeling tool to manage model variants. We demonstrate its use for two modeling use cases: (1) end-to-end analysis with focus on various portions of the overall process, and (2) a single solution scheme involving use of models at different detail levels. In Section 2, we describe these use cases and their modeling challenges in detail. In Section 3, we discuss application of these techniques on an illustrative example and the resulting implications. We present concluding remarks in Section 4.

#### **2. Problem Statement**

We examine two scenarios in which GDP is useful as a model management platform.


In the first case, complex value chains often result in a process being subdivided among major process sections, assigned to different modeling teams. Each of these modeling teams may develop specialized, detailed models that describe decision points and specifications relevant to their portion of the overall process. However, as sequential optimizations of process sections may yield a suboptimal overall result, coordination is necessary. Each section therefore needs to model the secondary impact of their decisions on the rest of the process. Ideally, the detailed models for each section could simply be linked together to produce a single optimization formulation. However, this formulation is often intractably large. Therefore, a less-detailed surrogate is often employed to model nonfocal portions of the process.

Consider an illustrative chemical production process consisting of three sections: procurement, production, and sales, displayed in Figure 1. We examine this process in more detail in Section 3. Multiple discrete options are available as decision variables for each section, denoted by the Boolean variables *Yi*, *Yj*, and *Yk* for each section, respectively. For procurement, these discrete decisions may describe selection among several available supply contracts from various vendors. For sales, there might exist several sales opportunities corresponding to different customers. In production, selection among various production modes and capital purchases are often key decisions. For our illustrative process, *Y*<sup>12</sup> = *True* describes the selection of the second procurement contract, for example.

In typical industrial practice, models at varying detail levels are often developed separately from each other, with interoperability suffering as a result. By making use of GDP, the choice of modeling detail can be integrated within a single framework, allowing related models to be developed in proximity to each other. Note that care should still be taken to define appropriate interconnections between the process sections so that relevant phenomena can be described (e.g., time dependence).

**Figure 1.** Example process illustrating the embedding of multiple detail levels within discrete options for each process section.

In the second case, a modeler may wish to use approximate models in a preliminary analysis to identify promising candidates among numerous alternatives, then drill down into progressively more detailed representations for remaining alternatives. For chemical processes, making assumptions that restrict temperatures, transport phenomena, or thermodynamic complexity can greatly simplify a model representation. However, variations in degrees of freedom and relevant physical phenomena should be subsequently revisited. In specialized simulation software, provisions for changing the thermodynamic assumptions of a chemical process are commonplace. However, the ability to do so is less common in equation-oriented optimization frameworks. Instead, a new model must frequently be developed at the desired complexity level. With GDP, the imposition or relaxation of these assumptions may be made by setting the value of Boolean variables. For a broader perspective, by imposing the implication that selection of unit *<sup>u</sup>* requires its modeling at a low detail level, *Yu* =⇒ *<sup>Y</sup>low <sup>u</sup>* , we can restrict consideration to the approximate models. Conversely, for a deeper perspective, we can consider only more rigorously modeled alternatives by imposing the use of high detail models *Yu* <sup>=</sup><sup>⇒</sup> *<sup>Y</sup>high <sup>u</sup>* . Note that at high fidelity, some alternatives may themselves involve selection among a discrete decision space. For example, selection of a distillation column in a chemical process may involve deciding on the number of trays. This discrete decision could itself be treated at various levels of modeling fidelity. Fortunately, with GDP, this decision can simply be nested within the higher-level selection as a nested disjunction [49].

#### **3. Case Study**

To demonstrate the principles of GDP as a tool for model management, we propose an illustrative end-to-end methanol synthesis process example adapted from literature [22]. As previously stated, the process consists of three sections: procurement, production, and sales (see Figure 2). Procurement must source the syngas from one of two different vendors, each of which offers a different purity-dependent cost schedule per unit feed. Production must select the optimal equipment configuration and operating conditions (temperatures, pressures, flows, and compositions) to convert syngas to methanol. This includes the discrete decision between single and two-stage compression for both the feed and recycle streams, as well as the choice between a higher-conversion, higher-cost reactor and a cheaper variant. Sales then contracts with one of two different customers, who are willing to pay a unit

price that depends on the product purity. Each of these sections may be modeled at a high, medium, and low level of detail. The production section superstructure appears in Figure 3, with unit 9 as the cheap reactor alternative with low conversion and unit 10 denoting the expensive reactor with high conversion. Amendments to the literature methanol process synthesis model are presented in Appendix C.

**Figure 2.** Simplified process diagram for the illustrative example. Three sections exist: procurement, production, and sales. In each section, decisions must be made in both discrete and continuous variable domains.

**Figure 3.** Methanol process flowsheet superstructure, adapted from [22], showing stream numbers in blue.

Below, we present a generic formulation for an end-to-end techno-economic analysis model and demonstrate how it may be adapted to our methanol synthesis example. In the following subsections, we lay out the relevant sets, parameters, decision variables, and functions before presenting the formulation. We then discuss manipulations of the formulation for both use cases from Section 2.

#### *3.1. Sets*

The chemical process involves a set of components *C*, which include the relevant raw materials, reaction intermediates, inerts, and products. For the methanol process, we have feed components *H*<sup>2</sup> and *CO*, inert *CH*4, and product *CH*3*OH*. The process is subdivided into three major process sections *I*: procurement, production, and sales. Of all possible process section alternatives *J*, a subset *Ji* ⊆ *J* is available for selection for each section *i* ∈ *I*. In the production section of the methanol process, two of these alternatives are single-stage feed compression, and two-stage feed compression. The process also involves a set of streams *K* that describe flows of material between process sections. Finally, for each process alternative, a set of modeling detail levels *l* ∈ *L* are available. We consider in our example three detail levels: "low", "medium", and "high."


#### *3.2. Variables*

The decision variables include characterization of each process stream *k* ∈ *K*: total molar flowrate *Fk* and component molar flowrate *fkc* for each component *c* ∈ *C*, as well as the stream temperature *Tk* and pressure *Pk*. A profit or cost (negative) contribution from each section *i* ∈ *I* is given by *zi*. This, in turn, may be influenced by the contribution *ζij* from selection of alternative *j* ∈ *Ji* for section *i* ∈ *I*. Other continuous state variables *x* may also be relevant for internal calculations within the alternatives. In the methanol case study, these variables include conversion rates in the reactors and shaft work required in the compressors. Finally, the Boolean variables *Y* govern selection among the process alternatives and modeling detail levels. *Yij* determines whether alternative *j* ∈ *Ji* is active for section *<sup>i</sup>* ∈ *<sup>I</sup>*. *<sup>Y</sup><sup>l</sup> <sup>i</sup>* determines whether process section *i* ∈ *I* is modeled at detail level *l* ∈ *Li*. Finally, for each process alternative *<sup>j</sup>* ∈ *Ji* in section *<sup>i</sup>* ∈ *<sup>I</sup>*, Boolean *<sup>Y</sup><sup>l</sup> ij* determines the modeling detail level *l* ∈ *Lij*.


#### *3.3. Functions*

The problem-specific variable relationships for the end-to-end process are represented by several functions. The globally relevant constraints *g*(*f* , *F*, *T*, *P*, *x*) describe variable relationships that must be satisfied regardless of discrete selections of the process alternatives or modeling detail levels. These include the linking constraints that equate stream flow properties between different process sections. That is, the exit stream from the procurement section should be equivalent to the inlet stream to the production section. The constraints *rij*(*f* , *F*, *T*, *P*, *x*) describe the relationships that are enforced regardless of the selected detail level when alternative *j* ∈ *Ji* is selected for section *i* ∈ *I*. For each section *<sup>i</sup>* ∈ *<sup>I</sup>*, the constraints *<sup>h</sup><sup>l</sup> i* (*f* , *F*, *T*, *P*, *x*) describe variable interactions at each detail level *l* ∈ *Li* that are relevant regardless of the selected process alternatives. These constraints include potential equality relationships that link different process alternatives in a section with each other. For each of these alternatives *<sup>j</sup>* ∈ *Ji*, the constraints *<sup>s</sup><sup>l</sup> ij*(*f* , *F*, *T*, *P*, *x*) describe the variable relationships at detail level *l* ∈ *Lij*. Here are included the kinetic calculations for the reactor conversion, or the shaft work calculation for the compressors. The cost functions are computed using *φ<sup>l</sup> i* (*f* , *F*, *T*, *P*, *x*, *ζ*) at the section level, and *ψ<sup>l</sup> ij*(*<sup>f</sup>* , *<sup>F</sup>*, *<sup>T</sup>*, *<sup>P</sup>*, *<sup>x</sup>*) for each process alternative. One common interpretation of *<sup>φ</sup><sup>l</sup> i* (*f* , *F*, *T*, *P*, *x*, *ζ*) is given in Equation (3), where the section cost is simply equal to the sum of the contributions *ζij* from each alternative *j* ∈ *Ji*, but more complex relationships are possible.


$$z\_i = \phi\_i^l(f, F, T, P, \mathbf{x}, \zeta) = \sum\_{j \in I\_i} \zeta\_{ij} \tag{3}$$

#### *3.4. Formulation*

*Yl*

The overall generic problem formulation is given in Problem (P1). The objective is to maximize the profit, denoted by *Z*, equal to the summation of profit (or negative cost) contributions from each section *i* ∈ *I*. In the methanol process, we consider revenue from the methanol sales, purchase costs from the syngas feed, utility costs for the heaters and coolers, electricity costs for the compressors, fuel credit for the purge stream, and annualized capital costs for equipment purchases.

$$\max Z = \sum\_{i \in I} z\_i \tag{P1.1}$$

$$\begin{array}{cc} \text{s.t.} & \text{g}(f\_{\prime}F\_{\prime}T\_{\prime}P\_{\prime}\mathbf{x}) \leq 0\\ & \text{ $\! $ } & \text{ $\! $ } \end{array} \tag{\text{P1.2}}$$

$$\bigvee\_{j \in I\_i} \begin{bmatrix} \chi\_{ij} \\ r\_{ij}(f, F, T, P, x) \le 0 \end{bmatrix} \tag{P1.3}$$

$$\bigvee\_{l \in \mathbb{L}\_i} \Bigg| \Bigg( h\_i^l(f, F, T, P, \mathbf{x}) \le 0 \\\\ z\_i = \phi\_i^l(f, F, T, P, \mathbf{x}, \zeta) \Big)$$

$$\begin{array}{c} \Upsilon\_{ij} \implies \bigvee\_{l \in \mathcal{L}\_{ij}} \left[ \begin{array}{c} \Upsilon\_{ij}^{l} \\ s\_{ij}^{l}(f, F, T, P, \mathbf{x}) \le 0 \\ \zeta\_{ij} = \Psi\_{ij}^{l}(f, F, T, P, \mathbf{x}) \end{array} \right] \end{array} \tag{P1.5} \qquad \forall i \in I, \forall j \in J\_{i} \tag{P1.5}$$

$$\bigvee\_{l \in L\_i} Y\_i^l \tag{P1.6}$$

$$Y\_{ij} \implies \bigvee\_{l \in L\_{ij}} Y\_{ij}^l \tag{P1.7} \tag{P1.7}$$

$$\stackrel{d}{\; \; ij} \implies \text{Y}\_{\text{ij}} \tag{1.8} \tag{1.9} \qquad \qquad \qquad \qquad \forall i \in I, \forall j \in J\_i, \forall l \in L\_{\text{ij}} \tag{P1.8}$$

$$\begin{aligned} \Omega(\mathbf{Y}) &= True \\ z\_i \in \mathbb{R} & \forall i \in I \\ z\_{ij} \in \mathbb{R} & \forall i \in I, \forall j \in J\_i \end{aligned} \tag{P1.10}$$

$$\begin{cases} \mathcal{T} & \text{if } \mathcal{T} \in \mathbb{R}^{|\mathcal{L}|} \\ \mathcal{T} \in \mathbb{R}^{|\mathcal{L}| \left| K \right|} & \text{(P1.12)} \\ F \in \mathbb{R}^{|K|} & \text{(P1.13)} \\ T \in \mathbb{R}^{|K|} & \text{(P1.14)} \\ P \in \mathbb{R}^{|K|} & \text{(P1.15)} \\ \mathcal{X} \in X \subseteq \mathbb{R}^{n\_{X}} & \text{(P1.16)} \\ \text{(P1.16)} & \text{(P1.16)} \\ \mathcal{T} & \text{(P1.16)} \\ \mathcal{T} & \text{(P1.16)} \end{cases}$$

$$Y \in \{True, False\}^p \tag{P1.17}$$

(P1.2) describes global constraints that are enforced independent of any alternative or detail selection. Disjunction (P1.3) governs the selection among the process alternatives *j* ∈ *Ji* for each section *i* ∈ *I*. Disjunction (P1.4) gives the detail level *l* ∈ *Li* at which the major process sections *i* ∈ *I* are modeled. Implication (P1.5) states that the selection of an alternative *j* ∈ *Ji* implies the choice of a detail level *l* ∈ *Lij*. For each section, the exclusive-OR relationship (P1.6) states that exactly one detail level *l* ∈ *Li* is used to model alternative-independent interactions. Similarly, for a selection of process alternative *j* ∈ *Ji*, Implication (P1.7) governs selection of exactly one detail level *l* ∈ *Lij* for modeling each alternative. Implication (P1.8) enforces that selection of a modeling detail level for an alternative implies that the alternative is selected. Other logical propositions are expressed using Ω(*Y*). Finally, the continuous variable definitions are given in lines (P1.10)–(P1.17).

Note that this illustrative example is meant to give a sense of the complexity that is possible to represent in a process design problem using GDP modeling techniques. GDP modeling easily supports augmentation of the model to consider, for example, methane reforming as another process section. Other process relationships and logical expressions are also possible to include, as the problem demands.

#### *3.5. Discussion*

The GDP model in Section 3.4 captures both the choice among discrete process alternatives as well as the level of modeling detail used to describe each alternative. The observant reader may note that the MINLP resulting from the reformulation of this GDP is more complex than simply modeling the entire process in high detail. While true, the power of GDP lies in the ability to systematically activate or deactivate entire blocks of related constraints. The intention of formulation (P1) is not to solve the monolithic GDP model, but rather to systematically generate models that trade off fidelity and tractability for different analyses from a single source of truth by imposing the relevant logical implications. As a result, Problem (P1) could be regarded as a parametric optimization in which the modeler prespecifies the values of *Y<sup>l</sup> <sup>i</sup>* to satisfy Equation (P1.6) for each section *i* ∈ *I* and provides logical implications to tie selection of an alternative *Yij* to selection of the desired detail level *Y<sup>l</sup> ij* for each alternative *j* ∈ *Ji*.

Once these decisions are made, model simplifications driven by logical inference are applied to generate a process model at the desired level of detail for each of its constituent sections. For example, the production team may want to adopt a simplified view of the procurement and sales sections while preserving a high-fidelity view of production section alternatives. To accomplish this, the logical statements in Equation (4) may be appended to the model.

$$\begin{aligned} \text{Y}\_{i} & \implies \text{Y}\_{i}^{low}, \qquad i \in \{\text{proc}urement, sales\} \\ \text{Y}\_{i} & \implies \text{Y}\_{i}^{high}, \qquad i \in \{\text{production}\} \\ \text{Y}\_{ij} & \implies \text{Y}\_{ij}^{low}, \qquad i \in \{\text{proc}urement, sales\}, \forall j \in \mathbf{J}\_{i} \\ \text{Y}\_{ij} & \implies \text{Y}\_{ij}^{high}, \qquad i \in \{\text{production}\}, \forall j \in \mathbf{J}\_{i} \end{aligned} \tag{4}$$

Due to the exclusive-OR-type relationship between different levels of modeling detail established by logical statements (P1.6) and (P1.7), this forces the implied level of detail to be selected for its corresponding alternative or section. Applying standard principles of logical inference [21], we arrive at the model shown in Problem (P2), with the sets *<sup>I</sup>low* = {*procurement*,*sales*} and *<sup>I</sup>high* = {*production*}.

$$\max Z = \sum\_{i \in I} z\_i \tag{P2.1}$$

$$\text{s.t.}\quad g(f, \mathbf{F}, T, P, \mathbf{x}) \le 0 \tag{\text{P.2}}$$

$$\bigvee\_{j \in I\_i} \begin{bmatrix} \chi\_{ij} \\ r\_{ij}(f, F, T, P, \mathbf{x}) \le 0 \end{bmatrix} \tag{P2.3}$$

$$\begin{aligned} Y\_i^{low} &= True\\ h\_i^{low}(f, F, T, P, \mathbf{x}) &\le 0\\ z\_i &= \phi\_i^{low}(f, F, T, P, \mathbf{x}, \zeta) \end{aligned} \qquad \forall i \in I^{low} \tag{\mathbb{P}2.4}$$

*Yl*

$$Y\_i^l = False \tag{2.5}$$

$$\forall i \in I^{low}, \forall l \in L\_i \; \{low\} \tag{P2.5}$$

$$\begin{aligned} Y\_i^{\text{high}} &= True\\ h\_i^{\text{high}}(f, F, T, P, \mathbf{x}) &\le 0\\ z\_i &= \phi\_i^{\text{high}}(f, F, T, P, \mathbf{x}, \zeta) \end{aligned} \qquad \forall i \in I^{\text{high}} \tag{12.6}$$

$$Y\_{ij}^{l} \implies \begin{cases} Y\_{ij}^{l} \implies & \forall i \in I^{l \text{igh}}, \forall l \in L\_{i} \\ & Y\_{ij}^{law} \\ \int\_{ij} s\_{ij}^{low}(f, F, T, P, x) \le 0 \\ \int\_{ij} \psi\_{ij}^{low}(f, F, T, P, x) \end{cases} \quad \forall i \in I^{low}, \forall j \in J\_{i} \tag{P2.8}$$

$$I\_{ij}^l = \text{False} \tag{2.2}$$

$$\forall i \in I^{low}, \forall j \in J\_{i\cdot}, \forall l \in L\_i \; (\; \{low\} \tag{P2.9})$$

$$\Psi\_{ij} \implies \left\{ \begin{aligned} \mathbf{Y}\_{ij}^{\text{high}} \\ s\_{ij}^{\text{high}}(f, F, T, P, x) \le 0 \\ \zeta\_{ij} = \psi\_{ij}^{\text{high}}(f, F, T, P, x) \end{aligned} \right\}\_{i \in I} \quad \forall i \in I^{\text{high}}, \forall j \in J\_i \tag{12.10}$$

$$\begin{aligned} \mathcal{Y}\_{ij}^{l} &= \text{False} & \quad \forall i \in I^{\text{tag}}, \forall j \in J\_{i\cdot} \forall l \in L\_{i} \; \{\} \\ \mathcal{Y}\_{ij}^{l} & \implies \mathcal{Y}\_{ij} & \quad \forall i \in I, \forall j \in J\_{i\cdot} \forall l \in L\_{i\cdot} \end{aligned} \tag{P2.11}$$

$$\begin{aligned} \Omega(Y) &= True \\ z\_i \in \mathbb{R} \end{aligned} \tag{\mathbb{P}2.13}$$

$$z\_{i\uparrow} \in \mathbb{R} \tag{2.15} \tag{2.16} \tag{2.16} \tag{2.17} \tag{2.17}$$

$$\begin{aligned} f &\in \mathbb{R}^{|C||K|} \\ F &\in \mathbb{R}^{|K|} \\ T &\in \mathbb{R}^{|K|} \\ P &\in \mathbb{R}^{|K|} \\ &\forall x \in X \subseteq \mathbb{R}^{n\_{\mathcal{X}}} \end{aligned} \tag{\mathbb{P}2.17}$$

$$\mathcal{Y} \in \{True, False \}^p \tag{\text{P2.21}}$$

Problem (P2) now describes the decision space of the overall process with a focus on the production section. Changing the logical implications in Equation (4), we can easily shift focus in model fidelity to other sections of the process. By imposing different logical relationships on the general GDP model and applying easily automated principles of logical inference, we are able to derive multiple model variants from a single source of truth. Different levels of detail can also be evaluated as a post-solve solution quality check. The modeler can hold constant the production section decisions and increase modeling detail in the other sections to examine impacts of their decision-making on other sections. While this type of analysis may be possible under other engineered frameworks, GDP offers the formalism of an end-to-end perspective that is tied together by mathematical theory.

#### *3.6. Solution Strategies*

After generating a variant such as Problem (P2), a solution approach may be selected to obtain the optimal decision values. As previously introduced, the BM and HR reformulations to MINLP are the most popular approaches [2], trading off formulation size versus tightness of the continuous relaxation. The BM formulation for Problem (P2) may be found in Appendix A. For BM, equality constraints in the disjuncts are replaced by their corresponding two inequalities to facilitate relaxation of the constraint. BM results in a smaller problem size, as it does not require the introduction of new variables. However, the use of the Big-M parameter *M* results in a looser continuous relaxation. Note that the relaxation may be improved by selecting unique values of *M* for each constraint [48]. For a tighter continuous relaxation, the HR formulation found in Appendix B may be used. HR requires additional disaggregated variables to be defined for each disjunct, so the problem size may be significantly larger than BM. However, the tighter HR formulation may require fewer iterations to converge. Once the reformulation to MINLP is made, the model can be sent to the user's solver of choice. Direct logic-based decomposition approaches [22,39] are also possible for solving GDP models, with implementations available in Pyomo.GDP via the GDPopt solver [46].

We solve the methanol synthesis model described in Problem (P2) and obtain a solution with annual profit of \$1.8 million using two-stage feed compression, the cheap reactor, and single-stage recycle compression, see Figure 4. We can now fix the production section discrete decisions and evaluate the solution at varying detail levels for the procurement and sales sections. The results are given in Table 1.

**Figure 4.** Solution flowsheet for Problem (P2), using two-stage feed compression, the cheap reactor, and single-stage recycle compression. At low procurement modeling detail, no feed selection decisions are made.

**Table 1.** Profit (1000 USD) at the fixed production section design of two-stage feed compression, the cheap reactor, and single-stage recycle compression compared to the profit achievable when the design is allowed to vary.


Notice that the solution profit tends to increase at higher modeling details for the feed and sales. This is due to the additional flexibility in adjusting coordinating purity levels across procurement, production, and sales. In the high-detail solutions, a higher purity syngas is purchased to enable production of a higher purity methanol product. The supplier and customer contracts are also selected to facilitate the purity decision. The procurement contract is selected with a higher base cost, but lower incremental cost for improved purity. Conversely, the customer contract is chosen where a higher purity is more valued. Thus, despite an increase in feed cost of \$4.3 million vs. \$3.4 million in the low-detail solution, the product revenue rises to \$10.6 million rather than \$7.7 million.

The Problem (P2) solution is also compared in Table 1 to the profit possible when the production section configuration is allowed to change at higher procurement and sales modeling detail levels. Here, we see that at higher levels of modeling detail, the single stage feed compressor becomes more advantageous, but only by the \$50 thousand margin that accounts for the difference in annualized

capital cost. In general, this analysis may not be possible, as the formulation for solving all decision degrees of freedom at a high level of modeling detail may be intractable. However, even at a medium level of sales modeling detail, it is possible to notice that the choice of single-stage feed compression may be relevant in the optimal production configuration.

As an illustration of the flexibility of GDP modeling, another solution approach that can be utilized is to emulate traditional process design strategies. GDP model (P1) is compatible with a design analysis akin to the hierarchical decomposition approach described by Douglas [50], by solving sequentially at different detail levels. First, the overall process (P1) is solved with a low detail level in iteration *iter* = 0, enforcing the implications in Equation (5).

$$\begin{aligned} Y\_i^0 &\implies Y\_i^{0,low}, & \forall i \in I\\ Y\_{ij}^0 &\implies Y\_{ij}^{0,low}, & \forall i \in I\_\prime \forall j \in J\_i \end{aligned} \tag{5}$$

The solution to this overview problem gives high-level decisions among the process alternatives and defines the following sets for the iteration *iter* = 1. Let *J*<sup>1</sup> *<sup>i</sup>* = {*<sup>j</sup>* ∈ *Ji* : *<sup>Y</sup>*<sup>0</sup> *ij* = *True*} denote the selected alternatives from the overview problem from iteration *iter* = 0. In the next iteration, we enforce for (P1) the implications in Equation (6) such that these alternatives are evaluated at a progressively higher detail level.

$$\begin{aligned} \mathcal{Y}\_i^1 &\implies \mathcal{Y}\_i^{1,mod}, \qquad \forall i \in I\\ \mathcal{Y}\_{ij}^1 &\implies \mathcal{Y}\_{ij}^{1,mod}, \qquad \forall i \in I, \forall j \in J\_i^1\\ \mathcal{Y}\_{ij}^A &\implies \mathcal{Y}\_{ij}^{0,low}, \qquad \forall i \in I, \forall j \in J\_i \; \mid J\_i^1 \end{aligned} \tag{6}$$

The algorithm would terminate at iteration *iter* = *N* when all selected alternatives are evaluated at a high detail level: *Yhigh ij* = *True*, ∀*i* ∈ *I*, ∀*j* ∈ *Ji*. Note that solutions at lower detail levels can be used to initialize the higher fidelity models for each alternative, aiding in algorithm robustness. The amount of backtracking done by the algorithm to evaluate other alternatives can be tuned by applying a penalty factor to alternatives modeled at lower detail levels.

We apply this algorithm to the methanol synthesis. In the base iteration (*iter* = 0), all alternatives are modeled at a low detail level. The solution gives a configuration of the cheap reactor with two-stage feed compression and single-stage recycle compression at a profit of \$1.1 million. The structure for this solution is identical to that shown in Figure 4. We increase the modeling detail level for these selections and solve iteration 1.

Solving the low-detail overview problem, we obtain a configuration using the cheap reactor, two-stage feed compression, and single-stage recycle compression, yielding a profit of \$1.1 million. We then increase the level of modeling detail for the selected alternatives. Solving again, we obtain a profit of \$0.4 million with the expensive feed alternative, the cheap reactor, two-stage feed compression, single-stage recycle compression, and the high-purity sales option. At this intermediate level of modeling detail, the solution profit decreases because physical constraints are more tightly enforced, but not as many optimization degrees of freedom are made available to the solver yet. For some analyses, it may be advantageous for intermediate detail levels to produce a monotonically tightening approximation of the high-detail representation, but we do not address that consideration in this work. At the next iteration (*iter* = 2), we increase the detail level again, obtaining now a profit of \$3.1 million, using the contract with supplier 1, the cheap reactor, two-stage feed compression, single-stage recycle compression, and the contract with customer 1. At this point, the algorithm terminates, as all selected alternatives are modeled at the highest possible level of detail.

#### **4. Conclusions**

In this paper, we present GDP not only as a mathematical modeling framework supporting advanced solution strategies, but also as a modeling tool to organize model variants. We demonstrate through an illustrative case study that GDP is a useful model abstraction that separates algebraic and logical relationships within a process design problem. From a single GDP model capturing discrete design alternatives as well as alternatives in modeling fidelity, variants can be generated to suit the modeling scope and focus for different process sections by specifying appropriate logical implications. As a result, model interoperability is improved, facilitating simplified post-solution validation analysis. Preservation of model logical structure in GDP also offers the modeler flexibility to reformulate the problem as an MINLP or to apply a variety of automated decomposition methods.

**Author Contributions:** Q.C. wrote the manuscript. I.G. supervised the project.

**Funding:** We graciously acknowledge funding from the U.S. Department of Energy, Office of Fossil Energy's Crosscutting Research Program through the Institute for the Design of Advanced Energy Systems (IDAES).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Big-M Reformulation (BM)**

The Big-M reformulation for GDP model in Problem (P2) is presented in Equations (A1)–(A31). All variables' domains for BM are assumed to be a subset of the non-negative real numbers.



$$y \in \{0, 1\}^p \tag{A31}$$

#### **Appendix B. Hull Reformulation (HR)**

*ylow*

*hlow*

*yl*

*h high*

*Yhigh*

*zi* = *φlow*

The Hull Reformulation model for Problem (P2) is presented in Equations (A32)–(A70). Note that in HR, the conditional nonlinear functions, Equations (A34), (A55), (A56), (A59), and (A60) use a perspective function formulation. To avoid the singularity at a zero value of the indicator binary variable, various reformulations have been proposed in literature [39,43,51] and are available in Pyomo.GDP.

$$\max Z = \sum\_{i \in I} z\_i \tag{A32}$$

$$\text{s.t.} \quad g(f, \mathbf{F}, T, P, \mathbf{x}) \le 0 \tag{A33}$$

$$y\_{ij}r\_{ij}\left(\frac{\hat{f}\_{ij}}{\mathcal{Y}\_{ij}}, \frac{\hat{\mathcal{F}}\_{ij}}{\mathcal{Y}\_{ij}}, \frac{\hat{\mathcal{T}}\_{ij}}{\mathcal{Y}\_{ij}}, \frac{\hat{\mathcal{P}}\_{ij}}{\mathcal{Y}\_{ij}}, \frac{\hat{\mathcal{R}}\_{ij}}{\mathcal{Y}\_{ij}}\right) \le 0 \tag{A34} \qquad \forall i \in I, j \in J\_i \tag{A34}$$

$$\sum\_{j \in I\_i} y\_{ij} \ge 1 \tag{A35}$$

$$f^{lb} y\_{ij} \le \hat{f}\_{ij} \le f^{ub} y\_{ij} \tag{A36} \\ \qquad \qquad \qquad \forall i \in I\_i \forall j \in J\_i \tag{A36}$$

$$\begin{array}{ccccc} F^{lb} y\_{ij} \le \mathring{F}\_{ij} \le F^{ub} y\_{ij} & & \forall i \in I, \forall j \in J\_i \\ \ldots & & \downarrow & \downarrow & \downarrow \\ & & \downarrow & \downarrow & \downarrow \end{array} \tag{A37}$$

$$\begin{aligned} T^{lb} y\_{ij} \le \hat{\Gamma}\_{ij} \le T^{ub} y\_{ij} & \forall i \in I, \forall j \in J\_i\\ P^{lb} y\_{ij} \le \hat{P}\_{ij} \le P^{ub} y\_{ij} & \forall i \in I, \forall j \in J\_i\\ x^{lb} y\_{ij} \le \pounds\_{ij} \le x^{ub} y\_{ij} & \forall i \in I, \forall j \in J\_i \end{aligned} \tag{A38}$$

$$f\_{\vec{i}\vec{j}} = \sum\_{\vec{j}\in I\_{\vec{i}}} \hat{f}\_{\vec{i}\vec{j}} \tag{A41} = \sum\_{\vec{j}\in I\_{\vec{i}}} \hat{f}\_{\vec{i}\vec{j}} \tag{A42}$$

$$F\_{\vec{ij}} = \sum\_{\vec{j} \in f\_i} \vec{\Gamma}\_{\vec{ij}} \tag{A42}$$
 
$$T\_{\vec{ij}} = \sum \Upsilon\_{\vec{ij}} \tag{A43} \tag{A44} \tag{A44}$$

$$P\_{ij} = \sum\_{j \in f\_i} \hat{P}\_{ij} \tag{A44}$$
 
$$\sum\_{i \in I} \tag{A45} \tag{A46}$$

$$\mathbf{x}\_{ij} = \sum\_{j \in f\_i} \mathbf{\hat{x}\_{ij}} \tag{A45} \tag{A46}$$
 
$$y\_i^{low} = 1 \tag{A46} \tag{A46} \\ \mathbf{\hat{x}} \in I^{low} \tag{A46}$$
  $\text{where } \mathbf{\hat{x}} \text{ represents} \tag{A46}$ 

*<sup>i</sup>* (*f* , *F*, *T*, *P*, *x*) ≤ 0 ∀*i* ∈ *I low* (A47) *<sup>i</sup>* (*f* , *F*, *T*, *P*, *x*, *ζ*) ∀*i* ∈ *I low* (A48)

$$\begin{aligned} \frac{d}{dt} &= 0 & \forall i \in I^{low}, \forall l \in L\_i \; (\;\{low\} \\ \frac{d}{dt} &= 1 & \forall i \in I^{high} \\ \frac{d}{dt} (f, F, T, P, x) &\le 0 & \forall i \in I^{high} \end{aligned} \tag{A50}$$

*yl*

*y*

*ylow*

$$z\_i = \phi\_i^{high}(f, F, T, P, \mathbf{x}, \mathbb{Q}) \tag{A52} \\ \tag{A52}$$

$$\frac{1}{i} = 0 \tag{4.5}$$

$$\forall i \in I^{high}, \forall l \in L\_i \; (\; \{hights\} \tag{A53}$$

$$\mathbf{j}\_{ij}^{low} = y\_{ij} \tag{4.54}$$
 
$$\forall i \in I^{low}, \forall j \in J\_i \tag{A54}$$

$$(y\_{ij}, \stackrel{low}{\underset{ij}{\colon}}\_{\stackrel{low}{\underset{ij}{\colon}}} \left( \frac{f\_{ij}}{y\_{ij}}, \frac{f\_{ij}}{y\_{ij}}, \frac{f\_{ij}}{y\_{ij}}, \frac{g\_{ij}}{y\_{ij}} \right) \le 0 \tag{A55} \\ \qquad \qquad \forall i \in I^{low}, \forall j \in J\_i \tag{A55}$$

$$\begin{aligned} \mathcal{Z}\_{ij} &= y\_{ij} \psi\_{ij}^{low} \left( \frac{\hat{f}\_{ij}}{y\_{ij}}, \frac{\hat{P}\_{ij}}{y\_{ij}}, \frac{\hat{T}\_{ij}}{y\_{ij}}, \frac{\hat{P}\_{ij}}{y\_{ij}}, \frac{\p\_{ij}}{y\_{ij}} \right) & \forall i \in I^{low}, \forall j \in J\_i \\\ y\_{ij}^{l} &= 0 & \forall i \in I^{low}, \forall j \in J\_i, \forall l \in L\_i \end{aligned} \tag{A56}$$

$$\begin{aligned} \stackrel{I}{\underset{ij}{=}} &= 0 & \quad \forall i \in I^{low}, \forall j \in J\_i, \forall l \in L\_i \; \backslash \; \{low\} \end{aligned} \tag{A57}$$
 
$$\stackrel{\text{high}}{=} \quad \forall i \in I^{high}, \forall j \in J\_i \tag{A58}$$

$$y\_{ij} s\_{ij}^{\text{high}} \left( \frac{\hat{\mathbf{f}}\_{ij}}{y\_{ij}}, \frac{\hat{\mathbf{f}}\_{ij}}{y\_{ij}}, \frac{\hat{\mathbf{f}}\_{ij}}{y\_{ij}}, \frac{\hat{\mathbf{f}}\_{ij}}{y\_{ij}}, \frac{\hat{\mathbf{x}}\_{ij}}{y\_{ij}} \right) \le 0 \tag{A59} \qquad \forall i \in I^{\text{high}}, \forall j \in J\_i \tag{A59}$$
 
$$\zeta\_{ij} = y\_{ij} \psi\_{ij}^{\text{high}} \left( \frac{f\_{ij}}{y\_{ij}}, \frac{\hat{\mathbf{f}}\_{ij}}{y\_{ij}}, \frac{\hat{\mathbf{T}}\_{ij}}{y\_{ij}}, \frac{\hat{\mathbf{P}}\_{ij}}{y\_{ij}}, \frac{\hat{\mathbf{x}}\_{ij}}{y\_{ij}} \right) \tag{A60} \qquad \forall i \in I^{\text{high}}, \forall j \in J\_i \tag{A60}$$

$$\begin{aligned} \mathcal{I}\_{ij} = y\_{ij} \psi\_{ij}^{high} \left( \frac{f\_{ij}}{y\_{ij}}, \frac{f\_{ij}}{y\_{ij}}, \frac{\Upsilon\_{ij}}{y\_{ij}}, \frac{\Psi\_{ij}}{y\_{ij}}, \frac{\mathfrak{s}\_{ij}}{y\_{ij}} \right) & \forall i \in I^{high}, \forall j \in J\_i \\\ y\_{ij}^{l} = 0 & \forall i \in I^{high}, \forall j \in J\_i, \forall l \in L\_i \end{aligned} \tag{A60}$$

$$\begin{aligned} \dot{\Omega}y &= b \\ z\_i \in \mathbb{R} \end{aligned} \tag{A62}$$

$$z\_{i\uparrow} \in \mathbb{R} \tag{4.64} \tag{4.65} \tag{4.66} \tag{4.66}$$

$$f \in \mathbb{R}^{|\mathcal{C}||K|} \tag{A65}$$
 
$$r \in \mathbb{R}^{|K|} \tag{A66}$$

$$F \in \mathbb{R}^{|K|} \tag{A66}$$

$$T \in \mathbb{R}^{|K|} \tag{A67}$$

$$T = \dots = |V| \tag{A67}$$

$$P \in \mathbb{R}^{|K|} \tag{A68}$$

$$\chi \in X \subseteq \mathbb{R}^{n\_{\chi}} \tag{A69}$$

$$y \in \{0, 1\}^p \tag{A70}$$

#### **Appendix C. Methanol Model**

The methanol synthesis model examined in this paper is adapted from the description in Example 3 of [22]. In this appendix, we describe the main differences between the literature model and the presented variant. We consider the original unit models to be the high-detail version, except for the feed and product description. The original feed model is considered to be the medium-detail description, and the original product sales model is considered the low-detail description. Linking constraints for stream flows and material balance equations are preserved from the original formulation.

#### *Appendix C.1. Feed Procurement*

#### Appendix C.1.1. High Detail

The high-detail feed model involves the choice between two different suppliers, with purity-dependent costs for the syngas feed.

Supplier 1

$$f\_{1,H2} = 0.7(F\_1 - f\_{1,CH\_4})\tag{A71}$$

$$f\_{1,CO} = 0.3(F\_1 - f\_{1,CH\_4})\tag{A72}$$

$$0.03F\_1 \le f\_{1,CH\_4} \le 0.20F\_1\tag{A73}$$

$$z\_1 = -425.712 + 194.973(\ln f\_{1,CH\_4} - \ln F\_1)F\_1 \tag{A74}$$

Supplier 2

$$f\_{2,H\_2} = 0.7(F\_2 - f\_{2,CH\_4})\tag{A75}$$

$$f\_{2,CO} = 0.3(F\_2 - f\_{2,CH\_4})\tag{A76}$$

$$0.03F\_2 \le f\_{2,CH\_4} \le 0.20F\_2\tag{A77}$$

$$z\_1 = -400 + 210(\ln f\_{2,CH\_4} - \ln F\_2)F\_2\tag{A78}$$

#### Appendix C.1.2. Medium Detail

The medium-detail feed model is the same as in [22], with the choice between a cheaper and a more expensive feed stream.

Cheap Feed

$$f\_{1,H\_2} = 0.6F\_1\tag{A79}$$

$$f\_{1,CO} = 0.25F\_1\tag{A80}$$

$$f\_{1,CH\_4} = 0.15F\_1\tag{A81}$$

$$z\_1 = -795.6F\_1\tag{A82}$$

$$f\_{1,H\_2} = 0.65F\_1\tag{A83}$$

$$f\_{1,CO} = 0.30F\_1\tag{A84}$$

$$f\_{1,CH\_4} = 0.05F\_1\tag{A85}$$

$$z\_1 = -1009.8F\_1\tag{A86}$$

#### Appendix C.1.3. Low Detail

The low-detail feed model simply defaults to the expensive feed option from the medium-detail case, with no other option available.

$$f\_{1,H\_2} = 0.65F\_1\tag{A87}$$

$$f\_{1,CO} = 0.30F\_1\tag{A88}$$

$$f\_{1,CH\_4} = 0.05F\_1\tag{A89}$$

$$z\_1 = -1009.8F\_1\tag{A90}$$

#### *Appendix C.2. Product Sales*

#### Appendix C.2.1. High Detail

Customer 1

$$0.85F\_{23} \le f\_{23,CH\_3OH} \le 0.98F\_{23} \tag{A91}$$

$$z\_3 = 1682.92 - 2275.6(\ln(f\_{23,H\_2} + f\_{23,CO} + f\_{23,CH\_4}) - \ln F\_{23})F\_{23} \tag{A92}$$

Customer 2

$$0.85F\_{23} \le f\_{23,CH\_3OH} \le 0.98F\_{23} \tag{A93}$$

$$z\_3 = 1700 - 2265(\ln(f\_{23,H\_2} + f\_{23,CO} + f\_{23,CH\_4}) - \ln F\_{23})F\_{23} \tag{A94}$$

Appendix C.2.2. Medium Detail

Low purity product

$$10.85F\_{23} \le f\_{23,CH\_3OH} \tag{A95}$$

$$z\_3 = 6000F\_{23} \tag{A96}$$

High purity product

$$0.95F\_{23} \le f\_{23,CH\_3OH} \tag{A97}$$

$$z\_3 = 8500F\_{23} \tag{A98}$$

Appendix C.2.3. Low Detail

$$\begin{aligned} 0.9F\_{23} &\le f\_{23,CH\_3OH} \\ z\_3 &= 7650F\_{23} \end{aligned} \tag{A99}$$

*Appendix C.3. Compressors*

$$0 \le P\_{\text{ratio}} \le 1.74 \\ \tag{A101}$$

$$0 \le W \le 50 \\ \tag{A102}$$

Appendix C.3.1. High Detail

$$T\_{out} = P\_{ratio} T\_{in} \tag{A103}$$

$$\mathcal{W} = \frac{0.72(P\_{\text{ratio}} - 1)T\_{\text{in}}F\_{\text{in}}}{(10)(0.75)(0.23077)}\tag{A104}$$

$$P\_{out}^{0.23077} = P\_{ratio} P\_{in}^{0.23077} \tag{A105}$$

Appendix C.3.2. Medium Detail

$$P\_{ratio} = \begin{cases} 1.3 & \text{if feed compressor} \\ 1.1 & \text{if recycle compressor} \end{cases} \tag{A106}$$

$$T\_{out} = \begin{cases} 1.3T\_{in} & \text{if feed compressor} \\ 1.1T\_{in} & \text{if recycle compressor} \end{cases} \tag{A107}$$

$$\mathcal{W} \ge \begin{cases} \frac{0.72(0.3)T\_{\text{in}}F\_{\text{in}}}{(10)(0.75)(0.23077)} & \text{if feed compressor} \\ 0.72(0.1)T\_{\text{in}}F\_{\text{in}} & \text{if recycle compressor} \end{cases} \tag{A108}$$

(10)(0.75)(0.23077) if recycle compressor *P*0.23077 *out* = 1.3*P*0.23077 *in* if feed compressor 1.1*P*0.23077 *in* if recycle compressor (A109)

Appendix C.3.3. Low Detail

$$P\_{ratio} = \begin{cases} 1.3 & \text{if feed compressor} \\ 1.1 & \text{if rece cycle compressor} \end{cases} \tag{A110}$$

$$T\_{out} = \begin{cases} 1.3T\_{in} & \text{if feed compressor} \\ 1.1T\_{in} & \text{if recycle compressor} \end{cases} \tag{A111}$$

$$W \ge 0.54F\_{in} \tag{A112}$$

$$P\_{out} = 3.8P\_{in} \tag{A113}$$

*Appendix C.4. Expansion Valve*

Appendix C.4.1. High Detail

$$T\_{\rm out} P\_{\rm in}^{0.23077} = T\_{\rm in} P\_{\rm out}^{0.23077} \tag{A114}$$

$$P\_{out} \le P\_{in} \tag{A115}$$

Appendix C.4.2. Medium Detail

$$
\begin{bmatrix} T\_{out} = 0.94 T\_{in} \\ 0.75 P\_{in} \le P\_{out} \le 0.9 P\_{in} \end{bmatrix} \vee \begin{bmatrix} T\_{out} = 0.9 T\_{in} \\ 0.6 P\_{in} \le P\_{out} \le 0.75 P\_{in} \end{bmatrix} \tag{A116}
$$

Appendix C.4.3. Low Detail

$$T\_{out} = 0.94 T\_{in} \tag{A117}$$

$$P\_{\rm out} = 0.77 P\_{\rm in} \tag{A118}$$

*Appendix C.5. Cooler*

Appendix C.5.1. High Detail

$$Q = 0.00306(35.0)(F\_{in}T\_{in} - F\_{out}T\_{out})\tag{A119}$$

Appendix C.5.2. Medium Detail

$$
\begin{bmatrix} Q = 1.61(T\_{\text{in}} - T\_{\text{out}}) \\ 10 \le F\_{\text{in}} \le 20 \end{bmatrix} \vee \begin{bmatrix} Q = 0.80(T\_{\text{in}} - T\_{\text{out}}) \\ 5 \le F\_{\text{in}} \le 10 \end{bmatrix} \vee \begin{bmatrix} Q = 0.54(T\_{\text{in}} - T\_{\text{out}}) \\ 0 \le F\_{\text{in}} \le 5 \end{bmatrix} \tag{A120}
$$

Appendix C.5.3. Low Detail

$$Q = 1.0(T\_{in} - T\_{out})\tag{A121}$$

*Appendix C.6. Heater*

Appendix C.6.1. High Detail

$$Q = 0.00306(35.0)(F\_{out}T\_{out} - F\_{in}T\_{in})\tag{A122}$$

Appendix C.6.2. Medium Detail

$$
\begin{bmatrix} Q = 1.8(T\_{\text{out}} - T\_{\text{in}}) \\ 10 \le F\_{\text{in}} \le 20 \end{bmatrix} \vee \begin{bmatrix} Q = 0.80(T\_{\text{out}} - T\_{\text{in}}) \\ 5 \le F\_{\text{in}} \le 10 \end{bmatrix} \vee \begin{bmatrix} Q = 0.11(T\_{\text{out}} - T\_{\text{in}}) \\ 0 \le F\_{\text{in}} \le 5 \end{bmatrix} \tag{A123}
$$

Appendix C.6.3. Low Detail

$$Q = 1.0(T\_{out} - T\_{in}) \tag{A124}$$

*Appendix C.7. Reactors*

Appendix C.7.1. High Detail

$$P^2 P\_{sq,inv} = 1\tag{A125}$$

$$T\ T\_{inv} = 1\tag{A126}$$

$$r = \mathcal{X}f\_{\text{in},H\_2} \tag{A127}$$

$$(F\_{\text{in}}T\_{\text{in}} - F\_{\text{out}}T\_{\text{out}}) (35.0) = 0.01 \Delta H\_{\text{rxn}} r \tag{A128}$$

$$\mathcal{X}\_{\text{eq}} = 0.415(1 - (26.25e^{-18T\_{inv})P\_{\text{sig},inv}}) \tag{A129}$$

$$\mathcal{X}F\_{\text{in}} = \begin{cases} \mathcal{X}\_{\text{eq}}'(1 - \varepsilon^{-10})(f\_{\text{in},H\_2} + f\_{\text{in},CO} + f\_{\text{in},CH\_3OH}) & \text{expression vector} \\ \mathcal{X}\_{\text{eq}}'(1 - \varepsilon^{-5})(f\_{\text{in},H\_2} + f\_{\text{in},CO} + f\_{\text{in},CH\_3OH}) & \text{cheap reactor} \end{cases} \tag{A130}$$

Appendix C.7.2. Medium Detail

$$T = 5.23\tag{A131}$$

$$P^2 P\_{sq,inv} = 1\tag{A132}$$

$$
\sigma = \mathcal{X} f\_{\text{in}, \mathcal{H}\_2} \tag{A133}
$$

$$\mathcal{X}\_{t\eta} = 0.415(1 - (26.25e^{-18/5.23}P\_{sq,inv})) \tag{A134}$$

$$\mathcal{A}F\_{\rm in} = \begin{cases} \lambda\_{\rm eq}' (1 - e^{-10}) (f\_{\rm in, H\_2} + f\_{\rm in, CO} + f\_{\rm in, CH\_3OH}) & \text{expansion reactor} \\ \lambda\_{\rm eq}' (1 - e^{-5}) (f\_{\rm in, H\_2} + f\_{\rm in, CO} + f\_{\rm in, CH\_3OH}) & \text{cheap reactor} \end{cases} \tag{A135}$$

Appendix C.7.3. Low Detail

$$T \ge 5.23\tag{A136}$$

$$\mathcal{X} = \begin{cases} 0.35 & \text{expression reactor} \\ 0.30 & \text{cheap reactor} \end{cases} \tag{A137}$$

$$P \geq 14.3 \tag{A138}$$

$$r = \mathcal{X}f\_{\text{in}\_{\text{}}H\_2} \tag{A139}$$

*Appendix C.8. Flash*

Appendix C.8.1. High Detail

$$(13.6333 - \ln(7500.6168 P\_{H\_2}^{\text{nip}}))(100T - 3.19) = 164.9\tag{A140}$$

$$(14.3686 - \ln(7500.6168 P\_{\text{CO}}^{\text{nap}}))(100T + 13.15) = 530.22\tag{A141}$$

$$(18.5875 - \ln(7500.6168 P\_{CH\_3OH}^{np}))(100T + 34.29) = 3626.55\tag{A142}$$

$$(15.2243 - \ln(7500.6168 P\_{CH\_4}^{\text{nop}}))(100T + 7.16) = 897.84\tag{A143}$$

$$\mathfrak{g}\_{H\_2}(\mathfrak{g}\_{\text{CO}}P\_{H\_2}^{\text{rap}} + (1 - \mathfrak{g}\_{\text{CO}})P\_{\text{CO}}^{\text{rap}}) = P\_{H\_2}^{\text{rap}}\mathfrak{g}\_{\text{CO}} \tag{A144}$$

$$\mathfrak{J}\_{H\_2}(\mathfrak{J}\_{CH\_3OH}P\_{H\_2}^{np} + (1 - \mathfrak{J}\_{CH\_3OH})P\_{CH\_3OH}^{npp}) = P\_{H\_2}^{npp}\mathfrak{J}\_{CH\_3OH} \tag{A145}$$

$$\xi\_{H\_2}(\mathbb{J}\_{CH\_4}P\_{H\_2}^{up} + (1 - \mathbb{J}\_{CH\_4})P\_{CH\_4}^{up}) = P\_{H\_2}^{up}\xi\_{CH\_4} \tag{A146}$$

$$f\_{21,H\_2} = \pounds\_{H\_2} f\_{20,H\_2} \tag{A147}$$

$$f\_{21, \text{CO}} = \text{\textsuperscript{x}}\_{\text{CO}} f\_{20, \text{CO}} \tag{A148}$$

$$f\_{21,CH\_3OH} = \mathbb{Z}\_{CH\_3OH} f\_{20,CH\_3OH} \tag{A149}$$

$$f\_{21,CH\_4} = \J\_{CH\_4} f\_{20,CH\_4} \tag{A150}$$

$$\text{PF22} = P\_{H\_2}^{\text{rap}} f\_{22,H\_2} + P\_{\text{CO}}^{\text{np}} f\_{22,\text{CO}} + P\_{\text{CH}\_3\text{OH}}^{\text{np}} f\_{22,\text{CH}\_3\text{OH}} + P\_{\text{CH}\_4}^{\text{np}} f\_{22,\text{CH}\_4} \tag{A151}$$

#### Appendix C.8.2. Medium Detail

$$P\_{H\_2}^{\text{np}} = 73.332\tag{A152}$$

$$P\_{\rm CO}^{\rm mp} = 64.232\tag{A153}$$

$$P\_{CH\_3OH}^{\text{vup}} = 3.722\tag{A154}$$

$$P\_{CH\_4}^{\text{vup}} = 60.125\tag{A155}$$

$$
\xi\_{H\_2} (73.332 \xi\_{CO} + 64.232 (1 - \xi\_{CO})) = 73.332 \xi\_{CO} \tag{A156}
$$

$$\mathcal{J}\_{H\_2}(73.332\mathcal{J}\_{CH\_3OH} + 3.722(1 - \mathcal{J}\_{CH\_3OH})) = 73.332\mathcal{J}\_{CH\_3OH} \tag{A157}$$

$$
\pounds\_{H\_2}(73.332\pounds\_{CH\_4} + 60.125(1 - \pounds\_{CH\_4})) = 73.332\pounds\_{CH\_4} \tag{A158}
$$

$$f\_{21,H\_2} = \mathcal{J}\_{H\_2} f\_{20,H\_2} \tag{A159}$$

$$f\_{21, \text{CO}} = \pounds\_{\text{CO}} f\_{20, \text{CO}} \tag{A160}$$

$$f\_{21,CH\_3OH} = \pounds\_{CH\_3OH} f\_{20,CH\_3OH} \tag{A161}$$

$$f\_{21,CH\_4} = \J\_{CH\_4}^{\times} f\_{20,CH\_4} \tag{A162}$$

$$\text{PF}\_{22} = 73.332f\_{22,H\_2} + 64.232f\_{22,CO} + 3.722f\_{22,CH\_3OH} + 60.125f\_{22,CH\_4} \tag{A163}$$

Appendix C.8.3. Low Detail

$$
\zeta\_{H\_2} = 0.99 \tag{A164}
$$

$$
\zeta\_{\rm CO} = 0.99 \,\tag{A165}
$$

$$
\xi\_{CH\_3OH} = 0.85\tag{A166}
$$

$$\mathcal{J}\_{CH\_4} = 0.99 \tag{A167}$$

$$f\_{21,H\_2} = 0.99 f\_{20,H\_2} \tag{A168}$$

$$f\_{21,CO} = 0.99 f\_{20,CO} \tag{A169}$$

$$f\_{21,CH\_3OH} = 0.85 f\_{20,CH\_3OH} \tag{A170}$$

$$f\_{21,CH\_4} = 0.99 f\_{20,CH\_4} \tag{A171}$$

$$T = 4 \tag{A172}$$

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Review*
