View source for Nonconvex Generalized disjunctive programming (GDP)

Author: Kaiwen Li (ChE345 Spring 2015) <br/>
Steward: Prof.Fengqi You, Dajun Yue
<br/>
=Introduction=
General disjunctive programming, GDP, is an alternative approach to represent the formulation of traditional Mixed-Integer Nonlinear Programming, solving discrete/continuous optimization problems. By using algebraic constraints, disjunctions and logic propositions, Boolean and continuous variables are involved in the GDP formulation. The formulation process of GDP problem are more intuitive, and the underlying logic structure of the problem can be kept, so that solution can be found more efficiently. GDP has been successfully applied to Process Design and Planning and Scheduling areas.<br/>
<br/>
However, funtions in GDP problem sometimes could be nonconvex. Due to nonconvexites, conventional MINLP algorithms are often trapped in suboptimal solutions. Thus, solutions to nonconvex GDP has been receiving increasing attention.<br/>

=Algorithm= 
In the GDP formulation, some functions might be nonconvex, which will lead to a nonconvex GDP problem. Traditional algorithms for MINLP such as Generalized Benders Decomposition (GBD), or Outer Approximation (OA) will fail to provide a global optimum because the solution of the NLP subproblem can be only a local optimum, while the cuts in master problem may not be valid either. Therefore, in order to get a global optimum, we need to introduce special algorithm for nonconvex GBD problems.<br/>
Most of the methods rely on spatial branch and bound method.

==Motivation==
In order to find the best design and operation of a system in many problems in engineering, a set of algebraic expressions are used to reduce the problem with continuous and discrete variables. Normally theses problems lead to Mixed-Integer Nonlinear Programming (MINLP). However, in order to represent the behavior of physical, chemical, biological, financial or social systems accurately, nonlinear expressions are often used and leads to a MINLP with a nonconvex solution space. Instead of reaching a global optimal, this may give rise to local solutions that are suboptimal. Meanwhile, most algorithms for nonconvex problems are just particular implementation of the spatial branch and bound framework. Therefore, in order to find the global optimum of a large-scaled nonconvex MINLP models in reasonable computational time, Grossmann and others proposed nonconvex generalized disjunctive programming method.

==General Formulation for GDP==
Consider the following Generalized Disjunctive Programming problem, which includes Boolean variables, disjunctions and logic propositions:<br/>
<math>min Z</math> <math>=</math> <math>\sum_{k\in K} c_k</math> <math>+ f(x)</math> <br/>
<br/>
s.t. <math>r(x)\le 0</math><br/>
<br/>
<math>\bigvee_{j\in J_k}</math> <math>\begin{bmatrix} Y_{jk} \\ g_{jk}(x) \le 0 \\ c_k = \gamma_{jk} \end{bmatrix}</math>, <math> k\in K</math><br/>
<br/>
<math>\Omega(Y) = True </math><br/>
<br/>
<math>0 \le x \le U</math><br/>
<br/>
<math>x \in R^n</math>, <math> c_k \in R^1</math>, <math> Y_{jk} \in {true, false}</math><br/>
<br/>

where, <math>f:R^n \rightarrow R^1</math> is a function of the continuous variables x in the objective function, <math>g:R^n \rightarrow R^1</math> bolongs to the set of global constraints, the disjuctions <math> k\in K</math> are composed of a number of terms <math> j\in J_k</math>, which is connected by the OR operator. <math>Y_{jk}</math> is a Boolean varibale, <math>g_{jk}(x) \le 0</math> is a set of inequalities, and <math>c_k</math> is a cost variable. When <math>Y_{jk}</math> is true, <math>g_{jk}(x) \le 0</math> and <math>c_k</math> are enforced.Also, <math>x</math>represents continuous variables, with lower and upper bounds.Each term in the disjunctions gives rise to a nonempty feasible region which is generally nonconvex. Also, <math>\Omega(Y) = True </math><br/> are logic propositions for the Boolean variables.<br>

==Overall Procedure==
The following flowchart(Fig.1) shows the overall procedure of the proposed two-level branch and bound algorithm.
[[File:Flowchart01_KL.png|600px]]<br/>
First, introduce convex underestimators in the non-convex GDP problem (P), and construct the underestimating problem (R). This convex GDP problem is then reformulated as the convex NLP problem (CRP) by using the convex hull relaxation of each disjunction, which generates valid lower bound. Initial upper bound is obtained in step 0, by sloving P-MIP problem. It is a MINLP reformulation of the nonconvex GDP by a standard MINLP method. Then, use upper bound for bound contraction to reduce the feasible region, which is solved as a Bound Contraction Problem (BCP) in step 1. Next in step 2, discrete branch and bound method is applied at the first level to slove problem (CRP). After fixing all Boolean variables, solve the corresponding nonconvex NLP problems for a upper bound by using spatial branch and bound at the second level. Then, problem (CRP) is solved with fixed discrete choice in step 3.

==Detailed Formulation and Models==
===Convex Relaxation of GDP===
The following reformulation shows the introduction of valid convex underestimating functions to change Problem (P) into a convex GDP problem.<br/>
<br/>
<math>min Z^R</math> <math>=</math> <math>\sum_{k\in K} c_k</math> <math>+ \bar{f}(x)</math> <br/>
<br/>
s.t. <math>\bar{r}(x)\le 0</math><br/>
<br/>
<math>\bigvee_{j\in J_k}</math> <math>\begin{bmatrix} Y_{jk} \\ \bar{g}_{jk}(x) \le 0 \\ c_k = \gamma_{jk} \end{bmatrix}</math>, <math> k\in K</math><br/>
<br/>
<math>\Omega(Y) = True </math><br/>
<br/>
<math>0 \le x \le U</math><br/>,<math>c_k \ge 0</math>
<br/>
<math>x \in R^n</math>, <math>c_k \in R^1</math>, <math> Y_{jk} \in \{ true, false\} </math><br/>
<br/>
where <math>\bar{f}</math>,<math>\bar{r}</math>,<math>\bar{g}</math> are valid convex underestimators so that <math>\bar{f}(x) \le f(x)</math>, <math>\bar{r}(x) \le 0</math>,<math>\bar{g}(x) \le 0</math> are satisfied if <math>r(x) \le 0</math>,<math>g(x) \le 0</math> (see fig.2)<br/>
[[File:Fig02_KL.png|450px]]<br/>
Considering problem (R) is a convex GDP, by replacing each disjunction by its convex hull we can relax problem (R), which generates the folloing convex NLP model:<br/>
<br/>
<math>min Z^L</math> <math>=</math> <math>\sum_{k\in K}\sum_{j\in J_k} \gamma_{jk}\lambda_{jk} + \bar{f}(x)</math> <br/>
<br/>
s.t. <math>\bar{r}(x)\le 0</math><br/>
<br/>
<math> x = \sum_{j \in J_k} \nu_{jk}, \sum_{j \in J_k} \lambda_{jk} = 1,  k\in K</math><br/>
<br/>
<math> 0 \le \nu_{jk} \le \lambda_{jk}U_{jk}, j \in J_k, k\in K</math><br/>
<br/>
<math> \lambda_{jk}\bar{g}_{jk}(\nu_{jk}/\lambda_{jk}) \le 0, j \in J_k, k\in K</math><br/>
<br/>
<math> A\lambda \le a </math>
<br/>
<math>0 \le x, \nu_{jk} \le U, 0 \le \lambda_{jk} \le 1, j \in J_k, k\in K, (CRP)</math><br/>
<br/>
where <math>\nu_{jk}</math> is the disaggregated continuous variable for the  <math>j </math>th term in the  <math>k </math>th disjunction and  <math>\lambda_{jk}</math> is the corresponding multiplier for each term <math>j \in J_k</math> in a given disjunction <math>k\in K</math><br/>
Given the problem (R) yields a lower bound, the problem (CRP) is a relaxation of problem (R).

===Global Upper Bound Subproblem===
This is step is to obtain a valid upper bound for problem (P) based on MINLP reformulation of (P) using big-M formulation.
<br/>
<math>min Z = \sum_{k\in K}\sum_{j\in J_k} \gamma_{jk}y_{jk} + f(x)</math> <br/>
<br/>
s.t. <math>r(x)\le 0</math><br/>
<br/>
<math>g_{jk}(x) \le M_jk(1-y_jk), j \in J_k, k\in K</math><br/>
<br/>
<math> x = \sum_{j \in J_k} y_{jk}= 1,  k\in K</math><br/>
<br/>
<math> Ay \le a </math>
<br/>
<math>0 \le x \le U, y_{jk} \in \{0,1\}, j \in J_k, k\in K, (P-MIP)</math><br/>
<br/>
where <math>M_jk</math> is the big-M parameter, which provides a valid upper bound to the violation of <math>g_jk \le 0</math> and <math>U</math> is an upper bound to <math>x</math>.

===Bound Contraction Procedure===
The bound contraction scheme is designed to tighten the lower and upper bound of a given continuous variable <math>x_i</math>, which will eliminate the nonoptimal subregions and accelerate the search. This will greatly reduce the difference between lower and upper bounds of the objective function and save computational work.<br/>
The contraction scheme is shown by the following NLP problem:<br/>
<br/>
<math>min/max   </math>  <math>x_i  </math> <br/>
<br/>
s.t. <math>\sum_{k\in K}\sum_{j\in J_k} \gamma_{jk}\lambda_{jk} + \bar{f}(x)r(x)\le GUB</math><br/>
<br/>
<math>\bar{r}(x) \le 0 </math><br/>
<br/>
<math> x = \sum_{j \in J_k}\nu_{jk}, \sum_{j \in J_k}\lambda_{jk}= 1,  k\in K</math><br/>
<br/>
<math> 0 \le \nu_{jk} \le \lambda_{jk}U_{jk},j \in J_K, k \in K</math><br/>
<br/>
<math> \lambda_{jk}\bar{g}_{jk}(\nu_{jk}/\lambda_{jk}) \le 0, j \in J_k, k\in K</math><br/>
<br/>
<math> Ay \le a </math>
<br/>
<math>0 \le \nu^{jk} \le U, 0 \le \lambda_{jk} \le 1, j \in J_k, k\in K, (BCP)</math><br/>
<br/>
As is shown in the following Fig.3, if the solution point <math>x_i</math> of problem (CRP) does not lie at its bound, then we solve the NLP problem (CRP) and update the upper and lower bounds based on the relative distance from the solution to each bound of to decide the direction of bound contraction. Bound contraction is applied to only the continuous variables. Not much reduction will be made if the contraction is not successful, and then move to next variable.<br/>
[[File:Fig03_KL.png|300px]]<br/>
===Branch and Bound on Boolean Variables===
This step is where branch and bound method is applied in the space of the terms of the disjunctions by solving the relaxed convex NLP problem (CRP) at each node. The rule is to select the variable <math>\lambda_{jk}</math> with the largest fractional value in the solution. By creating two child nodes with <math>\lambda_{jk}=1</math> and <math>\lambda_{jk}=0</math>, <math>Y_{jk}</math> is fixed in problem (R) respectively. The nodes number should be finite due to finite Boolean variables and the global lower bound increases monotonically. <br/>
Upper bound needs to be updated when there is a gap between the solution of this problem and the original nonconvex GDP problem (P). When the feasible solution to problem (R) is obtained and the gap between nonconvex term and convex underestimator in problem (P) and problem (R) is nonzero, all Boolean variables should be fixed and Spatial Branch and Bound Method should be introduced.<br/>

===Spatial Branch and Bound Method===
Since there is a specific wiki page introducing this method, there is no need to discuss the detailed algorithm here. In this case, by solving problem (CRP), the nonconvex term with the max gap is chosen. Then, choose the middle point of the variable bound as the branching point. Lastly, select the node with lowest objective value. The current node will be fathomed if the lowest objective value of the node is greater than GUB.
==Advantages and Disadvantages==
===Advantages===
Many application shows that tighter lower bounds can be obtained in bilinear and concave problems by applying basic Nonconvex GDP steps, which often leads to a significant reduction of the numbers of nodes in Spatial Branch and Bound method. This is a direct indication of achieving tightening bounds.<br/>
Meanwhile, GDPs are often reformulated as an MINLP, which enables us to take advantages of the existing MINLP solvers.

===Disadvantages===
The main problem about this method is that the question how to efficiently implement the strong relaxations within a spatial branch and bound framework in a large scale system still cannot be answered.

=Applications and Examples=
==Application: Solution Algorithm for nonconvex GDP problems==
===Step 0 Heuristic Search (Nonconvex MINLP)===
Use MINLP solvers (such as DICOPT++) to solve the nonconvex problem (P-MIP), and set GUB as the best upper bound and let <math>(Y^*,x*)</math> be the solution.<br/>
===Step 1 Bound Contraction (Convex NLP)===
# Initialize the Relative Distance from <math>x^*_i</math> to each bound<br/>
<math>RDL_i = \frac{x^*_i - x^L_i}{x^U_i - x^L_i}, RDU_i = \frac{x^U_i - x^*_i}{x^U_i - x^L_i}</math><br/>
# Increase iteration, solve problem (BCP)<br/>
Update bound and continue with <math>x^*_i</math> when contraction is greater than <math>SP_M</math>. Otherwise select next x.<br/>
# Return to 2 and repeat until max iteration.<br/>
===Step 2 Branch and Bound on Discrete Variables (Convex NLP)===
# set tolerance <math>\alpha</math> for difference in <math>Z^L</math> and GUB. Set <math>\varepsilon</math> for max gap.
# Solve problem (CRP) to obtain lower bound <math>Z^L</math>. Update lowest lower bound as GLB.
# Fathom, go to step 3 or Branch on the node.
===Step 3 Spatial Branch and Bound (Convex NLP)===
# Fix all <math>\lambda_{jk}</math> according to solution from step 2
# Solve problem (CRP) until no open node with <math>Z^L \le GUB - \alpha</math>
# Go to step 2

==Example==
===Optimal structure for process system===
The system net work is shown in the following Fig4.<br/>

[[File:Fig04_KL.png|thumb|left|450px|Fig4. Superstructure of the process]]
<br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/>

The formulation for this problem is shown as follows:<br/>
<math> min Z =-35P1A - 30P2B +10F1 + 8F2 + F4A + F4B + 4F5A + 4F5B + 2CF + 50CD</math><br/>
<br/>
<math>s.t.  F3A = 0.55F1 + 0.5F2</math><br/>
<br/>
<math>F3B = 0.45F1 + 0.5F2</math><br/>
<br/>
<math>P1A = F8A + F10A + F6A </math><br/>
<br/>
<math>P1B = F8B + F10B + F6B </math><br/>
<br/>
<math>P2A = F9A + F11A + F7A</math><br/>
<br/>
<math>P2B = F9B +F11B+F7B</math><br/>
<br/>
<math>F6A=E6F3A,F6B=F6F3B</math><br/>
<br/>
<math>F7A=E7F3A,F7B=E7F3B</math><br/>
<br/>
<math>E4+D5+D6+D7 = 1</math><br/>
<br/>
<math>P1A \ge 4.0P1B, P2B \ge 3.0P2A</math><br/>
<br/>
<math>P1A+P1B \le 15, P2A +P2B \le 18</math><br/>
<br/>

<math>\begin{bmatrix}
YF\\
F4A=E4F3A,F4B=E4F3B\\
2.5 \le F4A+F4B \le 25\\
F8A=0.85F4A, F8B=0.20F4B\\
F9A=0.15F4A,F9B=0.8F4B\\
CF=2
\end{bmatrix}</math>
<math>\or</math>
<math>\begin{bmatrix}
\lnot YF\\
F4A=F4B=0\\
F8A=F8B=0\\
F9A=F9B=0\\
E4=0\\
CF=0
\end{bmatrix}</math><br/>
<br/>

<math>\begin{bmatrix}
YD\\
F5A=E5F3A,F5B=E5F3B\\
2.5 \le F5A+F5B \le 25\\
F10A=0.975F5A, F10B=0.050F5B\\
F11A=0.025F5A,F11B=0.950F5B\\
CD=25
\end{bmatrix}</math>
<math>\or</math>
<math>\begin{bmatrix}
\lnot YD\\
F5A=F5B=0\\
F10A=F10B=0\\
F11A=F11B=0\\
E5=0\\
CD=0
\end{bmatrix}</math>
<br/>
<math>0 \le CF,CD; F1,F2 \le 25; 0 \le E4,E5,E6,E7 \le 1;</math><br/>
<br/>
<math>Y\in \{true,false\}</math><br/>
<br/>
The optimal solution is -510.08, and <math> F1^* = 8, F2^* = 25, P1^* = 15, P2^* = 18, E^* = (0.108,0.758,0,0.134)</math>, and <math>Y^*=(true,true)</math>, shown as the following Fig5.<br/>
[[File:Fig05_KL.png|thumb|left|450px|Fig5. Optimal Solution]]
<br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/>

Here is the detailed explanation of the method applied in this example. Bilinear terms are replaced by continuous variables and linear underestimators and overestimators are introduced to construct the relaxed convex GDP problem (R). In step 0, which is mentioned above, DICOPT++ finds the trivial solution of <math>0.0</math> ath first NLP. By solving 14 LP subproblems in step 1, the bound is reduced <math>54.17%</math> for continuous variables. In step 2, discrete branch and bound on problem (CRP) reaches a feasible solution to the relaxed GDP at the third node, where <math>Y^L=(1,1)</math> and objective function of <math>-661.56</math>. When switches to Spatial Branch and Bound in step 3, <math>Y^F=(1,1)</math> is fixed. Then after 27 nodes and gap being reduced, upper bound is found to be <math>-510.08</math>. Using the upper bound to update GUB value in step 2, a cut to exclude <math>Y = (true,true)</math> is added to the previous node and the problem (CRP) is resolved, which leads to a infeasible solution. Then by closing current node and backtrack, two more nodes are searched and fathomed by the GUB.

=Conclusion=
For nonconvex generalized disjunctive programming, specified algorithm can provide a global optimum more efficiently by tightening the lower bound and reduce nodes in Spatial Branch and Bound step, compared with conventional MINLP and GBD algorithms. This method address the increasing need in applications in engineering and other area, where nonlinearity and nonconvexity are essential due to the problem nature.

=References=
# Lee, Sangbum, and Ignacio E. Grossmann. "New algorithms for nonlinear generalized disjunctive programming." Computers & Chemical Engineering 24.9 (2000): 2125-2141.
# Lee, Sangbum, and Ignacio E. Grossmann. "A global optimization algorithm for nonconvex generalized disjunctive programming and applications to process systems." Computers & Chemical Engineering 25.11 (2001): 1675-1697.
# Lee, Sangbum, and Ignacio E. Grossmann. "Global optimization of nonlinear generalized disjunctive programming with bilinear equality constraints: applications to process networks." Computers & chemical engineering 27.11 (2003): 1557-1575.
# Grossmann, Ignacio E., and Juan P. Ruiz. "Generalized Disjunctive Programming: A framework for formulation and alternative algorithms for MINLP optimization." Mixed Integer Nonlinear Programming. Springer New York, 2012. 93-115.