{'fraction_non_alphanumeric': 0.045459997601055534, 'fraction_numerical': 0.01206669065611131, 'mean_word_length': 3.9922754491017964, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 34, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the Structure of Decision Diagram-Representable Mixed Integer Programs with Application to Unit Commitment

Over the past decade, decision diagrams (DDs) have been used to model and solve integer programming and combinatorial optimization problems. Despite successful performance of DDs in solving various discrete optimization problems, their extension to model mixed integer programs (MIPs) such as those appearing in energy applications has been lacking. More broadly, the question which problem structures admit a DD representation is still open in the DDs community. In this paper, we address this question by introducing a geometric decomposition framework based on rectangular formations that provides both necessary and sufficient conditions for a general MIP to be representable by DDs. As a special case, we show that any bounded mixed integer linear program admits a DD representation through a specialized Benders decomposition technique. The resulting DD encodes both integer and continuous variables, and therefore is amenable to the addition of feasibility and optimality cuts through refinement procedures. As an application for this framework, we develop a novel solution methodology for the unit commitment problem (UCP) in the wholesale electricity market. Computational experiments conducted on a stochastic variant of the UCP show a significant improvement of the solution time for the proposed method when compared to the outcome of modern solvers.


Introduction
Improved methods for scheduling and dispatching resources are necessary to accommodate the increases in renewable generation, distributed energy resources and storage in the electric grid. At the core of electric grid lies the Unit Commitment Problem (UCP), which determines the thermal unit schedules and generation dispatch in wholesale electricity markets. Unfortunately, the UCP is computationally difficult, making it challenging for available software despite all advances in computing power and technology. For this reason, innovative methodologies that produce optimal solutions to the UCP satisfying realistic yet complicated constraints are needed to manage the increasing amounts of variable renewable generation and distributed energy resources in the electric grid.
Over the past decade, a powerful alternative to traditional solution techniques was developed based on a new paradigm called Decision Diagrams (DDs), in which discrete optimization problems are reconstructed as network models with special structure. The network structure is exploited to speed up the solution process and improve the solution quality, especially for formulations that are poorly handled by conventional solvers. However, since DDs have been designed for modeling discrete problems, they have never been used to solve optimization problems that appear in the electric grid applications, as they naturally contain both discrete and continuous variables. In this paper, we generalize the concept of DDs to mixed integer programs (MIPs), and thereby, establish a novel framework that can be applied to the UCP with the purpose of improving the solution time and quality obtained from existing techniques.

Background on Decision Diagrams
DDs were introduced in 2006 (Hadzić and Hooker 2006) as a solution method for discrete optimization and combinatorial problems. DDs' principal idea is to model the solutions of the underlying discrete set through a directed acyclic graphical structure with a root and a terminal node. Each path of the graph from the root to the terminal node encodes a solution of the problem where arcs are labeled by values of decision variables on their relative domain. Structurally, DDs resemble branch-and-bound trees and dynamic programming graphs. The key difference, however, is that the size of DDs can be controlled by a width limit factor that mitigates the exponential growth rate intrinsic to branch-and-bound and dynamic programming. This control factor is attributed to a node-merging concept where nodes of the DD are merged to reduce its size at the price of expanding the solution set, thereby yielding a discrete relaxation. By exploiting the concept of relaxed DDs, a specialized branch-and-bound method is designed to successively refine relaxed DDs until proving optimality.
This innovative approach to restructure the graph of the solution set makes DDs a competitive alternative to traditional divide-and-conquer solution techniques such as branch-and-bound and dynamic programming. Various studies in the past decade have been devoted to showing remarkable improvements of solution time and quality achieved by DDs when compared to the outcome of modern solvers. DDs have found their way to a broad array of application areas, from healthcare to supply chain management to finance (Bergman and Cire 2018). Other directions of application include cutting plane theory (Davarnia and van Hoeve 2020), multi-objective optimization , post-optimality analysis (Serra and Hooker 2020), and integrated branch-and-bound (Gonzalez et al. 2020).
While DDs provide a powerful optimization tool, their applicability domain is limited to discrete problems. This limitation is due to a structural requirement that DDs contain a finite number of arcs whose labels represent variable values on their relative domain. Consequently, a successful extension of DDs to model MIPs has been lacking in the literature. Davarnia (2021) took the initiative in this direction by introducing a new concept called arc-reduction that enables building relaxed DDs for continuous nonlinear programs. Such a DD is used in a cut-generating oracle to obtain linear valid inequalities for the continuous feasible region of the underlying problem. When embedded inside an outer-approximation scheme, the resulting cutting planes can improve the bounds obtained from classical methods. This cut-generating framework can be viewed as an interface between DDs and traditional cutting plane methods. Such an intermediary role, however, does not utilize the full potential of DDs achieved through specialized branch-and-bound and refinement methods. In the present paper, we develop a framework that directly models MIPs through DDs.

Background on Unit Commitment Problem
The UCP is strongly NP-hard as shown in Bendotti et al. (2019). Over the past three decades, numerous exact and heuristic approaches have been proposed to solve different variants of the UCP and its extensions under both deterministic and stochastic settings, resulting in a rich literature.
Examples to solve stochastic/robust UCP include Benders decomposition (Li et al. 2007), two-stage stochastic programming (Blanco and Morales 2017), heuristic scenario reduction methods (Feng and Ryan 2016), two-stage adaptive robust optimization (Bertsimas et al. 2012), and algorithms based on constraint generation (Lorca et al. 2016). We encourage the interested reader to consult a survey by van Ackooij et al. (2018) on solution methods for the UCP under uncertainty.
Meanwhile, examples to solve deterministic variants include priority listing methods (Lee and Feng 1992), Tabu search (Rajan et al. 2003), Genetic algorithms (Maifeld and Sheble 1996), dynamic programming (Pang et al. 1981, Frangioni andGentile 2006), simulated annealing (Mantawy et al. 1998), Lagrangian relaxation (Baldick 1995), fuzzy systems (Saneifard et al. 1997), and MIP formulations. We refer the interested reader to Padhy (2004) for a detailed account on deterministic UCP. Other solution approaches include polyhedral analysis of UCP variants; see Rajan et al. (2005), Queyranne and Wolsey (2017), and Bendotti et al. (2018) for examples. Garver (1962) is among the pioneers to model the initial variants of the UCP by a MIP that uses three sets of binary variables: (i) on/off variables, (ii) start-up variables, and (iii) shut-down variables. Since then, many efforts have been made to improve MIP formulations either by making them easier to solve through reducing the number of variables, constraints and nonzeros, or by proposing stronger formulations with tighter LP relaxation. Ostrowski et al. (2011) propose a new MIP formulation with three sets of binary variables and show that the computational results are superior to those obtained by using fewer sets of 0-1 variables. Meanwhile, as explained by the authors, one can generate strong valid inequalities when using all three sets of 0-1 variables, which in turn leads to tighter LP relaxations. Subsequently, Morales-España et al. (2013) propose a stronger MIP formulation than those of Carrión and Arroyo (2006) and Ostrowski et al. (2011) with a fewer number of constraints and nonzeros. Interested reader is referred to Wu (2011) and Frangioni et al. (2008) for other examples of MIP models. In addition, a comprehensive analysis and comparison of different MIP formulations are provided by Knueven et al. (2020b).
Due to the substantial cost savings resulted from high quality solutions of the unit commitment, a host of UCP formulations and solution techniques have been proposed to solve real-world instances; see Bertsimas et al. (2012), Atakan et al. (2017), Franz et al. (2020), and Knueven et al. (2020a), for successful implementations.

Contributions
As discussed earlier, the existing DD-based solution methods are limited to discrete problems only.
In this paper, we introduce a framework that generalizes the concept of DDs to solve MIPs directly, which allows for taking advantage of the full arsenal of branch-and-bound and refinement techniques available for DDs. In particular, we propose a geometric concept based on decomposing a MIP solution set into rectangular formations of certain properties, which provides both necessary and sufficient conditions for a general optimization problem to be representable by DDs. The significance of this contribution is two-fold: (i) it notably expands the applicability domain of DDs by lifting the integrality barrier; and (ii) it addresses a fundamental open question in the community on which problem structures are amenable to a DD representation. As a consequence, we show that a bounded mixed integer linear program admits a DD representation when reformulated through a specialized Benders decomposition. This result opens a new vein of applications for DDs to solve a broader class of optimization problems. As an evidence for such applications, we develop a novel DD-based solution method to solve the UCP in the energy sector. Being the first DD-based solution technique for the UCP, our results exhibit a great potential in solving this challenging problem class more efficiently.
The remainder of the paper is organized as follows. We develop a rectangular decomposition technique for MIPs and establish its connection to constructing DDs with both integer and continuous variables in Section 2. We transition from the concept of rectangular decomposition to the UCP application through a specialized Benders decomposition technique presented in Section 3.
Section 4 is devoted to solving the UCP using the developed DD-based framework. Concluding remarks are given in Section 5.
Due to space restrictions, proofs are omitted from the main part of the paper and can be found in Appendix A. Further, a technical comparison of our framework on the integration of Benders decomposition and DDs with those in the literature is given in Appendix B.

Decision Diagrams Representation
In this section, we introduce the concept of rectangular decomposition which bridges a general MIP and DDs. But first, we give a brief overview on the basics of DDs for optimization. A comprehensive review can be found in .

Overview on Decision Diagrams
Define N := {1, . . . , n}. We refer to a DD by D = (U, A, l(·)), where U denotes the set of nodes and A represents the set of arcs, and function l : A → R indicates the label of arcs in A. The multi-graph induced by D is composed of n arc layers A 1 , A 2 , . . . , A n , and n + 1 node layers U 1 , U 2 , . . . , U n+1 .
Node layer U 1 contains a single root node r, and node layer U n+1 contains a single terminal node t. We define the width of a DD as the maximum number of nodes at any node layer, that is max{|U i | : i ∈ {1, . . . , n + 1}}.
DD D represents a set of points of the form x = (x 1 , . . . , x n ) with the following characteristics.
The label l(a) of each arc a ∈ A j , for j ∈ N , represents the value of x j . Each arc-sequence (path) from r to t encodes a specific value assignment to x associated with the labels of the arcs on the path. The collection of points encoded by all paths of D is referred to as the solution set of D, which is denoted by Sol(D).
Consider a bounded integer set P ⊆ Z n . Set P can be expressed by a DD D whose collection of all r-t paths encodes the solutions of P, i.e., P = Sol(D). Using this definition, we can model discrete problems with DDs. Consider a bounded integer program z * = max {f (x) | x ∈ P} where f : R n → R and P ⊆ Z n . To model the above integer program with a DD, we first construct a weighted DD, denoted by [D|w(·)], where (i) D represents an exact DD encoding solutions of P, and (ii) w(·) : A → R is a weight function associated with arcs of D so that for each r-t path P = (a 1 , . . . , a n ), its weight w(P) := n i=1 w(a i ) is equal to f (x P ), the objective value of the integral solution corresponding to P. Construction of the weight function is immediate when as we can simply define w(a i ) = f i (l(a i )) for a ∈ A i and i ∈ N . The above definition implies that z * is equal to the weight of a longest path from r to t in the corresponding weighted acyclic graph.

Rectangular Decomposition
DDs are often described as a union of arc-sequences (paths) from r to t. Since arcs on the path encode certain value assignments to variables, DDs are composed of a finite number of solution points. The key to extend DDs to model continuous variables is to view them as a union of nodesequences from r to t. From this perspective, multiple arcs with different label values can be considered simultaneously between two consecutive nodes. We will show later that, in an extreme case, the arcs between two nodes can virtually span an entire continuous interval between two values in the domain of variables. As a result, the node-sequence will encode a hyper-rectangle, as opposed to a single point encoded by the arc-sequence. This analogy lays the foundation for a decomposition technique, which we refer to as rectangular decomposition, that restructures the feasible region through a collection of hyper-rectangles with specific properties. This decomposition framework plays a key role in our ability to represent general sets via DDs.
For any I ⊆ N , let x I represent coordinates of x ∈ R n whose index belongs to I. Given a set P ⊆ R n , conv(P ) describes the convex hull of P , proj x I (P ) represents the projection of P onto the space of variables x I , and dim P indicates the dimension of conv(P ). Given a function f (x) : R n → R and a variable subset When we refer to extreme points of a bounded set P , denoted by X (P ), we mean extreme points of conv(P ).
Theorem 1. Consider a compact set P ⊆ R n , and select I ⊆ N . Assume that there exists a finite collection of compact sets P j I for j ∈ J, where J is an index set, such that (i) dim proj x i (P j I ) = 0, for i ∈ N \ I and j ∈ J, i.e., coordinate x i is fixed.
(ii) j∈J X (P j I ) ⊆ P ⊆ j∈J P j I .
(iii) For each j ∈ J, there exists a finite collection of hyper-rectangles of the form We say that set P admits a rectangular decomposition w.r.t. I, if there exists a finite collection of compact sets P j I that satisfy conditions (i)-(iii) of Theorem 1. We next illustrate this decomposition concept in an example.
Example 1. Consider a compact mixed integer set Set P is shown in Figure 1. Select I = {1, 2}. To use the result of Theorem 1, we first introduce sets It is easy to verify that conditions (i) and (ii) of Theorem 1 hold for these sets as x 3 is fixed in sets P 1 I and P 2 I , and {(0, 0, 0), (1, 0, 0), (0, 2, 0), (1, 2, 0), (1, 1, 1), (1, 2, 1), (2, 0, 1), (2, 1, 1)} ⊆ P ⊆ In view of Theorem 1, two levels of decomposition are performed. First, P is decomposed into sets P j I through fixing variables not in I. Second, each set P j I is decomposed into hyper-rectangles through a convex hull description. The second level can be viewed as a generalization of the standard extreme point decomposition theorem for maximizing convex functions (Rockafeller 1970), as the hyper-rectangles can be chosen to be the extreme points. We note, however, that the addition of the first level is necessary to account for the non-convex coordinates of the objective function. does not provide the desired equivalence as max{−x 2 2 | x ∈ P} = max{−x 2 2 | x ∈ X (P)}. Further, we will discuss later that the rectangular volume of the decomposition is key to build efficient DDs as compared to the case for extreme points. Next, we present a partial converse of Theorem 1 as a sufficient condition for a set to admit rectangular decomposition.
Proposition 1. Consider a compact set P ⊆ R n , and select I ⊆ N . Assume that there exists a finite set Q ⊆ R n such that max{f (x) | x ∈ P} = max{f (x) | x ∈ Q} for every function f (x) that is convex in x I . Then, P admits a rectangular decomposition w.r.t. I.
The above rectangular decomposition method has two important consequences in DDs context that will be demonstrated in the next two subsections.

Equivalence Class
In this section, we establish an equivalence relation between DDs that produce the same optimal value for a given objective function. This result sets the stage for modeling continuous variables via DDs.
Lemma 1. Consider a compact set P ⊆ R n , and select I ⊆ N . Assume that there exists a finite collection of hyper-rectangles of the form R j I = n k=1 [l j k , u j k ] for j ∈ J, where J is an index set, such that (i) l j i = u j i for all i ∈ N \ I and j ∈ J, The difference between the result of Lemma 1 and that of Theorem 1 stems from the levels of decomposition involved in the process. In Lemma 1, a set is directly decomposed into rectangular components, whereas in Theorem 1, a preliminary level of decomposition is performed. This additional level, through definition of sets P j I , makes it possible to decompose sets that cannot be directly represented through union of hyper-rectangles. For instance, it is easy to verify that set P of Example 1 cannot be directly decomposed into hyper-rectangles prescribed in Lemma 1, whereas it admits rectangular decomposition through the use of Theorem 1.
In view of Lemma 1, we say that the hyper-rectangles R j I that satisfy the lemma assumptions form an equivalence class w.r.t. I, as they share the optimal value of the objective function over any set enveloped between R j I and their extreme points. We next present a similar equivalence class for DD formulations. For notational convenience, we use label value l a for an arc a ∈ A as a shorthand for l(a). The DDs considered here are not limited by the unique arc-label rule, i.e., they can contain several arcs with equal label values at a layer. For uniformity of arguments, we refer to a virtual DD as one that contains, between two nodes, a collection of arcs with label values that span a closed continuous interval; see Davarnia (2021) for an example.
Consider a DD D = (U, A, l(·)). For any pair (u, v) of nodes of D, define A(u, v) to be the set of all arcs of D directed from u to v. Further, define l max (u,v) (resp. l min (u,v) ) to be the maximum (resp. minimum) label of the arcs in a non-empty set A(u, v).
(i) LetD be a DD constructed from D with a difference that, for any node pair (u, v) ∈ U i × U i+1 and any i ∈ I, only the arcs with label values l min (u,v) and l max (u,v) are maintained and the rest are removed.
(ii) LetD be a virtual DD constructed from D with a difference that, for any node pair (u, v) ∈ U i × U i+1 and any i ∈ I, the collection of arcs with label values spanning the interval [l min (u,v) , l max (u,v) ] is added between u and v.
Proposition 2 defines an equivalence class for DDs. In words, when there are multiple arcs between two consecutive nodes of a DD, as long as the objective function is convex in the variable corresponding to that arc layer, we can add or remove the middle arcs, and yet preserve the optimal value. This equivalence property enables us to model continuous intervals by translating virtual DDs into regular DDs. This property can also be useful from a computational perspective as the size of the DDs can be reduced through removing unnecessary arcs or through adding arcs that provide more efficient merging possibilities.

DD-Representable Sets
An important and fundamental question in DDs community concerns the characteristics of optimization problems that can be represented through a DD formulation.
Definition 1. We say that a compact set P ⊆ R n is DD-representable w.r.t. an index set I, if there exists a DD D such that max{f (x)|x ∈ P} = max{f (x)|x ∈ Sol(D)} for every function f (x) that is convex in x I .
In the above definition, the requirement that a set and its associated DD must have the same optimal value for every objective function-under appropriate convexity assumptions-is critical for DD-representability. There are several DD-based procedures, such as Lagrangian relaxation methods (Bergman et al. 2015), that use DDs with varying objective functions in different stages of the procedure, showing the necessity of properties that guarantee the conservation of optimal values throughout stages.
The question of whether a given set is DD-representable is straightforward when the set is bounded and discrete. In this case, every point of the set can be encoded by a unique r-t path of a DD, leading to a finite width limit. This argument, however, does not hold when the set contains continuous variables, as all solutions cannot be fully encoded by finitely many paths in a DD. In light of the rectangular decomposition technique introduced above, we next give necessary and sufficient conditions for a mixed integer set to be DD-representable.
Corollary 1. Consider a compact set P ⊆ R n , and select I ⊆ N . The set P is DD-representable w.r.t. I if and only if it admits a rectangular decomposition w.r.t. I.
As mentioned above, the conditions of Corollary 1 are immediately satisfied for bounded discrete sets. For mixed integer sets, there is an important class that also satisfies the conditions of Corollary 1 as given next.
Corollary 2. Let Q ⊆ R n+1 be a compact set and define the mixed integer set P The class of mixed integer sets that contain a single continuous variable plays an important role in optimization as it can be viewed as the bridge between a pure discrete case and a general mixed integer case; see e.g. the derivation of mixed integer rounding inequalities that is built upon this analogy (Nemhauser andWolsey 1990, Conforti et al. 2014b). The ability to model this class by DDs, therefore, opens pathways to a wide range of new application domains. The approach is to represent the continuous components of a general mixed integer set by a new variable and develop its corresponding DD, while accounting for the relation between the continuous components and the substitute variable separately. The Benders decomposition technique provides an ideal framework for such an approach. We will establish this framework in Section 3 and apply it to the unit commitment problem in Section 4.
In view of Corollary 2, we remark that when a set involves multiple continuous variables, the existence of a rectangular decomposition is not guaranteed. For example, the unit disk in R 2 described by x 2 1 + x 2 2 ≤ 1 does not admit a rectangular decomposition; hence it is not DD-representable.
We conclude this section by presenting a key role of the rectangular decomposition technique in determining the DD size, as the main factor in designing efficient DD-based solution methods. This implication follows from the fact that the rectangular decomposition of a set is not unique, and it can be achieved in different ways. For each decomposition, the resulting hyper-rectangles form an equivalence class for DDs representing them. This variety can be exploited to build DDs with smaller sizes that are computationally more efficient, as illustrated in the next example.
Example 2. Consider the set defined by proj (x 1 ,x 2 ) (P) of Example 1; see Figure 2. The goal is 2}. An immediate rectangular decomposition is obtained through rectangles R j 1 that represent extreme points {(0, 0), (0, 2), (1, 2), (2, 0), (2, 1)} of P. The DD D 1 that models these extreme points forms an equivalence class. The minimum width of D 1 (with the natural ordering of variables in layers) is three, as shown in Figure  It follows from Example 2 that to reduce the size of the DD representation of a given set, one may seek the solution to the combinatorial problem that finds the minimum number of hyperrectangles that achieve a rectangular decomposition of the set. This observation implies that a decomposition that contains higher-dimensional hyper-rectangles is generally preferred to one that is simply composed of extreme points (i.e., zero-dimensional rectangles) of the set, as the former can be modeled via fewer node-sequences in a DD.

Figure 3
The impact of equivalence class on DDs width. Numbers next to arcs represent the labels.

Benders Decomposition
As an application of the rectangular decomposition technique developed in Section 2.2, we next present a specialized Benders decomposition (BD) framework that uses DDs to solve MIPs.

Consider a bounded MIP
vectors. We use BD to define the master problem M = max{ax + z | (x, z) ∈ P × [l 0 , u 0 ]} where P ⊆ Z n is the projection of the feasible region of H onto x-space, and the bounds on z are induced from the boundedness of H. Subproblems are defined as S(x) = max{by | By ≤ c − Ax} for a given x. At each iteration of the algorithm, the output of the subproblems is either a feasibility cut of the form αx ≤ α 0 or an optimality cut of the form z + αx ≤ α 0 , where α is a row vector of matching dimension. These cuts are added to the master problem and the problem is resolved.
It follows from Corollary 2 that M admits a DD representation, forming an equivalence class.
Such a DD contains n layers corresponding to integer variables x and one (last) layer corresponding to the continuous variable z with arc labels representing a lower and an upper bound on this variable. Using this DD, we can find its longest r-t path, solve the subproblems for the solution encoding that path, and refine the DD with respect to the generated optimality and feasibility cuts.
The refinement operation separates the solutions violating the cuts from the DD solution set; see  for a detailed account on refinement techniques in DD. The main difficulty, however, is that the size of a DD that gives an exact representation of M grows exponentially with the problem size. It is then imperative to employ the notion of restricted and relaxed DDs in our algorithm to increase efficiency. The idea is to design DDs with a limited width that represent a restriction (resp. relaxation) of the underlying problem, and thereby providing a primal (resp. dual) bound. The construction of a restricted DD is straightforward, as it can be achieved by selecting any subset of the r-t paths of the exact variant that satisfies the width limit. The construction of a relaxed DD, however, requires a careful manipulation of the DD structure through a so-called node-merging operation in such a way that the resulting DD with a smaller width contains all feasible r-t paths of the exact variant, with a possible addition of some infeasible paths. Such a construction exploits the combinatorial structure of M, and hence is problem-specific. We give an instance for designing relaxed DDs for the unit commitment application in Section 4.
Algorithm 1 presents the full DD-based framework to solve H through BD. In this approach, referred to as DD-BD, we use the following definitions. Let C be the set of optimality and feasibility cuts, and define F C to be the feasible region described by constraints of C. Consider x = (x 1 , . . . ,x |x| ) to be a partial assignment of size |x| to variables x 1 , . . . , . . , |x|} to be the restriction of the master problem M after adding the cuts in C and fixing the partial assignmentx. We denote the case with an empty partial assignment and empty constraint set as input by M ∅ (∅) = M. We assume that oracles to build relaxed and restricted DDs for the master problem M are available. Using these oracles, one can simply construct relaxed and restricted DDs for M C (x) by fixing the path associated with the partial assignmentx and refining the DD with respect to the cuts in C. Given a DD D with n layers, we denote by EXACT(D) the last exact node layer of D, i.e., the last node layer that contains no merged nodes; see  for properties of the exact cut set operator.
Next, we show the finiteness and correctness of Algorithm 1, followed by some implementation remarks. solve S(x) to obtain feasibility/optimality cuts C and optimal value ρ 20 update C ← C ∪ C, refine D w.r.t C, and update its longest r-t path solution (x, z) and length w 21 until z = ρ;

22
For each node u ∈ EXACT(D), add the partial assignment encoding a longest r-u path of D toX Remark 2. The addition of cuts generated from the subproblems to C to be carried over to next iterations is not necessary for the correctness of Algorithm 1 as shown in the proof of Theorem 2.
The advantage of this addition is to avoid generating the same cuts at different iterations for next restricted or relaxed DDs. In fact, even invoking subproblems for relaxed DDs (in lines 18-21) can be skipped without invalidating correctness of the algorithm. The addition of this subroutine helps improve the dual bound obtained from the relaxed DDs, and thereby, trigger the pruning condition of line 17 faster.
The next example illustrates steps of Algorithm 1.
We intend to solve this problem using the DD-BD approach of Algorithm 1. The master problem Defining the dual vector π ∈ R 2 + , we obtain the dual of the above subproblem as follows.
Different iterations of solving the master problem of Example 3.
The feasibility cuts are of the formπ 1 (x 1 + x 2 ) −π 2 (0.1x 1 + 0.3) ≤ 0 whereπ is a recession ray of the dual subproblem. Similarly, the optimality cuts are of the form z +π 1 ( whereπ is a feasible solution of the dual problem. We create the exact DD D of Figure  Solving S(x) gives a feasibility cut 0.66x 1 +x 2 ≤ 1. We update the set of cuts C and refine D w.r.t to the feasibility cut. The output is the new DD depicted in Figure 4b. At this iteration, a longest path corresponding to the solution (x 1 , x 2 , z) = (0, 1, M ) gives the optimal solutionπ = (0, 6.66) and the optimal value ρ = 2 for the dual subproblem, generating the optimality cut z ≤ 0.66x 1 + 2 through solving S(x). Refining D for the last time w.r.t this optimality cut, the DD of Figure 4c is obtained.
It follows that the longest path has length w equal to 3.66 and is achieved by (x 1 , x 2 , z) = (1, 0, 2.66).
At this point, we update w * = w = 3.66 and (x * , z * ) = (x, z) = (1, 0, 2.66). Since D models an exact representation of M C (x), it follows from Remark 1 that we can skip creating relaxed DDs, which leads to terminating the algorithm with the optimal solution (x * , z * ) = (1, 0, 2.66) and the optimal value w * = 3.66.
While DDs have been used in conjunction with BD for mixed integer models in the literature as in van der Linden (2017), their layer construction has been limited to integer variables only, which has led to major computational shortcomings such as not admitting reduced forms-a critical feature for building efficient DDs in practice. Such limitations are not present in our approach, as it directly incorporates continuous variables within DD layers. We refer the reader to Appendix B for a detailed comparison.

Application to the Unit Commitment Problem
In this section, we provide an evidence of practicality for the DD-BD framework by applying it to the UCP-a well-known problem in electrid grid applications. We first present a common MIP formulation for the UCP. To apply our framework, we consider the standard BD formulation of the problem, and through a projective transformation, make it amenable to DD representation.
We close the section by comparing the computational results of the DD-BD approach with that of the standard BD formulation.

Standard Formulations
In the UCP variant we consider in this paper, we seek to find a power generation schedule that minimizes the total operational cost subject to typical UCP constraints such as minimum up/down times, generation capacity, demand satisfaction, ramping requirements, and spinning reserve. The typical MIP formulation of this problem is given in (1a) shuts down) at time j. Further, continuous variables p i j andp i j denote the production and maximum power available of unit i at time j.
The objective function (1a) of the UCP minimizes the total fixed operating costs c i f , production costs c i g , and logarithmic start-up costs q i j for all units over the planning time horizon. The fixed operating cost c i f is paid whenever unit i is up and working. The production cost c i g of unit i models the variable cost per unit of generated power. The start-up cost q i j of generator i at time j is a logarithmic function of the inactive duration of the generator to model the fact that: the colder a generator gets, the greater cost is incurred to start it again. This function is commonly linearized by discretizing the time horizon and evaluating the start-up costs K i k of unit i after it has been inactive for k consecutive time periods (Ostrowski et al. 2011, Tuffaha andGravdahl 2013). This linearization step is modeled in constraint (1b). Constraint (1c) represents the logical relation between the commitment variables x i j and variables y i j andȳ i j . constraints (1d) and (1e) model the requirement that the schedule of each generator i must satisfy a minimum down-time i ≥ 1 and a minimum up-time L i ≥ 1, i.e., if a generator i shuts down (resp. starts up) at time j, then it must be down (resp. up) for at least i (resp. L i ) time periods. In addition, each generator i has a start-up SU i , ramp-up RU i , shut-down SD i , and ramp-down RD i rates. These parameters bound the changes in production rate of generators in consecutive time periods. For instance, if unit i is down at time j, then its production is bounded above by SU i at time j + 1. Similarly, if unit i is working at time j with a generation output strictly greater than SD i , then it cannot be offline at time j + 1. These requirements are represented in constraints (1f) and (1g) as given in Tuffaha and Gravdahl (2013). It is commonly assumed that SU i ≤ RU i and SD i ≤ RD i for all i ∈ N . Further, each generator i has a minimum m i and a maximum M i production capacity that are captured in constraints (1h). Finally, constraints (1i) and (1j) guarantee that the total demand D j and the spinning reserve requirement R j are satisfied at each time period j.
The above MIP formulation is used as the core model for the two-stage stochastic UCP, where the first stage decides the commitment status of units, and the second stage determines the generation schedule to satisfy uncertain demand represented by a number of possible scenarios. Due to the L-shaped structure of this stochastic UCP, a BD method is commonly used to solve the problem; see (Zheng et al. 2013(Zheng et al. , 2014  The standard BD applied to (1a)-(1l) is composed of a master problem that contains binary variables together with the linearized cost variable, and the subproblems that are defined over continuous variables for fixed value assignments to binary variables. The master problem for this BD formulation is written as where Γ and −Γ are some valid bounds on z induced from the MIP formulation, and z represents the objective value of the following subproblem For the two-stage stochastic UCP described above, master problem (2) represents the first-stage model, and subproblems (3) are defined for each demand scenario realized at the second stage.
One way to obtain explicit numerical values for Γ bounds is to minimize/maximize the objective function of (3) over the LP relaxation of (1b)-(1l). Section 4.4 presents computational results for this stochastic problem.
When applying the BD algorithm, at each iteration, an optimal solution of the master problem (2) is found and then fed into the subproblem (3), which provides feasibility or optimality cuts to be added back to the master problem. Since the structure of the above BD formulation conforms to that used in the DD-BD approach, we adopt this model as the basis for our analysis. We show in the next section that through exploiting the DD structure, we can streamline the master problem representation, which results in a superior computational performance.

DD-BD: Master Problem Formulation
The first step in applying the DD-BD approach of Section 3 to the UCP is to construct a DD that represents the feasible region of the master problem (2). Since this set is defined over binary variables (x, y,ȳ) together with continuous variables (q, z), a typical DD would contain all five variable types in its arc layers. It turns out that it suffices to model variables x and z only in DD arc layers, while recording the value of other variables through state representations. This projection will allow for a direct application of Algorithm 1, and it will lead to a significant size reduction for DD models.
We present the DD structure for the single-unit case, i.e., n = 1. The extension to multi-unit case follows similarly by replicating the DD structure for other units, as they can be considered independently in the master problem. For this reason and to simplify notation, we drop the superscript representing unit number in variables and parameters in the sequel.
The construction of the DD representing (2) is given in Algorithm 2. In view of this algorithm, each node u ∈ U j for j ∈ {1, . . . , T + 1} contains a state value of the form (s + u , s − u ) where s + u and s − u record the number of time periods passed since the last start-up and shut-down of the unit at time j, respectively. An advantage of this state definition is its memoryless property that allows for a direct evaluation of the minimum up/down requirements as well as the logarithmic start-up cost at each time period without the need for backtracking to identify the unit status at previous periods; see the proof of Theorem 3 for a detailed analogy. The initial state value at the root node is set to (∞, ∞) to indicate that the unit is down at the start of the planning timeline and is ready to start up if decided to. The output of Algorithm 2 is a DD D = (U, A, l(.)) with T + 1 arc layers representing variables x j for j ∈ T , as well as the continuous variable z at the last layer. Each arc a ∈ A has a label l a that represents the value assignment of its corresponding variable, and a weight w a that captures the objective function rate of that assignment. The following example illustrates an exact DD for the master problem of an instance of the UCP.
Example 4. Consider a UCP instance with one generator over two time periods with the fixed cost c f , and the logarithmic start-up cost K k for k consecutive inactive time periods. Further, suppose that the minimum down-time of the generator is one and the minimum up-time is two.
The DD depicted in Figure 5 represents the feasible region of the associated master problem where −Γ and Γ are valid bounds for the continuous variable z. In this DD, the dashed and solid arcs represent arc labels with values of zero and one, respectively. In addition, the state value of each node and weight of each arc are given next to them, where K ∞ denotes the cost for the first time that the unit starts up.
The next result shows that the solution set of the equivalence class formed by the DD obtained from Algorithm 2 matches the projection of the feasible region of the master problem (2). In what follows, we refer to a solution (ẋ,ẏ,ẏ,q,ż) of (2) as non-redundant if there is no distinct point (ẋ,ẏ,ẏ,q,ż) of (2) with a strictly smaller objective value.
Theorem 3. Let D be a DD constructed by Algorithm 2 in relation to a UCP with the master problem (2). Then, a point (ẋ,ż) encodes an r-t path of lengthċ in the equivalence class formed by D, if and only if this point can be extended to a non-redundant solution (ẋ,ẏ,ẏ,q,ż) of (2) with objective valueċ.
Illustration of an exact DD in Example 4 It follows from Theorem 3 that refining the DD constructed by Algorithm 2 with respect to optimality and feasibility cuts produced from subproblems (3) results in a solution set equivalent to that of the projection of the non-redundant feasible region of (2) after addition of those cuts.
As a consequence, the proposed DD-BD approach converges to an optimal solution of the UCP problem. We will give details on the derivation of optimality and feasibility cuts in the next section.
Algorithm 2 demonstrates the construction of an exact DD representing the master problem (2), which can be used as the DD oracle required for Algorithm 1. For practical applications, however, the construction of such a DD could be computationally prohibitive. To mitigate this computational burden, as discussed in Section 3, it is necessary to use restricted and relaxed DDs with smaller width limits to improve efficiency of the algorithm. We next give details on designing a relaxed DD for (2).
At the heart of a relaxed DD lies a node-merging operation, through which multiple nodes at each layer of the DD are merged into one node, and thereby reducing the DD size. This node-merging operation must satisfy the condition that any r-t path of the exact DD remains an r-t path of the resulting DD with no greater length (for a minimization original problem.) As a byproduct, such a merging operation could create some infeasible r-t paths, hence providing a relaxation. The following definitions set the stage for building a relaxed DD for the UCP.
Algorithm 2: The construction of DD for the master problem of the UCP for n = 1 Data: parameters , L, Γ, K
create a node v ∈ U j+1 with state value (s + u + 1, s − u + 1), and an arc a ∈ A j connecting u to v with l a = 0 and w a = 0 create a node v ∈ U j+1 with state value (1, s − u + 1), and an arc a ∈ A j connecting u to v with l a = 1 and w a = c f + K s − u 7 else 8 create a node v ∈ U j+1 with state value (s + u + 1, s − u + 1), and an arc a ∈ A j connecting u to v with l a = 1 and w a = c f 9 if s + u ≥ L then 10 create a node v ∈ U j+1 with state value (s + u + 1, 1), and an arc a ∈ A j connecting u to v with l a = 0 and w a = 0 11 forall u ∈ U T +1 do 12 create two arcs a 1 , a 2 ∈ A T +1 connecting u to the terminal node with l a 1 = w a 1 = Γ, Definition 2. Consider the following modification of Algorithm 2: (i) each node u ∈ U j for j ∈ {1, . . . , T + 1} is assigned with state values (s + u , s − u , s = u ) with an additional component s = u whose value is initialized and updated similarly to that of s − u throughout the iterations of the algorithm; (ii) parameter K s − u is replaced with K s = u when calculating arc weights in line 6 of the algorithm.
We next show that the modified algorithm of Definition 2 combined with the node-merging operation of Definition 3 yield a relaxed DD for the master problem (2).
Proposition 3. LetD be the DD obtained from Algorithm 2 in relation to the master problem (2). LetD be a DD constructed from the modified algorithm of Definition 2 in combination with the node-merging operation of Definition 3. Then,D provides a relaxation ofD.

DD-BD: Subproblem Formulation
Through application of Algorithm 1 to the UCP, subproblem (3) is solved for fixed values of master variables to obtain feasibility and optimality cuts, with respect to which the restricted and relaxed DDs representing the master problem are refined. These cuts, however, are generated in the space of variable (x, y,ȳ, z). As a result, they cannot be directly used to refine the DDs representing the master problem as they only contain variables (x, z); see Section 4.2. To resolve this discrepancy, we next show that constraints (1f) and (1g) in the description of the UCP can be replaced with two constraints that do not contain variables (y,ȳ). Such a substitution in the subproblem (3) leads to feasibility and optimality cuts that are in the space of (x, z), which can be directly used to refine DDs of the master problem.
Proposition 4. Assume that SU i ≤ RU i and SD i ≤ RD i for all generators i ∈ N . Let E be the feasible region of the UCP described by (1b)-(1l). Define G to be the feasible region of the UCP where constraints (1f) and (1g) are respectively replaced with Then, E = G.
Substituting (1f) and (1g) with (4a) and (4b) in the description of the subproblem (3), we obtain the following dual problem When applying Algorithm 1, we solve the above dual problem for each given assignmentx. If the dual subproblem is unbounded (i.e., the subproblem is infeasible,) we obtain a ray (ψ,β,φ,π,γ,δ) and add the following feasibility cut to the master problem.

Computational Experiments
In this section, we assess the performance of the proposed DD-BD algorithm in solving the stochastic UCP instances for a short-term (day-ahead) power generation model. For our experiments, we consider benchmark instances from the problem class OR-LIB/UC in UnitCommitment.jl (Xavier et al. 2020). Since these instances lack some of the modeling elements in the UCP variant we consider in this paper, we set their associated parameter values as follows. We set the production cost c i g for unit i ∈ N as the average of production costs for that unit. We create the fixed operating cost c i f for each unit i ∈ N randomly from the interval [400, 1000]. We set RU i = SU i and RD i = SD i for each unit i ∈ N . We capture demand uncertainty by a set of scenarios Ξ that are generated randomly from the interval [0.75M ,M ] whereM denotes the total capacity over all generators.
All experiments were conducted on a machine running Windows 10, x64 operating system with Intel ® Core i7 processor (2.60 GHz) and 32 GB RAM. The code is written in Microsoft Visual Studio 2015 in C++. For the optimization parts, we use the Gurobi solver version 9.0.0.

Comparison with Decomposition Methods We first present computational results
that compare the outcome of the DD-BD method, as a decomposition technique, with other standard decomposition methods. In particular, we consider two BD approaches referred to as G-BD 1 and G-BD 2.
G-BD 1 employs the classical BD in a textbook fashion, where the master problem (2) is solved by Gurobi at its default settings, and the optimality/feasibility cuts are produced by solving the subproblems (3) for each demand scenario.
Since the traditional implementation of BD, as used in G-BD 1, is known to suffer from a slow convergence, we employ modern modeling and algorithmic practices in the literature that boost the performance of the G-BD. We refer to this improved approach as G-BD 2. The details of these improvements are as follows. (i) To accelerate the convergence rate associated with the feasibility cuts, we use the L 1 normalization technique proposed by by Fischetti et al. (2010) and Bonami et al. (2020). Consider the MIP H defined in Section 3. Assume that the subproblem S(x) is infeasible for a fixed solutionx obtained from the master problem. It follows that w * = min w,y {w | By + w1 ≤ c − Ax} = π * (c − Ax) > 0, where π * is the optimal dual vector of the normalized subproblem, 1 is a unit column vector, and the equality holds because of the strong duality property. Therefore, the inequality π * (c − Ax) ≤ 0 can be added to the master problem as a feasibility cut that excludesx.
(ii) During presolve, we generate heuristic cuts to tighten the initial LP relaxation of the master problem as follows. At each time period, we identify the scenario with the highest demand and add valid inequalities in the space of commitment variables to ensure the feasibility of constraints (1h) and (1j). (iii) We derive the so-called combinatorial cuts that exclude the current infeasible solution of the master problem by forcing a change of value in at least one commitment variable; see Codato and Fischetti (2006) and Rodríguez et al. (2021) for details about these cuts. (iv) Finally, we strengthen the logarithmic start-up cost constraints (1b) by reducing the coefficient of variables x i j−h as proposed in Silbernagl et al. (2015) and Knueven et al. (2020a).
For the DD-BD method, we set the width limit for the restricted/relaxed DDs to two and use the node-merging operation given in Definition 3 for the construction of relaxed DDs. When necessary during the refinement procedures, we allow the increase in the DD width to accommodate for the separation of current non-optimal solutions.
For these experiments, we fix the number of demand scenarios at 100. We consider four size categories with n ∈ {20, 50, 100, 150} for the number of generators, and solve five 5 instances for each problem size as given in the OR-LIB/UC. Table 1 contains the result of solving these instances with G-BD 1, G-BD 2 and DD-BD. For these runs, we fix the time limit at 7200 seconds. The first column of Table 1 indicates the number of generators. The second column shows the instance number. The columns under "G-BD 1" report the solution time (in seconds), the total number of generated feasibility and optimality cuts during the solution process. The subsequent columns contain the same quantities for G-BD 2 and DD-BD. As observed in Table 1, while the boosting methods implemented in G-BD 2 have substantially improved the solution time over G-BD 1, the DD-BD approach uniformly outperforms both G-BD alternatives across all size categories, rendering it as the most effective decomposition approach.

Comparison with Full Formulations
In this section, we compare the performance of the decomposition methods DD-BD and G-BD 2 with that of the full formulation (1a)-(1l) when solved by Gurobi at its default settings. Since it is well-known that decomposition methods are Table 1 Solution times (in seconds) of G-BD 1, G-BD 2, and DD-BD for instances of the stochastic UCP with |Ξ| = 100.  Tables 2 and 3 for 20 and 50 generators, respectively. The solution time for larger number of generators exceeds the time limit of 7200 seconds for all methods. In these tables, the first column shows the number of demand scenarios, and the second column contains the instance number. The next three columns under "G-BD 2" report the solution time, the number of feasibility and optimality cuts, respectively. Columns 6-8 contain similar results for the DD-BD method. The last column shows the solution time for the full formulation.
As evident in Tables 2 and 3, similarly to Table 1 We conclude this section by noting that the presented computational results are obtained through a basic implementation of the DD-BD without incorporating any computational enhancement methods for DDs such as variables ordering, DD reduction, dynamic width adjustment, primal heuristics, etc. As common in the literature, such basic DD implementations are often compared against modern solvers with the presolve/cuts/heuristics settings turned off to provide a fair comparison for the core solution methodology. Despite being at such disadvantage, the above results indicate that the proposed DD-BD method can still outperform the G-BD and the full UCP formulation solved with the solver's default settings. While this shows the potential of the DD-BD method in energy applications, its performance can be further improved using modern computational practices in the DD literature, which can be considered as a direction for future work.

Conclusion
In this paper, we propose a geometric decomposition that restructures a mixed integer set through a collection of hyper-rectangle formations. We show that this decomposition approach is the key to model optimization problems with both integer and continuous variables through decision diagrams. In this regard, we introduce a method, referred to as DD-BD, that extends the applicability domain of decision diagrams to MIPs. We evaluate the performance of DD-BD by applying it to the unit commitment problem in the electric grid market.
Proof of Proposition 1. Consider any pointx N \I ∈ proj x N \I (P), and define the indicator We can use δ(x|x N \I ) as the objective function in the relation max{f (x) | x ∈ P} = max{f (x) | x ∈ Q} as it is convex in x I . This follows from the fact that δ(x|x N \I ) = 0 when x N \I =x N \I , and it is −∞ otherwise, which is improper convex; see Rockafeller (1970). We conclude thatx N \I ∈ proj x N \I (Q). Using a similar argument for the reverse direction, we conclude that proj x N \I (P) = proj x N \I (Q). This also shows that proj x N \I (P) is finite as Q is finite by assumption. For any For any convex function f I (x I ) defined in the space of variables x I , and any pointx N \I ∈ proj x N \I (P), we have that where (8a) and (8d) follow from the definition of the indicator function above, and (8b) follows from the assumption. Since f I (x I ) is convex, the optimal value of (8a) is attained at an extreme pointx of P(x N \I ). For each such extreme point, there are infinitely many convex functionsf I (x I ) withx as their unique maximizer over P(x N \I ). For every one of these functions, according to the chain equalities (8a)-(8d), Q has a point that matches the optimal value. Since Q is finite, it must be thatx ∈ Q(x N \I ). This shows that the set of extreme points of P(x N \I ) is finite. Similarly, it can be shown that for any extreme pointẋ of conv(Q(x N \I )), we have thatẋ ∈ P(x N \I ). As a result, conv(P(x N \I )) = conv(Q(x N \I )) for everyx N \I ∈ proj x N \I (P). Now, construct sets P j I in Theorem 1 as conv(P(x N \I )) for everyx N \I ∈ proj x N \I (P), which yields a finite collection.
Condition (i) is satisfied as each set P j I is restricted at {x|x N \I =x N \I }. Condition (ii) holds since for any extreme point of P j I , there is a point of P by construction, and since any pointx ∈ P satisfiesx ∈ P(x N \I ) ⊆ conv(P(x N \I )). For condition (iii), rectangles R k j can be considered as the set of extreme points of P j I , which has been shown to be finite.
Proof of Lemma 1. The result is obtained as a special case of Theorem 1, where sets P j I and R k j coincide, i.e., P j I = R 1 j for all j ∈ J. The conditions and the result follow immediately.
Proof of Proposition 2. For a given DD D, define a node-sequence as an ordered set of connected nodes from the root to the terminal, i.e., u = (u 1 , u 2 , · · · , u n+1 ) where u i ∈ U i for i ∈ N ∪ {n + 1}. Let U be the collection of all node-sequences of D. For u ∈ U , define Viewing Sol(D) as a compact set P in Lemma 1, it is straightforward to verify that sets S I (u) satisfy the conditions for hyper-rectangles R j I . It also follows from the definition of (virtual) DDs that Sol(D) = u∈U S I (u) and Sol(D) = u∈U X (S I (u)). The result follows from Lemma 1.
Proof of Corollary 1. For the direct implication, assume that a given compact set P admits a rectangular decomposition w.r.t. I through sets P j I for j ∈ J. Then, the DD that encodes the finite collection of points j∈J X (P j I ) provides the desired DD representation because of Theorem 1. For the reverse implication, assume that a given compact set P is DD-representable w.r.t. I through a DD D. Since Sol(D) is finite, it follows from Proposition 1 that P admits a rectangular decomposition w.r.t. I.
Proof of Corollary 2. Since Q is bounded, there are finitely many pointsx ∈ proj x (P). For each such point, it follows from the definition of projection and the boundedness of Q that there exists an interval [lx, ux] such that (x;ȳ) ∈ P for everyȳ ∈ [lx, ux]. As a result, we can write that P = x∈proj x (P) Rx y where Rx y = {(x; y) ∈ R n+1 | x =x, y ∈ [lx, ux]}. Hyper-rectangles Rx y satisfy the conditions of Lemma 1, and hence admits a rectangular decomposition w.r.t. the index of variable y which is n + 1. The result follows from Corollary 1.
Proof of Theorem 2. First, we show that the Algorithm 1 terminates after a finite number of iterations. To this end, we argue that each loop in the algorithm is repeated for a finite number of iterations. The outer while loop is executed for each member of the partial assignment setX .
These partial assignments are generated based on longest r-u paths associated with nodes u in the exact cut set of relaxed DDs D, which contain a finite number of paths by definition; see line 22 of the algorithm. Each resulting partial assignmentx is different due to the structure of DDs that do not admit multiple paths with similar encoding values. As a result, an upper bound for the number of partial assignments that can be included inX is all possible partial assignments of the solution set of M for x variables. Because the discrete set P in the description of M is bounded by assumption, we conclude thatX contains a finite number of elements. For the repeat loop in lines 6-10, the goal is to obtain the optimal value of the solution set represented by D subject to the constraints of the subproblem S(.). It follows from the property of BD method that the resulting optimality and feasibility cuts correspond to the extreme points and extreme rays of the dual of the subproblem, which is independent of the choice of the fixed point x. Since the feasible region of this dual problem is a polyhedron, the set of its extreme points and rays are finite, which yield a finite set of optimality and feasibility cuts that can be added through the execution of this loop.
It is also known in the BD literature (Conforti et al. (2014a)) that the same optimality/feasibility cuts cannot be generated repeatedly through different iterations, since the added cuts remain in the description of the master DD due to refinement. A similar argument can be used to show that the number of iterations of the repeat loop in lines 18-21 is finite, thereby yielding the result.
Next, we show that the outputs (x * , z * ) and w * of Algorithm 1 give an optimal solution and optimal value of H, respectively. To this end, we first prove that (x * , z * ) is a feasible solution to H and w * is a lower bound for its optimal value. This solution is updated at line 12 of Algorithm 1 as the point encoding a longest r-t path of D, after being refined with respect to optimality and feasibility cuts generated from subproblems S(.). We write that ( where the second inclusion holds because D is a restricted DD associated with M C (x) for some C andx, and the third inclusion follows from the fact that a feasible solution to M C (x) is feasible to M by definition. Further, the BD structure implies that the point encoding a longest path obtained at the termination of the loop in line 6-10 satisfies the constraints of the subproblem S(.).
As a result, (x * , z * ) satisfies the constraints of both the master and the subproblem, hence being feasible to H. Using a similar argument, we obtain that w * is a lower bound for the optimal value of M subject to constraints in S(.), since w * is the length of the longest path in the restricted DD representing a restriction of M, after valid optimality and feasibility cuts are added based on the subproblem. Now, we show that (x * , z * ) is indeed an optimal solution to H. Assume by contradiction that there exists an optimal solution (x,z) of H with optimal valuew such that w > w * . There are three cases. (i) Assume thatx is added as a partial assignment to setX at some iteration of the algorithm. For the while loop wherex is selected at line 3, all x variables are fixed for the restricted DD D. Therefore,x is used as the input for the subproblem, i.e., we solve S(x), which yields the optimal valuez since (x,z) is an optimal solution of H. Since the weights on the arcs of the restricted DD D are set as the coefficients of variables in the linear objective function of H, the length of the r-t path of D encoded by (x,z) isw. If follows from the contradiction assumption thatw > w * in line 11 of Algorithm 1, which leads to updating the optimal solution and optimal value to (x,z) andw. As a result, (x * , z * ) and w * cannot be returned by the algorithm as the output, a contradiction. (ii) Assume thatx is not added as a partial assignment to setX because the algorithm terminates before such a partial assignment is reached in line 22. This implies that there must be a partial assignmentx ∈X withx i =x i for i = 1, . . . , j − 1 for some j ∈ N such that the relaxed DD D associated with M C (x) for some C is pruned without reaching line 22. The only possibility for this event is that the lengthw of the longest r-t path of D must satisfyw ≤ w * , violating the condition in line 17. However, because of the facts that D is a relaxed DD associated with M C (x), and that (x,z) must be a feasible solution to M C (x), we conclude thatw ≥w. Combining this inequality with that of the line above, we obtainw ≤ w * , a contradiction to the initial assumption on the optimal value of the problem. (iii) Assume thatx is not added as a partial assignment to setX because the longest r-u path chosen in line 22 deviates fromx for some node u of the the path associated withx. This implies that there must be a partial assignmentx ∈X withx i =x i for i = 1, . . . , j − 1 for some j ∈ N such that the relaxed DD D associated with M C (x) for some C has an exact cut set that contains node u in some layer k ≥ j. Letẋ be the solution encoding the longest r-u path of D that is chosen in line 22 in place ofx. It follows from definition of exact cut set that an extension ofẋ with componentṡ x i =x i for i = k, . . . , n andż =z is a feasible solution to H. However, the above assumption implies that the lengthẇ of the path encoded by (ẋ,ż) is no less thanw. Ifẇ >w, then we have reached a contradiction to the assumption thatw is the optimal value of H. Ifẇ =w, we may repeat the contradiction arguments for the new point (ẋ,ż) until we exhaust all possible replacements for the assumed optimal solution. The last such solution must either fall in case (i) or case (ii) above, reaching contradiction.
Proof of Theorem 3. First, we give a useful observation from Algorithm 2. It follows from this algorithm that state values (s + u , s − u ) at each node u ∈ U j records the number of time periods passed since the last start-up and shut-down of the unit at time j, as these values reset to 1 whenever there is change in the unit status over two consecutive periods, and are incremented by one if the unit status remains the same. These state values imply the status of the unit at the beginning of each period, i.e., the unit is down when s + u ≥ s − u , and it is up otherwise. For the direct implication, assume that (ẋ,ż) encodes an r-t path of lengthċ in the equivalence class formed by D. We construct an extended solution (ẋ,ẏ,ẏ,q,ż) of (2) with objective valueċ as follows. First, we show thatẋ j ∈ {0, 1} for all j ∈ T andż ∈ [−Γ, Γ], thereby satisfying domain constraints in (2). The definition of equivalence class implies that, for each j ∈ T , there exists a node pair (u, v) ∈ A j × A j+1 such thatẋ j ∈ [l min (u,v) , l max (u,v) ], i.e., the variable value belongs to the interval defined by the min and max label values of the arcs connecting nodes u and v; see Section 2.3.
It follows from the construction of D in Algorithm 2 that two consecutive nodes in layers 1, . . . , T cannot be connected with multiple arcs with different label values, since different arc labels lead to different state values for the head node. As a result, we must have l min (u,v) = l max (u,v) ∈ {0, 1}. The argument forż ∈ [−Γ, Γ] follows directly from the construction and label values assigned to the arcs at the last layer of D. To construct the extended solution, for each j ∈ T , defineẏ j = 1 iḟ x j−1 = 0 andẋ j = 1, andẏ j = 0 otherwise. Similarly, defineẏ j = 1 ifẋ j−1 = 1 andẋ j = 0, andẏ j = 0 otherwise. These definitions guarantee the satisfaction of (1c). Further, let u ∈ U j be the node at layer j of the path encoding (ẋ,ż), and defineq j = K s − u ifẋ j = 1 andẋ j−1 = 0, andq j = 0 otherwise. These value assignments satisfy (1b) because this constraint implies that when a unit changes status from down to up, i.e., x j = 1 and x j−1 = 0, then q j ≥ max k=1,...,s − u K k , as s − u represents the number of periods that the unit has been down consecutively before going up. Since the start-up cost function is assumed to be logarithmic, we have that max k=1,...,s − u K k = K s − u . Further, sincė q j takes the maximum value at equality, it yields a non-redundant solution. It remains to show that the constructed point satisfies constraints (1d) and (1e). Assume by contradiction that there exists j * ∈ T for which constraint (1d) is violated by (ẋ,ẏ,ẏ,q,ż). Note that (1d) models two requirements at time j * : (i) if the unit is down, then it could not have started up in the last L time periods; and (ii) if the unit is up, then it could not have started up more than once in the last L time periods. For the contradiction, assume first that condition (i) is violated, i.e.,ẋ j * = 0 and there existsj ∈ {j * − L + 1, . . . , j * } such thatẏj = 1. It follows from construction thatẋj = 1 anḋ xj −1 = 0. Let (u, v) ∈ Uj × Uj +1 be the nodes at layersj andj + 1 of the r-t path encoding (ẋ,ż).
We have s + u ≥ s − u because of the earlier argument on the state values. Since the arc connecting u to v on the r-t path is a 1-arc (ẋj = 1), we have s + v = 1; see line 6 of Algorithm 2. It follows from the algorithm steps that the only way for the unit to be able to shut down afterj is to satisfy the condition of line 9, i.e., s + h ≥ L for some node h in layersj + L, . . . , T . Sincej ≥ j * − L + 1 by definition, we obtain that j * <j + L, and henceẋ j * cannot be equal to zero, a contradiction.
Next, assume that condition (ii) above is violated for the contradiction assumption, i.e., the unit starts up more than once in periods j * − L + 1, . . . , j * . It is easy to verify that there exists a layer j ∈ {j * − L + 1, . . . , j * } for which condition (i) is violated. Therefore, an argument similar to that of condition (i) yields the desired contradiction. The contradiction for violating constraint (1e) is obtained similarly due to the symmetry in the problem structure. We conclude that (ẋ,ẏ,ẏ,q,ż) is feasible to (2). We next show that the length of the r-t path encoding (ẋ,ż) matches the objective value of (ẋ,ẏ,ẏ,q,ż). The proof follows from considering the contribution of three terms in the objective function of (2). First, the contribution of each variable assignmentẋ j = 1 to the objective function is c f , which is captured in the weight of the associated 1-arcs of the r-t path through Algorithm 2 in lines 6 and 8. Second, the contribution of the start-up status of the unit to the objective function isq j = K 1 s − u where u is the node at layer j of the path encoding (ẋ,ż), which is captured in line 6 of Algorithm 2. Third, the contribution ofż in the objective function is directly considered in the arc weight in the last layer of the equivalence class of D.
For the reverse implication, assume that (ẋ,ẏ,ẏ,q,ż) is a non-redundant solution of (2) with objective valueċ. We show that the projected point (ẋ,ż) encodes an r-t path of lengthċ in the equivalence class formed by D. Assume by contradiction that no such path exists in the equivalence class of D. Therefore, there exists an r-t path P of D with associated point (x,z) and a layer number j * ∈ T \ {T } such thatx j =ẋ j for j ∈ {1, . . . , j * } and the node u at layer j * + 1 of P does not have an outgoing arc with label valueẋ j * +1 . There are four cases. For the first case, assume thatẋ j * = 0 andẋ j * +1 = 0. It follows from the state value definitions that s + u ≥ s − u asẋ j * =x j * = 0. Line 4 of Algorithm 2 implies that u has an outgoing arc with label valueẋ j * +1 = 0, a contradiction. For the second case, assume thatẋ j * = 1 andẋ j * +1 = 0. Because of constraints (1c) and (1d), we must haveẋ j =x j = 1 for j = j * − L + 1, . . . , j * . Let v be the node at layer j * − L + 1 of P. Sincex j = 1 for j = j * − L + 1, . . . , j * , it follows from line 8 of Algorithm 2 that s + u = s + v + L − 1 ≥ L as s + v ≥ 1. As a result, the condition of line 9 of Algorithm 2 is satisfied at node u, which leads to creating a 0-arc as an outgoing arc of u, a contradiction. For the third case, whereẋ j * = 1 andẋ j * +1 = 1, the contradiction is obtained similarly to the first case due to the symmetry of the problem. For the fourth case, whereẋ j * = 0 andẋ j * +1 = 1, the contradiction is achieved similarly to the second case due to the symmetry of the problem. Finally, sinceż ∈ [−Γ, Γ], it must belong to the equivalence class of D as −Γ and Γ are used as label values of the pair of arcs connecting every node of layer T + 1 to the terminal node in D. We conclude that (ẋ,ż) encodes an r-t path of lengthċ in the equivalence class formed by D. To show that the objective valueċ is equal to the length of the r-t path encoding (ẋ,ż), we note thatq 1 where h is the node at layer j of the r-t path, since otherwise (ẋ,ẏ,ẏ,q,ż) would be a redundant solution of (2). The rest of the proof follows from an argument similar to that of the direct implication case given above.
Proof of Proposition 3. We need to show that for each pathP ofD with encoding point (x,z), there exists a pathP ofD with encoding point (x,z) such thatx =x,z =z, and w(P) ≥ w(P), where w(.) represent the length of the path. DefineP k to be the path composed of the first k arcs ofP for k ∈ T . We prove the result by induction on the size of k, i.e., we show that for eachP k with encoding pointx k , there exists a pathP k ofD with encoding pointx k such thatx k =x k and w(P k ) ≥ w(P k ). Note that it suffices to show the above result for k up to T , since each node at node layer T + 1 is connected to the terminal node by two arcs with labels −Γ and Γ in bothD andD, and therefore a matching arc can always be found for the desired path. For the base of induction, i.e., k = 1, it follows from Algorithm 2 thatD has two arcs with labels 0 and 1 going out of the root node. These two arcs remain inD with the same weights according to Definition 2. For the induction hypothesis, assume that, for some k ≥ 2, there exists a pathP k−1 ofD with encoding pointx k−1 such thatx k−1 =x k−1 and w(P k−1 ) ≥ w(P k−1 ). We show that there exists a pathP k of D with encoding pointx k such thatx k =x k and w(P k ) ≥ w(P k ). There are four cases. (i) Assume thatx k k−1 = 0 andx k k = 0. Note thatx k k−1 =x k−1 k−1 = 0 as they are defined over the same pathP. It follows from the induction hypothesis thatx k−1 k−1 =x k−1 k−1 = 0. Let u be the node at node layer k of the pathP k−1 . We consider two cases for this node. For the first case, assume that u is not a merged node. As a result, it has been created directly from the modified algorithm of Definition 2, which implies that s + u ≥ s − u . Following line 4 of Algorithm 2, there is a 0-arc going out of u inD. Adding this arc toP k−1 yields the desired pathP k , as the contribution to the objective function is 0 for this variable assignment in both pathsP k andP k . For the second case, assume that u is a merged node. Let {v 1 , . . . , v r } be the set of nodes that have been merged into u using the merging operation of Definition 3. Since u has an incoming arc with label 0, it means that one of the nodes v i in the above set must have an incoming arc with label 0. The modified algorithm implies that . . , r}. It follows from the merging equations that s + Therefore, an argument similar to that of the first case yields the desired result. (ii) Assume thatx k k−1 = 0 andx k k = 1. Let h be the node at node layer k ofP k . It follows from line 5 of Algorithm 2 that s − h ≥ and that the arc weight at layer k is c f + K s − Let u and v be the nodes at node layers k − s − h + 1 and k ofP k−1 , respectively. We consider two cases for node u. For the first case, assume that u is not a merged node. As a result, it has been created directly from the modified algorithm of Definition 2, which implies that s − u = s = u = 1. For the second case, assume that u is a merged node. Let {u 1 , . . . , u r } be the set of nodes that have been merged into u using the merging operation of Definition 3. Since u has an incoming arc with label 0 right after an arc with label 1, it implies that one of the nodes u i in the above set must have s − u i = s = u i = 1. It follows from the modified algorithm that s − u ≥ 1 and s = u ≤ 1, which holds for the first case as well. Using an argument similar to that given above for the nodes where the first inequalities follow from the state definitions through merging operation applied on successive node layers, and the second inequalities hold because of the relations above on state values at node u. As a result, the condition in line 5 of Algorithm 2 is satisfied at node v as meaning that it has an outgoing arc with label 1 and weight c f + K s = v ≤ c f + K s = h . Adding this arc toP k−1 yields the desired pathP k . (iii) Assume thatx k k−1 = 1 andx k k = 1. The result is obtained by using an argument similar to that of case (i) due to symmetry. (iv) Assume that x k k−1 = 1 andx k k = 0. An argument similar to that of case (ii) proves the result due to symmetry.
Proof of Proposition 4. We prove the equivalence by showing that for each feasible variable assignment, both sets of constraints impose the same restrictions. There are four cases for each i ∈ N and j ∈ T . (i) Assume that x i j−1 = x i j = 0. It follows from constraints (1c) and (1d) that y i j =ȳ i j = 0. In this case, the right-hand-side (RHS) of both inequalities (1f) and (4a), as well as (1g) and (4b) reduce to zero. (ii) Assume that x i j−1 = 0 and x i j = 1. It follows from constraint (1c) that y i j = 1 andȳ i j = 0. In this case, the RHS of both constraints (1f) and (4a) reduce to SU i . Further, inequality (1h) implies that p i j−1 = 0 as x i j−1 = 0. As a result, both constraints (1g) and (4b) become redundant in the description of E and G, since their LHS is a non-positive quantity −p i j ≤ 0 and their RHS are non-negative values as RD i ≥ 0 and RD i − SD i ≥ 0 by assumption. (iii) Assume that x i j−1 = 1 and x i j = 0. It follows from constraint (1c) that y i j = 0 andȳ i j = 1. In this case, the RHS of both constraints (1g) and (4b) reduce to SD i . Further, inequality (1h) implies that p i j = 0 as x i j = 0. As a result, both constraints (1f) and (4a) become redundant in the description of E and G, since their LHS is a non-positive quantity −p i j−1 ≤ 0 and their RHS are non-negative values as RU i ≥ 0 and RU i − SU i ≥ 0 by assumption. (iv) Assume that x i j−1 = x i j = 1. It follows from constraints (1c) and (1e) that y i j =ȳ i j = 0. In this case, the RHS of both (1f) and (4a) reduce to RU i , and the RHS of both (1g) and (4b) reduce to RD i .

Appendix B: Comparison of DD-BD Method with the Literature
In the literature, DDs have been used in conjunction with BD framework in different contexts. In this section, we give a detailed comparison between DD-BD approach of Section 3 and those used in the literature. In particular, we show that our framework generalizes other existing approaches, while addressing their shortcomings.
van der Linden (2017) proposes a BD framework for linear mixed integer programs that uses a similar approach to ours in defining the master problem over integer variables, and solving the subproblems for continuous variables. A DD is constructed to represent the solution set of the master problem, and is successively refined through cuts obtained from the subproblems. There is, however, a fundamental difference between our approach and that of van der Linden (2017) as they use separation for the feasibility cuts only since their proposed DD cannot contain any continuous variable. To incorporate the optimality cuts, the authors propose a so-called "costtuple" approach where, for each optimality cut, a cost parameter is calculated at each node of the DD to record the contribution of the variable assignments in relation to that cut. Such an approach requires the underlying DD to be non-reduced, limiting its practicality. In contrast, our proposed DD-BD approach provides a unified framework that treats both feasibility and optimality cuts similarly through separation as the continuous variable is directly incorporated in the DD layers.
To elaborate, we first give a brief review of the cost-tuple approach used in van der Linden (2017) to handle optimality cuts.
The cost-tuple approach is given in Algorithm 3 where, without loss of generality, the master problem of a MIP is defined as max{z | x ∈ P} with a bounded P ⊆ Z n . Consider D = (U, A, l(·)) to be the DD representing P, with a property that each node (except the terminal) has a unique incoming arc, and the terminal receives a unique arc from each node of the previous layer. Let inequalities of the form z ≤ α j x + α j 0 , for j ∈ J, represent the set of added optimality cuts at the current iteration of the BD algorithm. For each node u ∈ U, parameter r j (u) records the reward (cost, for a minimization variant) contribution of variable assignments at that node to the right-hand-side value of cut z ≤ α j x + α j 0 for j ∈ J. In this algorithm, the accumulated reward r(u) captures the contribution amount at the intersection of optimality cuts for each unique value assignment of variables, while the final reward r * computes the optimal value of max{z | x ∈ P, z ≤ α j x + α j 0 , ∀j ∈ J}. The main shortcoming of Algorithm 3 is due to the assumption that the DD nodes have a unique incoming arcs. This requirement leads to an exponential growth of the DD size to distinguish between each r-t path. As a consequence, it is mentioned in van der Linden (2017) that the costtuple algorithm is not applicable to reduced DDs. Reduction in DDs is referred to the operation of merging nodes that share a similar subtree; hence it is regarded as a crucial component in building efficient DDs for practical applications and keeping them within a width limit. The DD-BD approach proposed in the present paper, on the other hand, can be applied to reduced DDs, and the separation of merged nodes, if needed, will be performed naturally through the separation operation. We illustrate this difference in the following example.
For the DD-BD approach of Algorithm 1, we represent the master problem as M = max{z | x ∈ P, z ∈ [−M, M ]} where valid bounds are imposed on z variable. A restricted DD D 2 can be created for M using the equivalence class argument of Section 2.3 as given in Figure 8b. Next, we may simply refine the DD through a sequential separation with respect to optimality cuts (9a) and (9b), as shown in Figures 9a and 9b. It is easy to verify that a longest path in the resulting DD gives either point (0, 0) or (1, 0) with the optimal value z * = 0. Since these solutions give the optimal value of the problem, Algorithm 1 is terminated due to bound pruning conditions in line 17.
(a) D 1 without continuous variable (b) D 2 with continuous variable Separation of optimality cuts for the master problem of Example 5 We conclude this section by showing that the DD-BD approach of Algorithm 1 is a generalization of the cost-tuple approach of Algorithm 3, since the latter can be viewed as a special case of the former when the underlying DD is not reduced.
Proposition 5. Consider the master problem max{z | x ∈ P} of a MIP with a bounded P ⊆ Z n , and let D = (U, A, l(·)) be the DD representing P, with a property that each node (except the terminal) has a unique incoming arc, and the terminal receives a unique arc from each node of the previous layer. Let inequalities of the form z ≤ α j x + α j 0 , for j ∈ J, represent the set of optimality cuts to be added at the current iteration of the BD algorithm. DefineD = (Ū,Ā,l(·)) to be the DD constructed from D by adding a layer representing z variable as the last arc layer. Then, the reward value obtained from Algorithm 3 applied to D is equal to the optimal value obtained from Algorithm 1 when applied toD.
Proof. We first describe the relation between D andD. All nodes and arcs from layer one to layer n − 1 are similar in both DDs, i.e., A i =Ā i for i ∈ N \ {n} and U i =Ū i for i ∈ N . The difference is thatD has an extra layer that represents z variable. This layer is created as follows.
For each node u ∈Ū n , we create a node v ∈Ū n+1 , and connect them with an arc a ∈Ā n with the same label as the unique arc connecting the replica of node u in U n to the terminal node of D. Then, for each node v ∈Ū n+1 , we create two arcs with labels −M and M in layerĀ n+1 that connect v to the terminal node inD. These label values represent some valid lower and upper bounds for variable z. According to Algorithm 1,D is refined with respect to the optimality cuts, and then the longest path is found over the refined DD weighted by the objective function rates.
To perform the refinement, the solution set satisfying each optimality cut is modeled by a DD called threshold diagram; see  for an exposure to threshold diagrams. The last layer of such DDs corresponds to variable z and contains two arcs with labels representing a lower and upper bound for z. Since for each node v ∈Ū n+1 , there exists a unique arc-specified path P = (a 1 , . . . , a n ) ∈Ā 1 × . . . ×Ā n from the root node to v, it follows from the assumption that the threshold diagrams have the same structure asD, with the difference that the upper bound label on arcs connecting v to the terminal is computed as M j (v) = n i=1 α j i l a i + α j 0 for each optimality cut j ∈ J. As a result, refiningD with respect to the optimality cuts reduces to intersecting label value intervals associated with arcs at the last layer, which are of the form [−M, M j (v)] for each node v ∈Ū n+1 . This intersection yields the lower bound label −M , and the upper bound label M (v) = min j∈J M j (v) for pair of arcs connected to v in the last layer of the refined DD. Note that M (v) is equal to r(v) as computed in Algorithm 3. Finally, to find the longest path on the refinedD, all arcs in layers 1-n are assigned zero weight and arcs in the last layer are assigned weight 1 deduced from the objective function of the master problem. Therefore, the length of the longest r-t path inD is computed as w * = max v∈Ū n+1 M (v), which is equal to r * obtained from Algorithm 3. SinceD represents an exact DD for the master problem and has been refined with respect to optimality cuts, w * is returned as the output of Algorithm 1.