Optimization of a refinery scheduling process with column generation and a quantum annealer

Vol.:(0123456789) Optimization and Engineering https://doi.org/10.1007/s11081-021-09662-8 1 3 RESEARCH ARTICLE Optimization ofarefinery scheduling process withcolumn generation andaquantum annealer J.Ossorio‑Castillo1· F.Pena‑Brage1,2 Received: 22 July 2020 / Revised: 21 April 2021 / Accepted: 3 July 2021 © The Author(s) 2021 Abstract This study focuses on the optimization of a refinery scheduling process with the help of an adiabatic quantum computer, and more concretely one of the quantum annealers developed by D-Wave Systems. We present an algorithm for finding a global optimal solution of a MILP that leans on a solver for QUBO problems, and apply it to various possible cases of refinery scheduling optimization. We analyze the inconveniences found during the whole process, whether due to the heuristic nature of D-Wave or the implications of reducing a MILP to QUBO, and present some experimental results. Keywords Refinery scheduling· Quantum annealing· Mixed-integer programming· Column generation 1 Introduction Since the first quantum algorithms, the range of problems where quantum computers can be applied have grown over time. Efforts have progressed in two fronts: to design algorithms that solve practical problems and to have operational quantum machines. In the first front, previous works on this matter that aimed to solve real-life problems with the help of a quantum annealer include Bauckhage etal. (2020), Calude and Dinneen (2017) and Venturelli etal. (2015). These articles consider a wide range of problems: the broadcast time problem, the job-shop scheduling problem and the Part of this research was developed as an activity in the Joint Research Unit Repsol-ITMATI (Code File: IN853A 2014/03) which is funded by FEDER, the Galician Agency for Innovation (GAIN) and the Ministry of Economy and Competitiveness in the framework of the Spanish Strategy for Innovation in Galicia. * F. Pena-Brage [email protected] 1 Instituto Tecnolóxico de Matemática Industrial (ITMATI), SantiagodeCompostela, Spain 2 Univ. de Santiago de Compostela, SantiagodeCompostela, Spain J.Ossorio-Castillo, F.Pena-Brage 1 3 max-sum diversification, respectively. However, all of them solely have binary variables in their formulations. In the second front, the possibility of exploiting the advantages of quantum computation over classical computers has begun to take form. The standard gate-based quantum computing is more natural for computer scientists and is the one used in most textbooks(Nannicini 2020). Adiabatic quantum computers are equivalent and they are well suited for optimization problems(Aharonov etal. 2008). Various companies, and more especially D-Wave Systems (D-Wave 2016), have been researching a way to physically implement these new computing paradigms prophesied by Feynman (1982) and Born and Fock (1928). Although the debate of whether the D-Wave machine really exploits quantum phenomena continues generating controversy, the fact that these supercomputers are especially predisposed to solve quadratic optimization problems with binary variables is indisputable (Calude and Dinneen 2017; Syrichas and Crispin 2017; Venturelli etal. 2015). Our aim is to take a problem of industrial interest and to follow the whole process of adapt it in order to solve it with a quantum annealer. We have choosen the scheduling of the arrival of vessels to the harbor in an oil refinery, the unloading of its contents into tanks, and also to get the most of them during the whole refining process. This classical problem in the oil and gas industry presents the advantage that its simpler version as mixed-integer linear problem (MILP) still retains the practical meaning(Lee etal. 1996). Surveys of representative works on optimization in the oil and gas industry can be found in Khor and Varvarezos (2017) and extensively in Furman etal. (2017). Mathematically, this kind of problems usually take a huge amount of constraints and variables (Karuppiah etal. 2008; Lee etal. 1996; Mouret 2010), a feature that clearly complicates the finding of the global solution within a reasonable amount of time. The problem of finding the global optimum of a MILP is an NP-Hard problem, and appears in multiple scenarios and applications. Our contribution and the main objective behind this work is to solve a specific MILP, the aforementioned scheduling problem, defined by binary and continuous variables with the help of a quantum annealer developed by D-Wave (see D-Wave (2016) for more details). To take advantage of the capacity of such computer, a decomposition technique is presented in this work. The approach we have developed consists of a combination of columngeneration algorithms such as the Dantzig–Wolfe decomposition (Dantzig and Wolfe 1960) and a branch-and-price method that correctly obtains an integer solution for the binary variables (Gamrath 2010). It is precisely the part that calculates the new columns the one that requires the solution of a binary linear problem, known as BLP, ZOLP or 0-1LP, where we can take advantage of the capacity of such quantum annealer. Although the refinery problem formulations, the column generation and the decomposition techniques present in this work are not novel, our contribution resides in using those techniques for separating the real part from the binary part of those problems and solving the latter one with a quantum annealer. Other hybrid methods that exploit the complementary strengths of quantum and classical computers can be consulted in Ajagekar etal. (2020) and Ajagekar (2020). The embedding of the binary problem into D-Wave requires transforming it into a Quadratic Unconstrained problem. However, it cannot be directly solved even in 1 3 Optimization ofarefinery scheduling process withcolumn… that case, since the available number of qubits is quite limited and the topology of the Chimera graph is rather specific. In this work we study the use the Qbsolv library(Booth etal. 2017; D-Wave Systems Inc. 2017) to overcome this difficulty and its extra cost in iterations. Due to the requirement of having a global solution in several steps of the algorithm, we have also checked that the solution has this property. Finally, we have tested the program with reference problems within the field of oil and gas. This paper is organized as follows: in Sect.2, we present the scheduling problem from an applied point of view. In Sect.3, we describe the column generation technique that best suits our goal, and the branch-and-price method that completes the algorithm to decompose our problem in the mixed and binary parts. In Sect.4, we describe more specifically how to solve the BLP part of the algorithm with a quantum annealer. Finally, in Sect.5 we present some experimental results: we solve a small refinery problem with an actual quantum annealer, the D-Wave 2X processor based in the University of Southern California. 2 Overview ofarefinery scheduling process First, we introduce the main mathematical problem behind the optimization of the refinery scheduling process, which can generally be modeled as a MILP. As the main objective of this work consists in optimizing the scheduling process of a refinery with the help of a quantum annealer, we have to separate the real-valued variables of the problem from the integer ones. For this purpose we shall use the Dantzig–Wolfe decomposition, which we explain in the next section. The refinery model has to take into account the variables and parameters involved in the refinery scheduling operations, such as the unloading of the vessels and the charging and storing into the tanks. It also has to define the sets that include all physical units: vessels, tanks, resources, etc. The model we have used as a basis for our algorithm can be found in Lee etal. (1996). This is a classical problem of inventory management of a refinery that imports several types of crude oil which are delivered by different vessels. The problem involves bilinear equations due to mixing operations. However, the linearity in the form of a MILP is maintained by replacing bilinear terms with individual component flows. This exact linear reformulation is possible since this scheduling system involves only mixing operation without splitting operation. More details about this problem and other possible reformulations can be seen in Mouret (2010) and Vyskocil and Djidjev (2019). The scheduling problem consists of a multistage system composed of NV vessels, NS storage tanks, NC charging tanks and ND distillation units, with NC key components of crude oil, as illustrated in Fig.1. For 1 ≤ 𝜈 ≤ NV , the 𝜈 th ship arrives at time TA,𝜈 with a volume of crude at initial time of VV,𝜈,0 . The cost of unloading a vessel per time unit is CU , and the cost of waiting in the sea per time unit is CS . For 1 ≤ i ≤ NS , the ith storage tank has a volume of crude at initial time of VS,i,0 and its load can vary from a minimum of VS,i,min to a maximum of VS,i,max , whereas the crude transferred from the 𝜈 th ship to the ith storage tank can vary J.Ossorio-Castillo, F.Pena-Brage 1 3 from FVS,𝜈,i,min to FVS,𝜈,i,max . The inventory cost for the storage tanks per unit of time and volume is CIST . For 1≤j≤NC , the mixed crude in the jth charging tank has an initial volume of VC,j,0 and it can vary from VC,j,min to VC,j,max , having a demand of Dj . The crude transferred from the ith storage tank can vary from FSC,i,j,min to FSC,i,j,max . The inventory cost for the charging tanks per unit of time and volume is CICT . Finally, for 1≤l≤ND , the crude transferred from the jth charging tank to the lth distillation unit can vary from FCD,j,l,min to FCD,j,l,max . In a distillation unit the changeover of crude from a charging tank to another has a cost CC≥0 . Our problem considers a scheduling horizon discretized in S time units. Following the previous index notations and for each time unit t, 1≤t≤S , some continuous variables must be determined (see Fig.2). For the 𝜈 th ship, we need to know the time point when the unloading starts tU,𝜈 and ends tD,𝜈 . We also have to find out at time t, • the volume vV,𝜈,t that the vessel has, • whether or not the vessel is unloading crude, 0≤xW,𝜈,t≤1 , • the volumetric flow rate of crude fVS,𝜈,i,t from the 𝜈 th vessel to the ith storage tank, • the volume of crude vS,i,t and the concentration of the kth component pS,i,k in it, • the volumetric flow rate of crude fVC,i,j,t and the kth component fSC,i,j,k,t from this tank to the jth charging tank, • the volume of crude vC,j,t and the volume of the kth component wC,j,k,t in it, • the volumetric flow rate of crude fVD,j,l,t and • the kth component fCD,j,l,k,t from this tank to the lth distillation unit, and • whether or not there is a transition from the jth charging tank to the j′ th one, 0≤zj,j ′ ,l,t≤1 . Fig. 1 Constant parameters considered in the scheduling problem Fig. 2 Variables considered in the scheduling problem 1 3 Optimization ofarefinery scheduling process withcolumn… Besides, we must find out several binary variables: • xU,𝜈,t and xD,𝜈,t , with value 1 when the 𝜈 th vessel starts and completes unloading at time t, respectively. • dj,t,l , with value 1 when the jth charging tank is charging into the lth distillation unit at time t. Finally, our scheduling problem is written in terms of the following cost minimization problem. Problem (SP): Find the value of the previous set of variables that optimizes with respect to tU,𝜈,tD,𝜈,vS,i,t,vC,j,t and zj,j ′ ,l,t the following MILP: and, for 𝜈=1, ..., NV,i=1, ..., NS,j=1, ..., NC,l=1, ..., ND and t=1, ..., S , subject to the following constraints over binary variables: subject to the constraints over continuous variables: (1) min CU N V ∑ 𝜈=1 (tD,𝜈−tU,𝜈)+CS N V ∑ 𝜈=1 (tU,𝜈−TA,𝜈)+CIST N S ∑ i=1 S ∑ t=1 ( vS,i,t−vS,i,t−1 2 ) +CICT NC ∑ j=1 S ∑ t=1( vC,j,t−vC,j,t−1 2 ) + S ∑ t=1 NC ∑ j=1 NC ∑ j � =1 ND ∑ l=1 CCzj,j�,l,t, (2a) S ∑ t=1 xD,𝜈,t= 1, (2b) N D ∑ l=1 dj,l,t≤ 1, (2c) N C ∑ j=1 dj,l,t= 1, (3a) tD,𝜈 ≥ TA,𝜈, (3b) t D,𝜈 −tU,𝜈≥ V V,𝜈,0 max i F VS,𝜈,i,max , (3c) tU,𝜈+1≥TD,𝜈, (3d) v V,𝜈,t=VV,𝜈,0 − N S ∑ i=1 t ∑ m=1 fVS,𝜈,i,t , J.Ossorio-Castillo, F.Pena-Brage 1 3 and subject to the constraints over mixed variables: (3e) v V,𝜈,t=VV,𝜈,0 − N S ∑ i=1 t ∑ m=1 fVS,𝜈,i,t , (3f) FVS,𝜈,i,min xW,𝜈,t≤fVS,𝜈,i,t≤FVS,𝜈,i,max xW,𝜈,t, (3g) N S ∑ i=1 S ∑ t=1 fVS,𝜈,i,t=VV,𝜈 ,0 (3h) v S,i,t=VIS,i+ N V ∑ 𝜈=1 t ∑ m=1 fVS,𝜈,i,m− N C ∑ j=1 t ∑ m=1 fSC,i,j,m , (3i) VS,i,min ≤vS,i,t≤VS,i,max, (3j) v C,j,t=VC,j,0 + N S ∑ i=1 t ∑ m=1 fSC,i,j,m= N D ∑ l=1 t ∑ m=1 fCD,j,l,m , (3k) VC,j,min ≤vC,j,t≤VC,j,max, (3l) N D ∑ l=1 S ∑ t=1 fCD,j,l,t=Dj , (3m) v C,j,t=VC,j,0 + t ∑ m=1 ( N S ∑ i=1 fSC,i,j,m− N D ∑ t=1 fCD,j,l,m) , (3n) fSC,i,j,k,t=fSC,i,j,t ⋅ pS,i,k, (3o) fCD,j,l,tpS,j,k,min ≤ fCD,j,l,k,t ≤ fCD,j,l,tpS,j,k,max, (3p) vC,j,tpS,j,k,min ≤wC,j,k,t≤vC,j,tpS,j,k,max, (3q) 0 ≤ zj,j ′ ,l,t ≤ 1, (3r) zj,j � ,l,t ≥ dj � ,l,t+dj,l,t−1−1, 1 3 Optimization ofarefinery scheduling process withcolumn… 3 Column generation In this section, we show how to decompose the previous problem—or more generally, any MILP problem of its type—into its real and its binary part. We apply the Dantzig–Wolfe decomposition (for a more detailed explanation of this decomposition, see Chvatal 1983 or Dantzig and Wolfe 1960), and then apply to it a column generation algorithm. The algorithm evolves around a master problem, namely (MP), which is updated in every iteration. Both iterations and the stopping criteria depend on two subproblems, (SP 1 ) and (SP 2 ), which are also updated during every iteration. Let us write first a generic form for the MILP we want to solve that includes the problem (SP). From now on, A matrices and 𝐚 vectors will describe the equality constraints, while B matrices and 𝐛 vectors will do the same with inequality constraints. Likewise, subindices b and r describe the coefficients associated with the binary and continuous variables respectively and 𝐜 vectors define the objective function. The 𝐱 and 𝐳 vectors designate the binary and real variables, respectively. Expressed as a linear optimization problem, we have the following objective function (4a) S ∑ t=1 tx D,𝜈,t=tD,𝜈 , (4b) S ∑ t=1 tx U,𝜈,t=tU,𝜈 , (4c) x W,𝜈,t≤ t ∑ m=1 xU,𝜈,m , (4d) x W,𝜈,t≤ S ∑ m=t xD,𝜈,t , (4e) F SC,i,j,min ( 1− N D ∑ l=1 dj,l,t ) ≤fSC,i,j,t≤FSC,i,j,max ( 1− N D ∑ l=1 dj,l,t ), (4f) FCD,j,l,mindj,l,t ≤ fCD,j,l,t ≤ FCD,j,l,maxdj,l,t. (5) min 𝐱 , 𝐳 𝐜 r⋅ 𝐱+𝐜 b⋅ 𝐳, J.Ossorio-Castillo, F.Pena-Brage 1 3 subject to the following constraints: for 𝐱 ∈ℝ m 1 ≥0 , 𝐳 ∈ℤ m 2 ≥0 , 𝐚r,𝐛r∈ℝmr , 𝐚b,𝐛b∈ℝmb , 𝐚r,𝐛r∈ℝmm , Ar,Br∈ℝn r ×mr , Ab,Bb∈ℝn b ×mb , Am,Bm∈ℝn×mm , with n=nr+nb+nm being the number of original constraints in the problem and m=mr+mb the number of original variables. We call this feasible point the first proposal, denoted by the vector 𝐱(1) for the continuous variables and the vector 𝐳(1) for the binary ones. From now on, we will denote by k and k′ the number of continuous and binary proposals stored in every iteration. These proposals will be weighted respectively in the problem (MP) with the additional continuous variables 𝜆i and 𝜇j . The master problem can be defined in its generic form as follows: Problem (MP): The previous linear and continuous problem can be solved efficiently with the help of a specialized global linear solver. The only information we need each time we solve problem (MP) is the dual solution of its constraints, which will be represented (6a) Ab𝐳=𝐚b, (6b) Bb𝐳 ≤ 𝐛b, (6c) Ar𝐱=𝐚r, (6d) Br𝐱 ≤ 𝐛r, (6e) Amr𝐱+Amb𝐳=𝐚m, (6f) Bmr𝐱+Bmb𝐳 ≤ 𝐛m, min𝜆i,𝜇j k ∑ i=1 (𝐜r⋅𝐱(i))𝜆i+ k� ∑ j=1 (𝐜b⋅𝐳(j))𝜇j, subject to k ∑ i=1 (Amr𝐱(i))𝜆i+ k� ∑ j=1 (Amb𝐳(j))𝜇j=𝐚m , k ∑ i=1 (Bmr𝐱(i))𝜆i+ k� ∑ j=1 (Bmb𝐳(j))𝜇j≤𝐛m , k ∑ i=1 𝜆i=1, k� ∑ j=1 𝜇j=1, 𝜆 i ,𝜇 j ≥0. 1 3 Optimization ofarefinery scheduling process withcolumn… from now on as the vector 𝝅 for the equality constraints and the vector 𝝆 for the inequality constraints. The value of these vectors will be updated each time we solve (MP). It is in the next step of the algorithm where we can optionally use an adiabatic quantum computer, but it is also possible to accomplish it with the help of a global nonlinear solver. Two new subproblems are defined for this step: Problem (SP 1 ): Problem (SP 2 ): As can be seen, both subproblems (SP 1 ) and (SP 2 ) depend on the dual solutions obtained from problem (MP) in order to define their objective functions. Problem (SP 1 ) is linear and continuous and can be solved again as problem (MP) with a LP solver. Problem (SP 2 ), however, have binary variables and therefore have a stronger complexity. Its resolution with the help of an quantum annealer will be explained in Sect.4. For the time being, let us just assume that we have a black box that solves it globally. Remark An Adiabatic Quantum Computing (AQC) algorithm (Farhi et al. 2001) is guaranteed to converge to the global optimum of problems but the tempering time might grow exponentially and the temperature has to be zero, while a Quantum Annealing(McGeoch and Wang 2013) process is a physical implementation of AQC (can be considered as a subcase) with finite temperature implementation and no deterministic convergence guarantees. Since our algorithm requires global convergence, an AQC would be suitable for it. However our results were obtained with a machine that implements a quantum annealing process, so we had to check that in this particular case the global solution was reached, as it is explained in Sect.5. The algorithm starts solving an instance of the original (SP 1 ) problem after dropping the objective function. This way, we will only need to get a feasible point for the problem, instead of one of its minima, a much easier achievement than its optimization counterpart. In order to complete this task we use the local nonlinear optimization solver Knitro(Byrd etal. 2006). See Fig.3 for a general overview of the algorithm. After both subproblems are solved, it is time to check the first terminating condition of the algorithm. We have to examine the value of both objective functions from (SP 1 ) and (SP 2 ). If one or both of them are less than 0, the algorithm min 𝐱 (A mr 𝝅+B mr 𝝆) ⋅ 𝐱, subject to Ar𝐱=𝐚r, Br 𝐱≤𝐛 r. min 𝐳 (A mb 𝝅+B mb 𝝆) ⋅ 𝐳, subject to Ab𝐳=𝐚b, Bb 𝐳≤𝐛 b. J.Ossorio-Castillo, F.Pena-Brage 1 3 The main problem we have found while using Qbsolv libraries or solving certain subproblems with the D-Wave machine is that finding the global optimum is not guaranteed in every case, a necessary condition for the correct execution of our algorithm. In order to avoid problems generated by these limitations, we have checked the solution obtained by Qbsolv or D-Wave with the nonlinear solver BARON (Sahinidis 2014; Tawarmalani and Sahinidis 2005), which guarantees a global solution. We have checked that in our problems the result of the algorithm was correct and, regarding the time performance, the CPU time reported by the Qbsolv tool for decomposing the original problem averages, in our case, between 500 ms and 2 s, and the estimation of the D-Wave 2X computation time for each of the subproblems (including the connection times) averages 300 ms in the P4 problem. These times are coherent with the value of 491 ms of McGeoch and Wang (2013) for an original D-wave Two and the range of 10–300 ms for a D-Wave 2000Q, presented in Chiscop etal. (2020). 6 Conclusions andfuture work In this paper we have described a possible technique for solving an optimization problem with continuous and binary variables, drawing upon a column generation scheme supported by the Dantzig–Wolfe decomposition. In this manner, thanks to breaking down the problem in its real and its binary parts, it is necessary to solve several subproblems, each of them of either a continuous or binary nature, instead of a mixed one. Although the algorithm can be used for any generic MILP, our interest resided in a certain type of optimization problems: the ones involved in the processes of a refinery, and more concretely in the ones associated with the correct scheduling of the arrival and discharge of vessels in the harbor. In problems with a different structure, this algorithm could not be that beneficial and the average number of iterations in any part of the algorithm may explode, but the complexity in the worst case remains to be seen. A possible way of palliating this is to decompose the MILP in various binary subproblems, not just one, using Dantzig–Wolfe, maintaining a certain grouping of variables and thus improving the performance of the constraint enforcement process. Other possible difficulties that may arise while solving our problems are related to the QUBO formulations. The Chimera graph present in the D-Wave chips have some serious limitations either related to the problem size or to the connectivity of its variables, proper of a technology that is currently in its early steps. In the meantime, in order to solve a QUBO problem with the D-Wave machine it is necessary to embed its graph into the Chimera graph, but this process could result in an increase in the number of calls to the quantum annealer, as explained in Sect.4. Those limitations in current hardware may alter the perception of the possibilities of using this technique, but we hope that future advancements in quantum computers will make it more competitive. Acknowledgements The authors would like to thank the Information Sciences Institute at the University of Southern California, for the possibility of using their D-Wave 2X machine, and especially Itay Hen for 1 3 Optimization ofarefinery scheduling process withcolumn… all the help and support granted during the whole process. We would also like to express our gratitude to the former employees of Repsol that supported this work. Finally, we would like to thank the referees for their helpful comments. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. References Aharonov D, van Dam W, Kempe J, Landau Z, Lloyd S, Regev O (2008) Adiabatic quantum computation is equivalent to standard quantum computation. SIAM Rev 50(4):755–787 Ajagekar AS (2020) Quantum computing for process systems optimization and data analytics. Master’s thesis, Cornell University, USA Ajagekar A, Humble T, You F (2020) Quantum computing based hybrid solution strategies for largescale discrete-continuous optimization Problems. Comput Chem Eng 132:106630 AMPL Optimization Inc. The AMPL book. Example files. http:// ampl. com/ resou rces/ theamplbook/ examp lefiles/. Accessed 22 Nov 2017 Bauckhage C, Sifa R, Wrobel S (2020) Adiabatic quantum computing for max-sum diversification. In: Proceedings of the 2020 SIAM international conference on data mining, SIAM, pp 343–351 Booth M, Reinhardt SP, Roy A (2017) Partitioning optimization problems for hybrid classical/quantum execution. 14-1006A-A. Tech. rep., D-Wave Systems, Inc Born M, Fock V (1928) Beweis des adiabatensatzes. Z Phys A Hadrons Nucl 51(3):165–180 Byrd RH, Nocedal J, Waltz RA (2006) Knitro: an integrated package for nonlinear optimization. Springer, Boston, pp 35–59 Calude CS, Dinneen MJ (2017) Solving the broadcast time problem using a D-wave quantum computer. In: Advances in unconventional computing, Springer, pp 439–453 Chiscop I, Nauta J, Veerman B, Phillipson F (2020) A hybrid solution method for the multi-service location set covering problem. In: International conference on computational science (ICCS) Chvatal V (1983) Linear programming. Macmillan, New York Dantzig GB, Wolfe P (1960) Decomposition principle for linear programs. Oper Res 8(1):101–111 D-Wave (2016) Programming with QUBOs. Release 2.3, 09-1002A-B. Tech. rep., D-Wave Systems Inc D-Wave Systems Inc. (2017) Qbsolv 2.0.4. https:// github. com/ dwave syste ms/ qbsolv Farhi E, Goldstone JS, Gutmann JL, Lundgren A, Preda D (2001) A quantum adiabatic evolution algorithm applied to random instances of an NP-complete problem. Science 292(5516):472–475 Feillet D (2010) A tutorial on column generation and branch-and-price for vehicle routing problems. 4OR Q J Oper Res 8(4):407–424 Feynman RP (1982) Simulating physics with computers. Int J Theor Phys 21(6):467–488 Furman K, El-Bakry A, Song JH (2017) Optimization in the oil and gas industry. Optim Eng 18(1):1–2 Gamrath G (2010) Generic branch-cut-and-price. Master’s thesis, Zuse Institute Berlin Karuppiah R, Furman KC, Grossmann IE (2008) Global optimization for scheduling refinery crude oil operations. Comput Chem Eng 32(11):2745–2766 Khor C, Varvarezos D (2017) Petroleum refinery optimization. Optim Eng 18(4):943–989 Lee H, Pinto JM, Grossmann IE, Park S (1996) Mixed-integer linear programming model for refinery short-term scheduling of crude oil unloading with inventory management. Ind Eng Chem Res 35(5):1630–1641 McGeoch C, Wang C (2013) Experimental evaluation of an adiabatic quantum system for combinatorial optimization. In: CF ’13: proceedings of the ACM international conference on computing frontiers, ACM, vol 23, pp 1–11 J.Ossorio-Castillo, F.Pena-Brage 1 3 Mouret S (2010) Optimal scheduling of refinery crude-oil operations. Ph.D. thesis, Carnegie Mellon University Nannicini G (2020) An introduction to quantum computing, without the physics. SIAM Rev 62(4):936–981 Sahinidis NV (2014) BARON 14.3.1: global optimization of mixed-integer nonlinear programs. User’s manual Syrichas A, Crispin A (2017) Large-scale vehicle routing problems: quantum annealing, tunings and results. Comput Oper Res 87:52–62 Tawarmalani M, Sahinidis NV (2005) A polyhedral branch-and-cut approach to global optimization. Math Program 103:225–249 Venturelli D, Marchand DJJ, Rojo G (2015) Quantum annealing implementation of job-shop scheduling. arXiv: 1506. 08479 Verstichel J, Kinable J, De Causmaecker P, Berghe GV (2015) A combinatorial benders’ decomposition for the lock scheduling problem. Comput Oper Res 54:117–128 Vyskocil T, Djidjev H (2019) Embedding equality constraints of optimization problems into a quantum annealer. Algorithms 12(4):77 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.