DisjunctiveProgramming.jl: Generalized Disjunctive Programming Models and Algorithms for JuMP

We present a Julia package, DisjunctiveProgramming.jl, that extends the functionality in JuMP.jl to allow modeling problems via logical propositions and disjunctive constraints. Such models can then be reformulated into Mixed-Integer Programs (MIPs) that can be solved with the various MIP solvers supported by JuMP. To do so, logical propositions are converted to Conjunctive Normal Form (CNF) and reformulated into equivalent algebraic constraints. Disjunctions are reformulated into mixed-integer constraints via the reformulation technique specified by the user (Big-M or Hull reformulations). The package supports reformulations for disjunctions containing linear, quadratic, and nonlinear constraints.


Introduction
The modeling of systems with discrete and continuous decisions is commonly done in algebraic form with mixed-integer programming (MIP) models.When the problems can be defined by purely linear constraints and a linear objective function, they are referred to as mixed-integer linear programs (MILP).When nonlinearities arrise in either the feasible space or the objective function, they are called mixed-integer nonlinear programs (MINLP).A more systematic approach to modeling such systems is to use Generalized Disjunctive Programming (GDP) [6,15], which generalizes the Disjunctive Programming paradigm proposed by Balas [3].GDP enables the modeling of systems from a logic-based level of abstraction that captures the fundamental rules governing such systems via algebraic constraints and logic.This formulation is useful for expressing problems in an intuitive way that relies on their logical underpinnings without needing to introduce mixed-integer constraints.GDP models are often easier to understand as related constraints are grouped into disjuncts that describe clearly defined subsets of the feasible space.The models obtained via GDP can be reformulated into the pure algebraic form best suited for the application of interest.It is also often possible to exploit the explicit logical structure of a GDP model to provide tighter relaxations than corresponding MIP models, which may improve convergence speed and robustness for solutions via advanced solution algorithms [6].Within the optimization community, there is a high volume of ongoing research that relies on GDP to formulate models for a variety of applications.Due to the combinatorial nature of system design problems, the GDP paradigm has been applied to the synthesis of complex processes and networks [21,27], the planning and optimal control of energy systems [9], and the modeling of chemical synthesis under uncertainty [8].These and numerous other applications of GDP illustrate the benefit of having a robust package for GDP that removes much of the overhead associated with developing and testing GDP models.Although packages with GDP capabilities exist for Pyomo [7] and GAMS [26], having such a package available in Julia can greatly accelerate research in optimization, where packages like JuMP.jl [10] are gaining significant traction.This paper provides background on the GDP paradigm, and the techniques for reformulating and solving such models.It then presents the package DisjunctiveProgramming.jl as an extension to JuMP.jl for creating models for optimization that follow the GDP modeling paradigm and can be solved using the vast list of supported solvers [10].A case study demonstrates the use of the package for chemical process superstructure optimization.

Generalized Disjunctive Programming
The GDP form of modeling is an abstraction that uses both algebraic and logical constraints to capture the fundamental rules governing a system.The two main reformulation strategies to transform GDP models into their equivalent MIP models are the Big-M reformulation [23,24] and the Hull reformulation [20], the latter of which yields tighter models at the expense of larger model sizes [14].

Model
The general notation for a GDP problem is given below, Here f (x) is the objective function to be minimized over the continuous variable x, g(x) represents the global constraints, h(x) represents the disjunct-specific constraints, and Y is the Boolean variable governing each disjunction.In this notation, there are k disjunctions with i disjuncts in each.Constraints h ik (x) ≤ 0, are applied only if the Boolean indicator variable for the respective disjunct, Y ik , is denoted as being active (i.e., true) [6].The set of logical constraints, Ω(Y ), describe the logical relationships between Proceedings of JuliaCon 1(1), 2021 Fig. 1.Feasible solution space for example disjunction the selections of indicator variables.These can take the form of logical propositions or constraint programming expressions.
In the case of a linear objective function and constraint set, the GDP model can be written as,

Solution Technique: Reformulation to Mixed-Integer Program
The simplest example of a linear GDP system is given below, where Y i is a Boolean indicator variable that enforces the constraints in the disjunct (Ax ≤ b or Cx ≤ d) when true, For visualization purposes and without loss of generality, the simple linear example is used to illustrate the Big-M and Hull reformulations.Figure 1 illustrates the feasible space of a simple linear GDP with one disjunction and two continuous variables, x 1 and x 2 .The rectangle on the left is described by the constraints in the left disjunct, Ax ≤ b.The rectangle on the right is defined by the constraints in the right disjunct, Cx ≤ d.The non-overlapping nature of the two regions is supported by the exclusive-OR relationship above.

Big-M Reformulation.
The Big-M reformulation for this problem is given below, where M is a sufficiently large scalar that makes the particular constraint redundant when its indicator variable is not selected (i.e., y i = 0).Note that the Boolean variables, Y i , are replaced by binary variables, y i .When the integrality con-Fig.2. Relaxed solution space using Big-M Reformulation straint y 1 , y 2 ∈ {0, 1} is relaxed to 0 ≤ x 1 , x 2 ≤ 1, the resulting feasible region can be visualized by projecting the relaxed model onto the x 1 , x 2 plane.This results in the region encapsulated by the dashed line in Figure 2. It should be noted that the relaxed feasible region is not as tight as possible around the original feasible solution space.The choice of the large M value determines the tightness of this relaxation, and the minimal value of M for the optimal relaxation can be found through interval arithmetic when the model is linear.For nonlinear models, the tightest M can be obtained by solving the maximization problem {max h ik (x) : x ∈ X}.An alternate method for tight Big-M relaxations is given in [24].
2.2.2 Hull Reformulation.The Hull reformulation is given below.This formulation requires lifting the model to a higherdimensional space.When projected to the original space, the continuous relaxation of the model is tighter than its Big-M equivalent [15].The reformulation relaxation can be visualized by the region encapsulated by the dashed line in Figure 3.Note that this reformulation provides a tighter relaxation than the Big-M reformulation in Figure 2. Also note that describing the geometry of this relaxation is more complex than the Big-M relaxation, which is made possible by the increased number of constraints and variables in the model.
Once the logic propositions are converted to CNF, each clause can be converted into an algebraic constraint with the following equivalence, where the set represents the subset boolean variables present in the clause, and the set J represents the subset of boolean variables present in the clause in negated form, Alternate approaches exist for converting propositional logic statements into CNF, which involve preserving clause satisfiability rather than clause equivalence.These approaches prevent exponential size increase in clauses and yield logically consistent results [18].

Constraint
Programming.Selection constraints analogous to those used in Constraint Programming (CP) can also be included in Ω(Y ).These constraints are of the form "allow exactly n elements in a list of Boolean variables to be true."This type of constraint overcomes the limitations of the Boolean exclusive-OR (∨) operator, which can only enforce that an odd number of elements in a list of Booleans be true.Other CP-like constraints can be ob-tained by replacing "exactly" with "at most" or "at least".These constraints are reformulated as follows, Exclusive-OR constraints as the one given in Eq. (??) are more generally modeled as exactly(1, {Y 1 , Y 2 }).

Other Solution Techniques
2.4.1 Disjunctive branch and bound.The disjunctive branch and bound method closely mirrors the standard branch and bound approach for the solution of mixed-integer programming problems [14].A search tree is initialized by solving the continuous relaxation of the Big-M or Hull reformulation of the original GDP to obtain a lower bound on the optimum.Branching is then done on the disjunction with an indicator binary variable closest to 1. Two nodes are created at this point: one where the respective indicator Boolean variable is fixed to true (the disjunct is enforced) and another where it is fixed to f alse (the disjunct is removed from the disjunction).Each node is reformulated and solved to obtain a candidate lower bound.If the solution to a node results in a feasible solution that satisfies all integrality constraints, the solution is an upper bound on the optimum.Any non-integral solutions that exceed an upper bound are pruned from the search tree.The process is repeated until the lower and upper bounds are within the desired tolerance.

Logic-based outer approximation.
Logic-based outer approximation is another algorithm which mirrors a standard technique for solving mixed-integer nonlinear programming problems [11].This approach starts by identifying a set of reduced Non-Linear Programming (NLP) sub-problems obtained by fixing Boolean variables in the different disjunctions such that each disjunct is selected at least once across the set of sub-problems (set covering step).Each sub-problem is solved to obtain an upper bound and a feasible point, about which the objective and constraints of the original GDP are linearized, and solve the resulting problem (via direct reformulation to MILP or via disjunctive branch and bound) to find a lower bound.If the lower and upper bound solutions have not converged, the Boolean variables from the previous solution are fixed and the resulting NLP is solved to find a potentially tighter upper bound solution.The procedure is repeated until convergence is obtained.

Hybrid cutting planes.
The cutting planes method is an algorithm for tightening the relaxed solution space of a problem reformulated with Big-M before solving it by adding additional constraints which remove parts of the relaxed space that are disjoint from the actual feasible solution space.These "cuts" to the relaxed solution space are derived from the tighter, hull relaxation of the problem.This algorithm provides a middle-ground for the tradeoff between the complexity and corresponding computational expense of the Hull reformulation with the less tight Big-M reformulation.[25].

DisjunctiveProgramming.jl
The following section describes the features of the DisjunctiveProgramming.jl package and illustrates its syntax with an example from the chemical processing industry for superstructure optimization.The use of nested disjunctions is also shown.

Features
DisjunctiveProgramming.jl allows for defining JuMP models with disjunctions that are directly reformulated via Big-M or Hull methods via the @disjunction macro or add_disjunction!function.Big-M values can be specified either for the entire disjunction, for each disjunct, or for each constraint in each disjunct.Alternately, if the constraints are linear, the code can use the variable bounds to perform interval arithmetic on each constraint to determine the tightest possible M value to use [2].For nonlinear GDP constraints, the epsilon approximation formulation for the perspective function in the Hull reformulation is used [12].Users can specify an epsilon tolerance value to use.Perspective functions are generated by relying on manipulation of symbolic expressions via Symbolics.jl[13].Logical propositions can be added to JuMP models using expressions involving the disjunction indicator variables and the standard Boolean operators (⇔, ⇒, ∨, ∧, and ¬) in an @proposition macro.These are automatically converted into CNF and added as integer algebraic constraints to the model.The constraint programming constraints can also be added using the choose!function.The expressions are also automatically reformulated into integer algebraic constraints.Nesting of disjunctions is also supported.

Example
To illustrate the syntax in DisjunctiveProgramming.jl (Version 0.4.1),consider the simple superstructure optimization problem for the chemical process given in Figure 4.In this problem a chemical plant with two candidate reactor technologies (R 1 and R 2 ) must be designed.If the second reactor technology is chosen, a separation system must also be installed, for which two separation technologies (S 1 and S 2 ) are available.The GDP model seeks to maximize the final product flow (F 7 ), while discounting for reactor and separator installation costs (C R and C S , respectively), subject to the nested disjunction and the global mass balances on the stream flows F i .The system variables are the flows on each stream i (F i ) and the installation costs, with their respective bounds given in problem formulation.The fixed cost and process yield parameters are given by γ and β, respectively.
The above system can be modeled and reformulated via the Big-M reformulation using DisjunctiveProgramming.jl.The resulting JuMP model is then optimized using the HiGHS open-source MILP solver [17] as shown below.

¦ ¥
(2) Define the inner (nested) disjunction for the separation technologies in the superstructure using the @disjunction macro.§ ¤ # d e f i n e c o n s t r a i n t s in l e f t YS d i s j u n c t Y S 1 _ d i s j u n c t = D i s j u n c t C o n s t r a i n t ( YS [ 1 ] ) @ c o n s t r a i n t ( m , F [ 5 ] = = β [ : S1 ] * F [ 4 ] , Y S 1 _ d i s j u n c t ) @ c o n s t r a i n t ( m , CS = = γ [ : S1 ] , Y S 1 _ d i s j u n c t ) # d e f i n e c o n s t r a i n t s in r i g h t YS d i s j u n c t Y S 2 _ d i s j u n c t = D i s j u n c t C o n s t r a i n t ( YS [ 2 ] ) @ c o n s t r a i n t ( m , F [ 5 ] = = β [ : S2 ] * F [ 4 ] , Y S 2 _ d i s j u n c t ) @ c o n s t r a i n t ( m , CS = = γ [ : S2 ] , Y S 2 _ d i s j u n c t ) # d e f i n e d i s j u n c t i o n ( s p e c i f y p a r e n t d i s j u n c t ) @ d i s j u n c t i o n ( m , YS , Y S 2 _ d i s j u n c t )

¦ ¥
(3) Define the outer disjunctions for the reaction pathway selection.§ ¤ # d e f i n e c o n s t r a i n t s in l e f t YR d i s j u n c t Y R 1 _ d i s j u n c t = D i s j u n c t C o n s t r a i n t ( YR [ 1 ] ) @ c o n s t r a i n t ( m , F Add the selection logical constraints using the choose!function.The first constraint enforces that only one reactor is selected The second constraint enforces that the separation system be defined only if the second reactor § ¤ @ c o n s t r a i n t ( m , YR in E x a c t l y ( 1 ) ) @ c o n s t r a i n t ( m , YS in E x a c t l y ( YR [ 2 ] ) ) ¦ ¥

Future Work
Since the package is currently limited to reformulating GDP models into Mixed-Integer Programming JuMP models, future work includes developing logic-based solvers to allow optimizing GDP models directly.This will include GDP solvers that implement the other solution techniques described in Section 2.4.Future work also involves extending the list of suported reformulation techniques to include the P-Split [19] and True-False [1] reformulations.

Related Work
The popular Python package Pyomo [5,16] is widely used for optimization development and includes an extension for generalized disjunctive programming [7].GAMS [4] is a widely used optimization modeling language with support for GDP under the GAMS EMP solver that uses LogMIP [26].Research is also being conducted to integrate modern process simulation technology, such as Aspen, within the GDP paradigm [22].

Conclusion
DisjunctiveProgramming.jl is an extension to JuMP for creating models for optimization that are formulated according to the generalized disjunctive programming paradigm.The package provides several options for reformulations including the Big-M and Hull relaxations.This package can be used to model problems, reformulate them, and optimize them using existing mathematical programming infrastructure in JuMP.This can be useful for industrial and academic applications of GDP, such as superstructure optimization.The capabilities of this package allow for this modeling paradigm to be exploited using Julia's efficient dynamicallytyped systems for rapid development, building, and testing of optimization models.

Fig. 3 .
Fig. 3. Relaxed solution space using Hull Reformulation 2.3 Logic constraint reformulation 2.3.1 Propositional Logic.The logic propositions within the set of decision selection relationships, Ω, must be converted into conjunctive normal form (CNF) to enable reformulating a GDP model as a MIP model.This means that each clause within the set of propositions must be formulated into a conjunction of disjunctions.This process can be accomplished by following the simplifying rules of propositional logic (De Morgan's laws).For boolean variables A, B, and C the following rules are used for converting to CNF (in the order given below),

( 1 )
Create the JuMP model and define the model variables and global constraints (mass balances).§ ¤ u s i n g D i s j u n c t i v e P r o g r a m m i n g , H i G H S # c r e a t e m o d e l m = G D P M o d e l ( H i G H S .O p t i m i z e r ) # a d d v a r i a b l e s to m o d e l @ v a r i a b l e ( m , 0 < = F [ i = 1 : 7 ] < = 10 ) @ v a r i a b l e ( m , 0 < = CS < = C S m a x ) @ v a r i a b l e ( m , C R m i n < = CR < = C R m a x ) # a d d l o g i c a l v a r i a b l e s to m o d e l @ v a r i a b l e ( m , YR [ 1 : 2 ] , L o g i c a l V a r i a b l e ) @ v a r i a b l e ( m , YS [ 1 : 2 ] , L o g i c a l V a r i a b l e ) # a d d g l o b a l c o n s t r a i n t s to m o d e l @ c o n s t r a i n t s ( m

( 5 )
Add the objective function and optimize.§ ¤ @ o b j e c t i v e ( m , M a x , F [ 7 ] -CS -CR ) o p t i m i z e ! ( m ) ¦ ¥