# Reformulation and Convex Relaxation Techniques for Global Optimization

Leo Liberti

<figure><img src="image_1.png" /></figure>

## To cite this version:

Leo Liberti. Reformulation and Convex Relaxation Techniques for Global Optimization. Operations Research [math.OC]. Imperial College London, 2004. English. NNT: . hal-00163562

HAL Id: hal-00163562

https://hal.science/hal-00163562v1

Submitted on 14 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Reformulation and Convex Relaxation Techniques
for Global Optimization

LEO SERGIO LIBERTI

15th March 2004

A thesis submitted for the degree of Doctor of Philosophy of the University of London
and for the Diploma of Imperial College

Department of Chemical Engineering and Chemical Technology

Imperial College London

South Kensington Campus

London SW7 2AZ
Abstract

Many engineering optimization problems can be formulated as nonconvex nonlinear programming problems (NLPs) involving a nonlinear objective function subject to nonlinear constraints. Such problems may exhibit more than one locally optimal point. However, one is often solely or primarily interested in determining the globally optimal point. This thesis is concerned with techniques for establishing such global optima using spatial Branch-and-Bound (sBB) algorithms.

A key issue in optimization is that of mathematical formulation, as there may be several different ways in which the same engineering problem can be expressed mathematically. This is particularly important in the case of global optimization as the solution of different mathematically equivalent formulations may pose very different computational requirements. Based on the concept of reduction constraints, the thesis presents a set of graph-theoretical algorithms which automatically reformulate large sparse nonconvex NLPs involving linear constraints and bilinear terms. It is shown that the resulting exact reformulations involve fewer bilinear terms and have tighter convex relaxations than the original NLPs. Numerical results illustrating the beneficial effects of applying such automatic reformulations to the well-known pooling and blending problem are presented.

All sBB algorithms rely on the construction of a convex relaxation of the original NLP problem. Relatively tight convex relaxations are known for many categories of algebraic expressions. One notable exception is that of monomials of odd degree, i.e. expressions of the form $x^{2n+1}$ where $n \ge 1$, when the range of the variable $x$ includes zero. These occur often (e.g. as cubic or quintic expressions) in practical applications. The thesis presents a novel method for constructing the convex envelope of such monomials, as well as a tight linear relaxation of this envelope.

Finally, the thesis discusses some of the software engineering issues involved in the design and implementation of codes for sBB, especially in view of the large amounts of both symbolic and numerical information required by these codes. A prototype object-oriented software library, *ooOPS*, is described.
Acknowledgments

After nearly four years of research, it is hard to think whom *not* to thank in this Ph.D. thesis.
Besides, I feel I owe everyone who has been around me during my stay in London a gesture of
thanks. However, this would result into a short biography of my last four years rather than a
Ph.D. acknowledgment.

I naturally want to thank my supervisor Costas Pantelides for his support, both technical and
financial. I would also like to thank: Panagiotis Tsiakis and Benjamin Keeping for produc-
ing the *ooMILP* software, on which my *ooOPS* software is heavily based; Fabrizio Bezzo
for listening to many of my woes and lamentations about my latest theoretical ideas leading
to nothing, and helping whenever he could; Gerard Gorman for our lengthy talks about C++,
coding and debugging techniques; and Panagiotis Karamertzanis for crucial help with the bu-
reaucratic pitfalls of thesis submission. A lot of other people, both in the CPSE and elsewhere,
have helped me in my research, but the list is just too long to include it here. As too long would
be the list of all people who have extended their moral support and friendship to me during these
years: so please excuse me if you are not on this page. You are in my memory nonetheless.

Finally, I would like to thank my family and three dearest friends, Jacopo, Giovanna and
Cecilia, for their love and support in times of need.
<table>
  <tbody>
    <tr>
      <td>2.1.2</td>
      <td>Separable problems</td>
      <td>35</td>
    </tr>
    <tr>
      <td>2.1.2.1</td>
      <td>Separation of bilinear forms</td>
      <td>35</td>
    </tr>
    <tr>
      <td>2.1.2.2</td>
      <td>Global solution of separable box-constrained problems</td>
      <td>36</td>
    </tr>
    <tr>
      <td>2.1.3</td>
      <td>Linear problems</td>
      <td>37</td>
    </tr>
    <tr>
      <td>2.1.3.1</td>
      <td>Reformulating quadratic binary problems to linear binary problems</td>
      <td>37</td>
    </tr>
    <tr>
      <td>2.1.4</td>
      <td>Convex problems</td>
      <td>38</td>
    </tr>
    <tr>
      <td>2.1.5</td>
      <td>Binary problems</td>
      <td>39</td>
    </tr>
    <tr>
      <td>2.1.5.1</td>
      <td>Reformulating discrete problems to binary problems</td>
      <td>39</td>
    </tr>
    <tr>
      <td>2.1.5.2</td>
      <td>Reformulating binary problems to continuous problems</td>
      <td>40</td>
    </tr>
    <tr>
      <td>2.1.6</td>
      <td>Concave problems</td>
      <td>41</td>
    </tr>
    <tr>
      <td>2.1.6.1</td>
      <td>Reformulating binary problems to concave problems</td>
      <td>41</td>
    </tr>
    <tr>
      <td>2.1.6.2</td>
      <td>Reformulating bilinear problems to concave problems</td>
      <td>41</td>
    </tr>
    <tr>
      <td>2.1.6.3</td>
      <td>Reformulating complementarity problems to concave problems</td>
      <td>42</td>
    </tr>
    <tr>
      <td>2.1.6.4</td>
      <td>Reformulating max-min problems to concave problems</td>
      <td>42</td>
    </tr>
    <tr>
      <td>2.1.7</td>
      <td>D.c. problems</td>
      <td>43</td>
    </tr>
    <tr>
      <td>2.1.7.1</td>
      <td>Reformulating continuous functions to d.c. functions</td>
      <td>43</td>
    </tr>
    <tr>
      <td>2.1.8</td>
      <td>Factorable problems</td>
      <td>44</td>
    </tr>
    <tr>
      <td>2.1.8.1</td>
      <td>Reformulation of factorable problems to separable form</td>
      <td>45</td>
    </tr>
    <tr>
      <td>2.1.9</td>
      <td>Smith's standard form</td>
      <td>46</td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>Exact reformulations</td>
      <td>47</td>
    </tr>
    <tr>
      <td>2.2.1</td>
      <td>Equality/inequality constrained problems</td>
      <td>47</td>
    </tr>
    <tr>
      <td>2.2.2</td>
      <td>Dimensionality reduction</td>
      <td>47</td>
    </tr>
    <tr>
      <td>2.3</td>
      <td>Convex relaxations</td>
      <td>49</td>
    </tr>
    <tr>
      <td>2.3.1</td>
      <td>αBB convex relaxation</td>
      <td>49</td>
    </tr>
  </tbody>
</table>
<table><tr><td>2.3.2</td><td>Smith's convex relaxation</td><td>53</td></tr><tr><td>2.3.3</td><td>BARON's convex relaxation</td><td>53</td></tr><tr><td>2.4</td><td>Linear relaxations</td><td>56</td></tr><tr><td>2.4.1</td><td>Reformulation-linearization technique</td><td>56</td></tr><tr><td>2.5</td><td>Other reformulations</td><td>58</td></tr><tr><td>2.6</td><td>Conclusion</td><td>59</td></tr><tr><td><b>3</b></td><td><b>Reduction constraints for sparse bilinear programs</b></td><td><b>60</b></td></tr><tr><td>3.1</td><td>Basic concepts</td><td>62</td></tr><tr><td>3.2</td><td>Fundamental properties</td><td>65</td></tr><tr><td>3.3</td><td>An algorithm for the identification of valid reduction constraints</td><td>67</td></tr><tr><td>3.3.1</td><td>Valid reduction constraint sets</td><td>67</td></tr><tr><td>3.3.2</td><td>A graph-theoretical representation of linear constraints</td><td>68</td></tr><tr><td>3.3.3</td><td>Efficient identification of dilations</td><td>69</td></tr><tr><td>3.4</td><td>A detailed example</td><td>73</td></tr><tr><td>3.4.1</td><td>Valid reduction constraints by multiplication with variable z<sub>1</sub></td><td>74</td></tr><tr><td>3.4.2</td><td>Valid reduction constraints by multiplication with variable z<sub>2</sub></td><td>76</td></tr><tr><td>3.4.3</td><td>Valid reduction constraints by multiplication with variable z<sub>3</sub></td><td>77</td></tr><tr><td>3.4.4</td><td>Valid reduction constraints by multiplication with variables z<sub>4</sub>, ..., z<sub>6</sub></td><td>78</td></tr><tr><td>3.4.5</td><td>The reformulated NLP</td><td>78</td></tr><tr><td>3.5</td><td>Computational results</td><td>80</td></tr><tr><td>3.6</td><td>Generalization of the graph-theoretical algorithm</td><td>85</td></tr><tr><td>3.7</td><td>Concluding remarks</td><td>88</td></tr><tr><td><b>4</b></td><td><b>A convex relaxation for monomials of odd degree</b></td><td><b>90</b></td></tr></table>
<table><tr><td>4.1</td><td>Statement of the problem</td><td>91</td></tr><tr><td>4.2</td><td>The tangent equations</td><td>92</td></tr><tr><td>4.3</td><td>The roots of Q<sup>k</sup>(x) and their uniqueness</td><td>93</td></tr><tr><td>4.3.1</td><td>Bounding the roots of Q<sup>k</sup>(x)</td><td>94</td></tr><tr><td>4.4</td><td>Nonlinear convex envelopes</td><td>96</td></tr><tr><td>4.5</td><td>Tight linear relaxation</td><td>97</td></tr><tr><td>4.6</td><td>Comparison to other relaxations</td><td>98</td></tr><tr><td>4.6.1</td><td>Reformulation in terms of bilinear products</td><td>99</td></tr><tr><td>4.6.2</td><td>Underestimation through α parameter</td><td>100</td></tr><tr><td>4.7</td><td>Computational results</td><td>101</td></tr><tr><td>4.8</td><td>Conclusion</td><td>102</td></tr></table>

5 Spatial Branch-and-Bound algorithm with symbolic reformulation 105

<table><tr><td>5.1</td><td>Overview of spatial Branch-and-Bound algorithms</td><td>106</td></tr><tr><td>5.2</td><td>Smith's sBB algorithm</td><td>107</td></tr><tr><td>5.2.1</td><td>Choice of region (step 2)</td><td>107</td></tr><tr><td>5.2.2</td><td>Convex relaxation (step 3)</td><td>108</td></tr><tr><td>5.2.2.1</td><td>Reformulation to standard form</td><td>108</td></tr><tr><td>5.2.2.2</td><td>Convexification</td><td>109</td></tr><tr><td>5.2.3</td><td>Branching (step 7)</td><td>110</td></tr><tr><td>5.2.4</td><td>Bounds tightening</td><td>111</td></tr><tr><td>5.2.4.1</td><td>Optimization-based bounds tightening (step 1)</td><td>111</td></tr><tr><td>5.2.4.2</td><td>Feasibility-based bounds tightening (step 2)</td><td>112</td></tr><tr><td>5.3</td><td>Improvements to Smith's sBB algorithm</td><td>112</td></tr><tr><td>5.3.1</td><td>Avoiding the introduction of slack variables</td><td>112</td></tr></table>
<table><tr><td>5.3.2 Avoiding unnecessary local optimizations</td><td>113</td></tr><tr><td>5.3.2.1 Branching on added variables</td><td>114</td></tr><tr><td>5.3.2.2 Branching on original variables</td><td>114</td></tr><tr><td>5.4 The <i>ooOPS</i> software framework for optimization</td><td>114</td></tr><tr><td>5.4.1 Main features of <i>ooOPS</i></td><td>116</td></tr><tr><td>5.4.2 Object classes in <i>ooOPS</i></td><td>118</td></tr><tr><td>5.4.3 Typical usage scenario for <i>ooOPS</i></td><td>119</td></tr><tr><td>5.4.4 Solver components for use within <i>ooOPS</i></td><td>120</td></tr><tr><td>5.5 An sBB solver for the <i>ooOPS</i> framework</td><td>121</td></tr><tr><td>5.5.1 Overview of the sBB code</td><td>121</td></tr><tr><td>5.5.2 The convexifier</td><td>122</td></tr><tr><td>5.5.3 Storing the list of regions</td><td>123</td></tr><tr><td>5.6 Concluding remarks</td><td>125</td></tr><tr><td><b>6 Concluding remarks</b></td><td><b>127</b></td></tr><tr><td><b>Bibliography</b></td><td><b>130</b></td></tr><tr><td><b>A <i>ooOPS</i> Reference Manual</b></td><td><b>142</b></td></tr><tr><td>A.1 Introduction</td><td>142</td></tr><tr><td>A.2 Fundamental concepts</td><td>143</td></tr><tr><td>A.2.1 Object classes</td><td>143</td></tr><tr><td>A.2.2 Multidimensional sets in <i>ooOPS</i></td><td>143</td></tr><tr><td>A.2.3 Constants, variables, constraints and objective function</td><td>144</td></tr><tr><td>A.2.4 Nonlinear expressions and constants</td><td>145</td></tr></table>
<table><tr><td>A.2.5</td><td>Notes</td><td>146</td></tr><tr><td>A.3</td><td>General software issues</td><td>147</td></tr><tr><td>A.3.1</td><td>Software constants</td><td>147</td></tr><tr><td>A.3.2</td><td>The si, sd and li Argument Types and the IntSeq Auxiliary Function</td><td>147</td></tr><tr><td>A.4</td><td>The MINLP object class</td><td>148</td></tr><tr><td>A.4.1</td><td>MINLP instantiation: the NewMINLP function</td><td>148</td></tr><tr><td>A.4.2</td><td>MINLP construction methods</td><td>148</td></tr><tr><td>A.4.2.1</td><td>Method NewContinuousVariable</td><td>148</td></tr><tr><td>A.4.2.2</td><td>Method NewIntegerVariable</td><td>149</td></tr><tr><td>A.4.2.3</td><td>Method NewConstraint</td><td>150</td></tr><tr><td>A.4.2.4</td><td>AddLinearVariableSliceToConstraintSlice</td><td>151</td></tr><tr><td>A.4.2.5</td><td>Method NewConstant</td><td>152</td></tr><tr><td>A.4.2.6</td><td>Method NewConstantExpression</td><td>153</td></tr><tr><td>A.4.2.7</td><td>Method NewVariableExpression</td><td>153</td></tr><tr><td>A.4.2.8</td><td>Method BinaryExpression</td><td>154</td></tr><tr><td>A.4.2.9</td><td>Method UnaryExpression</td><td>154</td></tr><tr><td>A.4.2.10</td><td>AssignExpressionSliceToConstraintSlice</td><td>155</td></tr><tr><td>A.4.2.11</td><td>Method NewObjectiveFunction</td><td>156</td></tr><tr><td>A.4.2.12</td><td>AddLinearVariableSliceToObjectiveFunction</td><td>157</td></tr><tr><td>A.4.2.13</td><td>Method AssignExpressionToObjectiveFunction</td><td>157</td></tr><tr><td>A.4.3</td><td>MINLP modification methods</td><td>158</td></tr><tr><td>A.4.3.1</td><td>Method SetVariableValue</td><td>158</td></tr><tr><td>A.4.3.2</td><td>Method SetVariableBounds</td><td>159</td></tr><tr><td>A.4.3.3</td><td>Method SetConstraintBounds</td><td>160</td></tr></table>
<table><tr><td>A.4.3.4</td><td>Method SetConstantValue</td><td>160</td></tr><tr><td>A.4.3.5</td><td>Method SetConstantSliceValue</td><td>161</td></tr><tr><td>A.4.3.6</td><td>Method SetKeyVariable</td><td>161</td></tr><tr><td>A.4.4</td><td>Structured MINLP Information Access Methods</td><td>162</td></tr><tr><td>A.4.4.1</td><td>Method GetVariableInfo</td><td>162</td></tr><tr><td>A.4.4.2</td><td>Method GetConstraintInfo</td><td>163</td></tr><tr><td>A.4.4.3</td><td>Method GetObjectiveFunctionInfo</td><td>164</td></tr><tr><td>A.4.4.4</td><td>Method GetProblemInfo</td><td>165</td></tr><tr><td>A.4.5</td><td>Flat MINLP Information Access Methods</td><td>165</td></tr><tr><td>A.4.5.1</td><td>General Properties of the Flat MINLP</td><td>165</td></tr><tr><td>A.4.5.2</td><td>Method GetFlatMINLPSize</td><td>166</td></tr><tr><td>A.4.5.3</td><td>Method GetFlatMINLPStructure</td><td>167</td></tr><tr><td>A.4.5.4</td><td>Method GetFlatMINLPVariableInfo</td><td>168</td></tr><tr><td>A.4.5.5</td><td>Method SetFlatMINLPVariableBounds</td><td>169</td></tr><tr><td>A.4.5.6</td><td>Method GetFlatMINLPVariableValues</td><td>170</td></tr><tr><td>A.4.5.7</td><td>Method SetFlatMINLPVariableValues</td><td>170</td></tr><tr><td>A.4.5.8</td><td>Method GetFlatMINLPConstraintInfo</td><td>171</td></tr><tr><td>A.4.5.9</td><td>EvalFlatMINLPNonlinearObjectiveFunction</td><td>172</td></tr><tr><td>A.4.5.10</td><td>Method EvalFlatMINLPNonlinearConstraint</td><td>172</td></tr><tr><td>A.4.5.11</td><td>GetFlatMINLPObjectiveFunctionDerivatives</td><td>173</td></tr><tr><td>A.4.5.12</td><td>Method GetFlatMINLPConstraintDerivatives</td><td>174</td></tr><tr><td>A.4.5.13</td><td>Method GetFlatMINLPNoPartitions</td><td>175</td></tr><tr><td>A.4.5.14</td><td>Method GetFlatMINLPPartition</td><td>175</td></tr><tr><td>A.4.6</td><td>Standard Form MINLP Information Access Methods</td><td>176</td></tr></table>
<table>
  <tr>
    <td>A.4.6.1</td>
    <td>Method GetSFNumberOfVariables</td>
    <td>176</td>
  </tr>
  <tr>
    <td>A.4.6.2</td>
    <td>Method GetSFVariableInfo</td>
    <td>177</td>
  </tr>
  <tr>
    <td>A.4.6.3</td>
    <td>Method GetSFObjFunVarIndex</td>
    <td>177</td>
  </tr>
  <tr>
    <td>A.4.6.4</td>
    <td>Method GetSFNumberOfLinearConstraints</td>
    <td>178</td>
  </tr>
  <tr>
    <td>A.4.6.5</td>
    <td>Method GetSFLinearBounds</td>
    <td>178</td>
  </tr>
  <tr>
    <td>A.4.6.6</td>
    <td>Method GetSFLinearStructure</td>
    <td>179</td>
  </tr>
  <tr>
    <td>A.4.6.7</td>
    <td>Method GetSFMatrix</td>
    <td>179</td>
  </tr>
  <tr>
    <td>A.4.6.8</td>
    <td>Method GetSFNumberOfNonlinearConstraints</td>
    <td>180</td>
  </tr>
  <tr>
    <td>A.4.6.9</td>
    <td>Method GetSFNonlinearConstraint</td>
    <td>180</td>
  </tr>
  <tr>
    <td>A.4.6.10</td>
    <td>Method UpdateSolution</td>
    <td>181</td>
  </tr>
  <tr>
    <td>A.5</td>
    <td>MINLP solvers and systems</td>
    <td>182</td>
  </tr>
  <tr>
    <td>A.5.1</td>
    <td>Introduction</td>
    <td>182</td>
  </tr>
  <tr>
    <td>A.5.1.1</td>
    <td>MINLP solver managers and MINLP systems</td>
    <td>182</td>
  </tr>
  <tr>
    <td>A.5.1.2</td>
    <td>Algorithmic parameters for MINLP solvers</td>
    <td>182</td>
  </tr>
  <tr>
    <td>A.5.2</td>
    <td>The ssolpar and sstat argument types</td>
    <td>182</td>
  </tr>
  <tr>
    <td>A.5.3</td>
    <td>NewMINLPSolverManager</td>
    <td>183</td>
  </tr>
  <tr>
    <td>A.5.4</td>
    <td>MINLP solver managers</td>
    <td>184</td>
  </tr>
  <tr>
    <td>A.5.4.1</td>
    <td>Method GetParameterList</td>
    <td>184</td>
  </tr>
  <tr>
    <td>A.5.4.2</td>
    <td>Method SetParameter</td>
    <td>185</td>
  </tr>
  <tr>
    <td>A.5.4.3</td>
    <td>Method NewMINLPSystem</td>
    <td>185</td>
  </tr>
  <tr>
    <td>A.5.5</td>
    <td>MINLP systems</td>
    <td>186</td>
  </tr>
  <tr>
    <td>A.5.5.1</td>
    <td>Method GetParameterList</td>
    <td>186</td>
  </tr>
  <tr>
    <td>A.5.5.2</td>
    <td>Method SetParameter</td>
    <td>186</td>
  </tr>
  <tr>
    <td>A.5.5.3</td>
    <td>Method GetStatistics</td>
    <td>187</td>
  </tr>
</table>
<table><tr><td>A.5.5.4</td><td>Method Solve</td><td>187</td></tr><tr><td>A.5.5.5</td><td>Method GetSolutionStatus</td><td>188</td></tr><tr><td>A.6</td><td>Auxiliary interfaces</td><td>188</td></tr><tr><td>A.6.1</td><td>The convexification module</td><td>188</td></tr><tr><td>A.6.1.1</td><td>Convexifier manager instantiation: the function NewConvexifierManager</td><td>189</td></tr><tr><td>A.6.1.2</td><td>Method GetConvexMINLP</td><td>189</td></tr><tr><td>A.6.1.3</td><td>Method UpdateConvexVarBounds</td><td>190</td></tr><tr><td>A.6.1.4</td><td>Methods of the convex MILP</td><td>191</td></tr><tr><td>A.6.2</td><td>The FlatExpression interface</td><td>192</td></tr><tr><td>A.6.2.0.1</td><td>Method GetKind</td><td>192</td></tr><tr><td>A.6.2.1</td><td>The FlatConstantExpression interface</td><td>193</td></tr><tr><td>A.6.2.1.1</td><td>Method GetValue</td><td>193</td></tr><tr><td>A.6.2.2</td><td>The FlatVariableExpression interface</td><td>193</td></tr><tr><td>A.6.2.2.1</td><td>Method GetVarIndex</td><td>193</td></tr><tr><td>A.6.2.3</td><td>The FlatOperatorExpression interface</td><td>194</td></tr><tr><td>A.6.2.3.1</td><td>Method GetOperator</td><td>194</td></tr><tr><td>A.6.2.3.2</td><td>Method GetOperand</td><td>194</td></tr><tr><td>A.6.2.4</td><td>Usage of FlatExpression interface</td><td>195</td></tr><tr><td>A.7</td><td>Implementation restrictions</td><td>196</td></tr><tr><td>A.8</td><td>An Example of the use of ooOPS</td><td>196</td></tr><tr><td>A.8.1</td><td>Creating the MINLP</td><td>198</td></tr><tr><td>A.8.1.1</td><td>Creating variables</td><td>198</td></tr><tr><td>A.8.1.2</td><td>Creating constraints</td><td>198</td></tr><tr><td>A.8.1.3</td><td>Adding variables to constraints</td><td>199</td></tr></table>
<table><tr><td>A.8.1.4</td><td>Objective function</td><td>200</td></tr><tr><td>A.8.1.5</td><td>Objective function coefficients</td><td>200</td></tr><tr><td>A.8.1.6</td><td>Modifying the variable bounds</td><td>201</td></tr><tr><td>A.8.1.7</td><td>Modifying the constraint bounds</td><td>201</td></tr><tr><td>A.8.1.8</td><td>Creating the nonlinear parts</td><td>201</td></tr><tr><td>A.8.1.9</td><td>Assigning expressions to constraints and objective function</td><td>202</td></tr><tr><td>A.8.2</td><td>MINLP solution</td><td>202</td></tr><tr><td>A.8.2.1</td><td>Creating an MINLP solver manager object</td><td>203</td></tr><tr><td>A.8.2.2</td><td>Creating an MINLP system</td><td>203</td></tr><tr><td>A.8.2.3</td><td>Solving the MINLP</td><td>203</td></tr><tr><td>A.8.3</td><td>Accessing the Solution of the MINLP</td><td>203</td></tr><tr><td>A.8.3.1</td><td>Obtaining information on the variables</td><td>203</td></tr><tr><td>A.8.3.2</td><td>Obtaining information on the objective function</td><td>204</td></tr></table>
List of Figures

<table>
  <tbody>
    <tr>
      <td>1.1</td>
      <td>Fathoming via upper bound computation.</td>
      <td>27</td>
    </tr>
    <tr>
      <td>2.1</td>
      <td>The Hilbert space-filling curve in 2-dimensional space.</td>
      <td>48</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Linearity embedded in a bilinear constraint.</td>
      <td>64</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Algorithm for the construction of an augmenting path emanating from a linear constraint node λ in graph G<sub>l</sub>.</td>
      <td>71</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Algorithm for identification of set of valid reduction constraint set L for variable z<sub>l</sub>.</td>
      <td>72</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>Generation of valid reduction constraints for multiplier variable z<sub>1</sub>.</td>
      <td>74</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Generation of valid reduction constraints for multiplier variable z<sub>2</sub>.</td>
      <td>76</td>
    </tr>
    <tr>
      <td>3.6</td>
      <td>Generation of valid reduction constraints for multiplier variable z<sub>3</sub>.</td>
      <td>77</td>
    </tr>
    <tr>
      <td>3.7</td>
      <td>Unified bipartite graph for problem (3.34) in the generalized algorithm.</td>
      <td>87</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Convex envelope of x<sup>2k+1</sup>.</td>
      <td>91</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>The case when c > b.</td>
      <td>92</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Tight linear relaxation of x<sup>2k+1</sup>.</td>
      <td>98</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>Convex relaxation of x<sup>3</sup> by reformulation to bilinear product.</td>
      <td>100</td>
    </tr>
    <tr>
      <td>4.5</td>
      <td>Convex relaxation of x<sup>3</sup> by the α method.</td>
      <td>101</td>
    </tr>
    <tr>
      <td>4.6</td>
      <td>Graphical description of simple test problem for k = 1 (in 2D).</td>
      <td>102</td>
    </tr>
  </tbody>
</table>
4.7 Graphical description of simple test problem for $k = 1$ (in 3D). . . . . . . . . . 103

5.1 If the locally optimal solution in $[x^L, x^U]$ has already been determined to be at $x^*$, solving in $[x^L, x']$ is unnecessary. . . . . . . . . . . . . . . . . . . . . . . . . . 115

5.2 The region tree. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
# List of Tables

<table><tbody><tr><td>3.1 Test problem statistics.</td><td>83</td></tr><tr><td>3.2 Numerical results.</td><td>84</td></tr><tr><td>4.1 Numerical values of the roots of Q<sup>k</sup>(x) for k = 1, .., 10 (to 10 significant digits).</td><td>93</td></tr><tr><td>4.2 Summary of linear relaxations for z = x<sup>2k+1</sup>, x ∈ [a, b], a &lt; 0 &lt; b.</td><td>98</td></tr><tr><td>4.3 Numerical results from the simple test problem.</td><td>104</td></tr><tr><td>A.1 Variable sets</td><td>197</td></tr><tr><td>A.2 Parameters appearing in the RTN formulation</td><td>197</td></tr></tbody></table>
# Chapter 1

## Introduction

Optimization is a branch of mathematics that studies the problem of finding the best choice among a set of entities satisfying some imposed requirements.

In general, the problem is formulated in terms of finding the point $x$ in a space set $\Omega$ (called the feasible region) where a certain function $f: \Omega \to T$ (called the objective function), attains a minimum or a maximum. $T$ is any ordered set.

### 1.1 Basic definitions

In this section we shall introduce some basic definitions which will be used throughout this thesis.

The set $\Omega \subseteq \mathbb{R}^n$ is convex if, for all $x, y \in \Omega$ and for all $\lambda \in [0, 1]$, the vector $\lambda x + (1 - \lambda)y$ is also in $\Omega$. The intersection of an arbitrary collection of convex sets is convex. Let $S \subseteq \mathbb{R}^n$. The intersection of all the convex subsets of $\mathbb{R}^n$ containing $S$ is called the convex hull of $S$. The convex hull of a finite subset $\{v_1, \dots, v_m\}$ of $\mathbb{R}^n$ consists of all the linear combinations $\sum_{i=1}^m \lambda_i v_i$ where $\lambda_i \ge 0$ for all $i \le m$ and $\sum_{i=1}^m \lambda_i = 1$. A function $f: \mathbb{R}^n \to \mathbb{R}$ is convex if and only if it is defined on a convex set $\Omega$ and is such that for all $x, y \in \Omega$ and for all $\lambda \in [0, 1]$ we have

$$f(\lambda x + (1 - \lambda)y) \le \lambda f(x) + (1 - \lambda)f(y).$$

A function $f$ is concave if $-f$ is convex.
A function $f$ is pseudo-convex if, for all $x_1, x_2$ such that $f(x_1) < f(x_2)$, we have $\nabla f(x_2)(x_1 - x_2) < 0$.

A function $f$ is quasi-conovex if all its sublevel sets $S_\alpha = \{x \mid f(x) \le \alpha\}$ are convex.

A function $f: \Omega \to \mathbb{R}$ is a d.c. function if it is a difference of convex functions, i.e. there exist convex functions $g, h: \Omega \to \mathbb{R}$ such that, for all $x \in \Omega$, we have $f(x) = g(x) - h(x)$. Let $C, D$ be convex sets; then the set $C\setminus D$ is a d.c. set<fnref n="1" />. It can be shown that a set $M$ for which there exist convex functions $g, h: \mathbb{R}^n$ such that $M = \{x \mid g(x) \le 0 \land h(x) \ge 0\}$ is a d.c. set. D.c. functions and sets have many interesting properties. See Section 2.1.7 for a more thorough discussion.

Let $f(x)$ be a nonconvex function. A convex relaxation of $f(x)$ is a convex function $\underline{f}(x)$ such that, for all $x$, we have $\underline{f}(x) \le f(x)$. Likewise, a concave relaxation of $f(x)$ is a concave function $\bar{f}(x)$ such that for all $x$ we have $\bar{f}(x) \ge f(x)$. Let $\underline{F}$ be the set of all convex relaxations of $f$ and $\bar{F}$ the set of all concave relaxations. The convex envelope of $f$ is $\max\{g(x) \mid g \in \underline{F}\}$; the concave envelope of $f$ is $\min\{g(x) \mid g \in \bar{F}\}$.

Given an inequality $g(x) \le 0$, the set $\{x \mid g(x) \le 0\}$ is convex if $g(x)$ is a convex function. Notice that, if $g(x)$ is concave, then $g(x) \ge 0$ is a convex inequality.

## 1.2 Classification of optimization problems

The field of optimization can be subdivided into many categories according to the different possible formulations of the problem.

* Classification with respect to the type of optimization decisions.

- **Continuous optimization.** The set $\Omega$ is a subset of a Euclidean space.

- **Integer optimization.** $\Omega$ is finite or countable.

- **Mixed integer optimization.** $\Omega = U \times V$ where $U$ is a subset of a Euclidean space and $V$ is finite or countable.

* Classification with respect to the type of objective function and constraints.

- **Constrained optimization.** $\Omega$ is defined by a set of constraints $C$ to be satisfied.

<footnote n="1">¹I.e. a difference of convex sets; here $C\setminus D$ is taken to be all the elements in $C$ which are not in $D$.</footnote>
- **Unconstrained optimization.** The set of constraints $C$ to be satisfied is empty. However, problems with variable ranges as their only constraints are sometimes also classified as “unconstrained problems”.

- **Linear optimization.** The function $f$ and the constraints in $C$ are linear functions of $x$.

- **Nonlinear optimization.** The function $f$ and the constraints in $C$ may be nonlinear functions of $x$.

- **Convex optimization.** The function $f$ is convex and feasible region $\Omega$ is convex. In particular, a problem with convex constraints belongs to this class.

- **Concave optimization.** The function $f$ is concave (usually, but not always, the constraints in $C$ are convex).

- **D.c. optimization.** The function $f$ is d.c. and the feasible region $\Omega$ is a d.c. set.

- **Nonconvex optimization.** The function $f$ and the constraints in $C$ may be general nonconvex functions.

* Classification with respect to the type of solution obtained.

- **Local optimization.** The solution $x$ is such that there is a subset $V \subseteq \Omega$, such that $x \in V$, where $\forall y \in V$ ($f(x) \le f(y)$). The point $x$ is called a *local minimizer*, or a *local optimum* of $f$ with respect to $V$. It is nearly always the case that $V$ is a neighbourhood of $x$ in the topological sense.

- **Global optimization.** The solution $x$ is such that $\forall y \in \Omega$ ($f(x) \le f(y)$). $x$ is called a *global minimizer*, or a *global optimum* of $f$.

One of the most interesting and most challenging classes of problems is global nonlinear constrained optimization. Problems belonging to this class are called nonlinear programs (NLPs).

## 1.3 Algorithms for global optimization of NLPs

Algorithmic methods for globally solving NLPs are usually divided into two main categories: deterministic and nondeterministic (or stochastic). In deterministic optimization, the solution
method never employs random choices, and convergence theorems do not use probability arguments. In stochastic optimization<fnref n="2" /> algorithms, the solution methods employ random<fnref n="3" /> choices. Convergence is proved through arguments based on probability.

### 1.3.1 Formulation of the NLP

In this thesis, we are interested in constrained NLPs of the following form:

$$ \left. \begin{array}{ll} \min_{x} & f(x) \\ \alpha & \leq g(x) \leq \beta \\ a & \leq x \leq b \end{array} \right\} \qquad (1.1) $$

where $x \in \mathbb{R}^n$ are the (continuous) decision variables, $\alpha, \beta$ are the lower and upper bounds of the constraints, and $a, b$ are the lower and upper bounds of the variables. The function $f: \mathbb{R}^n \to \mathbb{R}$ is called the objective function and $g: \mathbb{R}^n \to \mathbb{R}^m$ are the constraints of the problem. We limit the discussion to the minimization of the problem; the maximization is equivalent to solving $\min_x -f(x)$ subject to the same constraints.

Notice that formulation (1.1) is slightly unusual because of the form of the constraints $\alpha \leq g(x) \leq \beta$. The most common constraint formulation in the literature is $h(x) = 0 \land g(x) \leq 0$. We employ the former because it makes the software representation of an optimization problem more compact.

### 1.3.2 A brief history of global optimization

Generic optimization problems have been important throughout history in engineering applications. The first significant work in optimization was carried out by Lagrange in 1797 [66]. Nevertheless, before the introduction and extensive usage of electronic computers, the requirement for global optimization was not even an issue given the tremendous amount of computational effort it requires. Global optimality of solutions was guaranteed only when locally optimizing convex functions, a rather limited class of problems.

The methods that were first used in global optimization were deterministic techniques, mostly based on the divide-and-conquer principle. This was introduced in the late 1950s with the ad-

<footnote n="2">²Albeit an abuse of terminology, the term “nondeterministic optimization” is sometimes also used.</footnote>

<footnote n="3">³The term “random” is used here in the special sense of pseudo-random.</footnote>
vent of the first electronic computers into the research environment. In contrast with the com-
putational principles employed before the introduction of computers, whereby people would
try to minimize the amount of actual computations to perform with respect to the theoretical
framework, the divide-and-conquer principle applies to the intrinsic iterative nature of methods
used when working with electronic computers: keep the complexity of the theoretical struc-
tures involved to a minimum, and rely on intensive computation to explore the solution space.
Thus, instead of trying to locate a minimum by solving a set of equations by symbolic/algebraic
methods, one would try to construct a sequence of approximate solutions which would converge
to the true solution, by dividing the problem into smaller subproblems. One typical algorithm
which embodies the divide-and-conquer principle is the Branch-and-Bound algorithm (BB).
Because of the nature of the algorithm, where the subproblems are produced by branching a
problem entity (e.g. variable) into its possible instances, the BB algorithm applies very well to
cases where problem entities are discrete in nature. Thus, the first applications of BB to global
optimization problems were devoted to discrete problems such as the Travelling Salesman Prob-
lem (TSP).

The first paper concerning continuous global optimization with a BB (deterministic) tech-
nique dates from 1969 [34]. In the 1970s and 1980s, work on continuous or mixed-integer
deterministic global optimization was scarce. Most of the papers published in this period dealt
either with applications of global optimization to very specific cases, or with theoretical results
concerning convergence proofs. One notable exception was the work of McCormick [80] who
considered symbolic transformations of problems: his methods are such that they can, in theory,
be carried out automatically by a computer.

A major reason for the slow pace of progress in continuous global optimization is that it is
computationally very expensive, and it was not until the 1990s that computer hardware with the
necessary power became available. Toward the end of the 1980s, however, the first elementary
textbooks began to appear on the subject. They explained the basics of local optimization (La-
grange multipliers and Karush-Kuhn-Tucker conditions) and some of the early techniques in
deterministic [92] and stochastic [131] global optimization. Early topics in deterministic global
optimization include convex optimization [94], Branch-and-Bound techniques restricted to par-
ticular classes of problems (e.g. concave minimization problems [53]) and some theoretical and
complexity-related studies [135, 97]. Stochastic algorithms based on adaptive random search
appeared between the 1970s and early 1980s [74, 77]. It is worth noting that the first deter-
ministic global optimization technique that was able to deal with generic nonconvex continuous
NLPs was Interval Optimization [48, 49]). Unfortunately, it often has slow convergence due to
the fact that Interval Arithmetic generally provides very wide intervals for the objective function
value.

Toward the end of the 1980s, many articles on deterministic global optimization started to appear in the literature. They were mostly concerned with iterative methods applied to particular classes of problems [89, 81, 54, 56, 134, 63, 82], but there were also theoretical studies [50, 55] and parallel implementation studies and reports [90, 42, 33]. In the area of stochastic optimization, some of the first algorithms employing “tunnelling” [145] began to appear in the same period.

Since the beginning of the 1990s, the optimization research community has witnessed an explosion of papers, books, algorithms, software packages and resources concerning deterministic and stochastic global optimization. In the early 1990s, most of the articles were still concerned with applications of global optimization or algorithms which perform well on particular classes of problems. It is significant that one of the first and most widely used book in global optimization [59], first published in 1990 and successively re-published in 1993 and 1996, does not even mention general nonconvex NLPs in the generic form (1.1). The same is true even for a “survey book” ([57]) which appeared in 1995.

The first method that was able to deal directly with the generic nonconvex NLPs in the form (1.1) was Ryoo and Sahinidis’ Branch-and-Reduce algorithm which appeared in a paper published in May 1995 [99, 100]. Shortly afterward, Floudas’ team published their first article on the αBB branch-and-bound method [13] which was then thoroughly explored and analysed in several subsequent papers by Adjiman [7, 6, 3, 4, 5, 8]. The first αBB variant that addressed problems in the form (1.1) appeared in 1997 [3].

One notable limitation of the αBB algorithm is that it relies on the functions being twice differentiable in the continuous variables. Since the inception of the αBB algorithm, a number of Branch-and-Select algorithms geared toward the most generic nonconvex MINLP formulation appeared in the literature, like Smith and Pantelides’ symbolic reformulation approach [114, 115, 116], Pistikopoulos’ Reduced Space Branch-and-Bound approach [32] (which only applies to continuous NLPs), Grossmann’s Branch-and-Contract algorithm [146] (which also only applies to continuous NLPs) and Barton’s Branch-and-Cut framework [61]. Within the Branch-and-Select strategy, we can also include modern Interval Analysis based global optimization methods [137, 84, 136].

Branch-and-Select is by no means the only available technique for deterministic global optimization, but it seems to be the method which least relies on the problem having a particular structure: this is an advantage insofar as one needs to implement software which solves prob-
lems in the general form (1.1). For example, d.c. optimization and the Extended Cutting Planes (ECP) method [143, 142] are structurally different from Branch-and-Select approaches; however, both make use of a problem which is already in a special form<fnref n="4" />. It must be said, however, that “special forms” are sometimes general enough to accommodate large classes of problems. For instance, the ECP method solves mixed-integer NLPs (MINLPs) in the following form: $\min\{f(z) | g(z) \le 0 \wedge Az \le a \wedge Bz = b \wedge z \in \mathbb{R}^n \times \mathbb{Z}^m\}$, where $f$ is a pseudo-convex function, $g$ is a set of pseudo-convex functions, $A, B$ are matrices and $a, b$ are vectors. Although not all functions are pseudo-convex, the latter do form quite a large class. It is the pseudo-convexity condition that allows the ECP method to derive a new valid cutting plane at each iteration.

In the stochastic global optimization field, recent advances include Differential Evolution [119], Adaptive Lagrange-Multiplier Methods [140], Simulated and Nested Annealing [96], Ant Colony Simulation [78, 25], Quantum Dynamics in complex biological evolution [52], Quantum Thermal Annealing [68], Ruin and Recreate Principle [105] and Tabu Search [23].

### 1.3.3 Two-phase global optimization algorithms

Most methods for global optimization are based on two phases [104]: a global search phase and a local search phase. The reason for this is that there are very efficient local optimization methods which can be employed in the local search phase, but these methods are not generally capable of determining a global optimum of a problem in its most generic form (1.1). The local phase is nearly always deterministic in nature whereas the global phase may be either deterministic or stochastic.

In the global phase, the problem is limited to a particular region of space where one can prove (or just hopes) that a local search will find the global optimum. This procedure is iterative in nature so that the search space can be explored exhaustively.

Local optimization of nonconvex problems is an NP-hard problem [93]. Because all iterative methods for global optimization rely on a local search phase, under the most general conditions global optimization of NLPs is an NP-hard problem. The only non-iterative method for global optimization, based on Gröbner bases [47], is also NP-hard [26].

This research focuses on deterministic algorithms for global optimization. In particular, the Branch-and-Select family of algorithms promises a good performance in terms of computational

<footnote n="4">⁴I.e., objective function and constraints must be pseudo-convex. This form is needed to guarantee global optimality.</footnote>
time, it can be applied to problems in very general form, and is deterministic in nature.

## 1.4 The branch-and-select strategy

Branch-and-Select algorithms include well-known techniques such as Branch-and-Cut and Branch-and-Bound, and are among the most effective methods for the deterministic solution of global optimization problems. They started out as divide-and-conquer type algorithms to solve combinatorial optimization problems like the Travelling Salesman Problem [9] but were very soon applied to continuous and mixed-integer nonlinear optimization [34, 59, 92]. In this section, we present a general theory of such algorithms (based on material from [133]), together with the necessary convergence proofs.

Branch-and-Select algorithms can be used to solve the widest possible class of optimization problems (1.1) to global optimality, even when objective function and constraint values are provided by black-box procedures. In fact, they can be designed to rely on no particular feature of the problem structure. However, in their basic form and without any acceleration devices (such as pre-processing steps or improved algorithmic steps), they tend to be very inefficient. Furthermore, their performance may depend very much on the formulation of the problem, meaning that a different algebraic form of the same equations might lead to quite different algorithmic performance.

Let $\Omega \subseteq \mathbb{R}^n$. A finite family of sets $\mathcal{S}$ is a *net* for $\Omega$ if it is pairwise disjoint and it covers $\Omega$, that is, $\forall s, t \in \mathcal{S} (s \cap t = \emptyset)$ and $\Omega \subseteq \bigcup_{s \in \mathcal{S}} s$. A net $\mathcal{S}'$ is a *refinement* of the net $\mathcal{S}$ if there are finitely many pairwise disjoint $s'_i \in \mathcal{S}'$ such that $s = \bigcup_i s'_i \in \mathcal{S}$ and $s \notin \mathcal{S}$. In other words, if $\mathcal{S}'$ is a refinement of $\mathcal{S}$, it has been obtained from $\mathcal{S}$ by finitely partitioning some set $s$ in $\mathcal{S}$ and then replacing $s$ by its partitions.

Let $\mathcal{S}_n$ be an infinite sequence of nets for $\Omega$ such that, for all $i \in \mathbb{N}$, $\mathcal{S}_i$ is a refinement of $\mathcal{S}_{i-1}$, and let $M_n$ be an infinite sequence of subsets of $\Omega$ such that $M_i \in \mathcal{S}_i$. $M_n$ is a *filter* for $\mathcal{S}_n$ if $\forall i \in \mathbb{N} (M_i \subseteq M_{i-1})$. Let $M_\infty = \bigcap_{i \in \mathbb{N}} M_i$ be the *limit* of the filter.

We will now present a general framework for Branch-and-Select algorithms that globally solve the generic problem $\min\{f(x) | x \in \Omega\}$. Let $\gamma \in \mathbb{R}$. Given any net $\mathcal{S}$ for $\Omega \cap \{x | f(x) < \gamma\}$, we consider a selection rule that determines:

1. a distinguished point $\omega(M)$ for every $M \in \mathcal{S}$,
2. a subfamily $\mathcal{R}$ of qualified members of $\mathcal{S}$,

3. a distinguished member $M^*(\mathcal{S})$ of $\mathcal{R}$,

such that, for all $x \in \bigcup_{s \in \mathcal{S} \setminus \mathcal{R}} s$, we have $f(x) \ge \gamma$. Basically, the selection rule rejects regions which cannot contain the global optimum by choosing the qualified members; for each of these members, it calculates a distinguished point (for example, by applying a local optimization procedure to the specified region) and a distinguished member for further net refinement. The prototype algorithm below also relies on a selection rule used to partition a region.

1. (Initialization) Start with a net $\mathcal{S}_1$ for $\Omega$, set $x_0 = \emptyset$ and let $\gamma_0$ be any upper bound for $f(\Omega)$. Set $\mathcal{P}_1 = \mathcal{S}_1$, $k = 1$, $\sigma_0 = \emptyset$.

2. (Evaluation) Let $\sigma_k = \{\omega(M) \mid M \in \mathcal{P}_k\}$ be the set of all distinguished points.

3. (Incumbent) Let $x_k$ be the point in $\{x_{k-1}\} \cup \sigma_k$ such that $\gamma_k = f(x_k)$ is lowest.

4. (Screening) Determine the family $\mathcal{R}_k$ of qualified members of $\mathcal{S}_k$ (in other words, reject the unqualified members, i.e. those that can be shown not to contain a solution better than $\gamma_k$).

5. (Termination) If $\mathcal{R}_k = \emptyset$, terminate. The problem is infeasible if $\gamma_k \ge \gamma_0$; otherwise $x_k$ is the global optimum.

6. (Selection) Select the distinguished member $M_k = M^*(\mathcal{S}_k) \in \mathcal{R}_k$ and partition $M_k$ according to a pre-specified branching rule. Let $\mathcal{P}_{k+1}$ be the partition of $M_k$. In $\mathcal{R}_k$, replace $M_k$ by $\mathcal{P}_{k+1}$, thus obtaining a new refinement net $\mathcal{S}_{k+1}$. Set $k \leftarrow k + 1$ and go back to Step 2.

A Branch-and-Select algorithm is convergent if $\gamma^* = \inf f(\Omega) = \lim_{k \to \infty} \gamma_k$. A selection rule is exact if:

1. the infimum objective function value of any region that remains qualified during the whole solution process is greater than or equal to the globally optimal objective function value, i.e.

$$ \forall M \in \bigcap_{k=1}^{\infty} \mathcal{R}_k \ (\inf f(\Omega \cap M) \ge \gamma^*) $$

2. the limit $M_\infty$ of any filter $M_k$ is such that $\inf f(\Omega \cap M_\infty) \ge \gamma^*$.
In this theoretical set-up, it is easy to show that a Branch-and-Select algorithm using an exact selection rule converges.

### 1.4.1 Theorem

*A Branch-and-Select algorithm using an exact selection rule converges.*

*Proof.* Suppose, to get a contradiction, that there is $x \in \Omega$ with $f(x) < \gamma^*$. Let $x \in M$ with $M \in \mathcal{R}_n$ for some $n \in \mathbb{N}$. Because of condition (1) above, $M$ cannot remain qualified forever; furthermore, unqualified regions may not, by hypothesis, include points with better objective function values than the current incumbent $\gamma_k$. Hence $M$ must necessarily be split at some iteration $n' > n$. So $x$ belongs to every $M_n$ in some filter $\{M_n\}$, thus $x \in \Omega \cap M_\infty$. By condition (2) above, $f(x) \ge \inf f(\Omega \cap M_\infty) \ge \gamma^*$. The result follows. □

It is worth noting that in this algorithmic framework does not provide a guarantee of a finite convergence. Consequently, most Branch-and-Select implementations make use of the concept of $\varepsilon$-optimality, rather than the usual definition of optimality. Recall that $x^* \in \Omega$ is a global optimum if, for all $x \in \Omega$, we have $f(x^*) \le f(x)$. Given $\varepsilon > 0$, $\bar{x} \in \Omega$ is $\varepsilon$-globally optimal if there exist bounds $m \le f(x^*) \le M$ such that $f(\bar{x}) \in [m, M]$ and $M - m < \varepsilon$. By employing this notion and finding converging lower and upper bounds sequences to the incumbent at each step of the algorithm, it is easier to ensure a finite termination with an $\varepsilon$-global optimum. For a theoretically proven finite termination, we would need some additional regularity assumptions (see for example [106, 11]); $\varepsilon$-optimality, however, is sufficient for most practical purposes.

### 1.4.1 Fathomning

Let $S_k$ be the net at iteration $k$. For each region $M \in S_k$, find the lower bound $l(M)$ of $f(M \cap \Omega)$ and define $\omega(M) = l(M)$. A set $M \in S_k$ is qualified if $l(M) \le \gamma_k$. The distinguished region is usually selected as the one with the lowest $l(M)$ (this may vary in some implementations).

The algorithm is accelerated if one also computes an upper bound $u(M)$ of $f(M \cap \Omega)$ and uses it to bound the problem from above by rejecting any $M$ for which $l(M)$ exceeds the best upper bound that has been encountered so far. In this case, we say that the rejected regions have been *fathomed* (see Fig. 1.1). This acceleration device has become part of most Branch-and-Bound algorithms in the literature. Usually, the upper bound $u(M)$ is computed by solving the problem to local optimality in the current region: this also represents a practical alternative to the implementation of the evaluation step (step (2) in the algorithm of Section 1.4), since we
can take the distinguished points $\omega(M)$ to be the local solutions $u(M)$ of the problem in the current regions. With a good numerical local solver, the accuracy of $\epsilon$-global optimum is likely to be better than if we simply use $l(M)$ to define the distinguished points.

<figure><img src="image_3.png" /><figcaption>Figure 1.1: Fathoming via upper bound computation.</figcaption></figure>

The convergence of the Branch-and-Bound algorithm is ensured if every filter $\{M_k | k \in K\}$ contains an infinite nested sequence $\{M_k | k \in K_1\}$ such that:

$$ \lim_{\substack{k \to \infty \\ k \in K_1}} l(M_k) = \gamma_k. \qquad (1.2) $$

To establish this, we will show that under such conditions the selection procedure is exact. Let $\{M_k | k \in K\}$ be any filter and $M_\infty$ its limit. Because of equation (1.2), $\inf f(\Omega \cap M_k) \ge l(M_k) \to \gamma^*$, hence $\inf f(\Omega \cap M_\infty) \ge \gamma^*$. Furthermore, if $M \in \mathcal{R}_k$ for all $k$ then $\inf f(\Omega \cap M) \ge l(M) \ge l(M_k) \to \gamma^*$ as $k \to \infty$, i.e. $\inf f(\Omega \cap M) \ge \gamma^*$. Thus the selection procedure is exact and, by theorem 1.4.1, the Branch-and-Bound algorithm converges.

## 1.5 Mathematical formulation and convex relaxation for non-convex NLPs

This thesis is concerned with techniques for the determination of global optima of nonconvex NLPs of the general form (1.1). We are particularly interested in improving the performance of branch-and-bound algorithms of the type outlined in Section 1.4. As indicated by the literature review presented in Section 1.3.2 this has been an area of much active research worldwide, especially over the past decade.
Albeit themselves relatively new, branch-and-bound algorithms for global optimization have many similarities to those which have been routinely used for the solution of mixed integer linear programming (MILP) problems for more than four decades – indeed, both classes of algorithms<fnref n="5" /> conform to the general Branch-and-Select framework described in Section 1.4. It is, therefore, instructive to draw analogies between these two closely related types of algorithm when trying to identify directions for fruitful research. In particular, the experience of applying branch-and-bound algorithms to MILPs points out the extreme importance of mathematical formulation. More specifically, the same engineering problem (e.g. process scheduling) may often be expressed mathematically in two or more different ways. Although all of these formulations are mathematically completely equivalent, the computational effort required for their solution is very different, often varying by several orders of magnitude. This indicates that the issue of formulation is likely to be central to devising efficient global optimization solution methods, perhaps much more so than, for instance, algorithmic details such as rules for selecting branching variables.

Better mathematical formulations for specific NLPs may sometimes be devised on the basis of insight and intuition. However, it is clear that methods which could *automatically* reformulate wide classes of nonconvex NLPs to a better behaved form would be highly desirable.

A topic that is closely related to that of mathematical formulation is that of the convex relaxation of a given nonconvex NLP. This is essential for all sBB algorithms as a means for establishing a lower bound on the objective function in any given region. The analogous convex relaxation for the case of Branch-and-Bound algorithms for MILPs is the “relaxed” linear programming (LP) problem, i.e. one obtained from the MILP by dropping the integrality requirement on (some of the) integer variables. This convex relaxation is universally used by MILP solution algorithms, and therefore, a “good” mathematical formulation for an MILP is usually considered as one which exhibits the smallest gap between the optimal objective function of the original MILP and that of its LP relaxation.

Interestingly, a significant difference between branch-and-bound algorithms for MILPs and those for nonconvex NLPs is that, in the latter (sBB) case, there is no single convex relaxation of a nonconvex NLP. Consequently, albeit still closely related, the issues of formulation and convex relaxation are no longer identical. For example, it is one important issue how to formulate an engineering optimization problem so that it exhibits the minimum number of nonconvexities –

<footnote n="5">⁵In this thesis, we shall follow the established practice of referring to branch-and-bound algorithms for global optimization as *spatial branch-and-bound (sBB) algorithms* in order to distinguish them from those traditionally used for MILPs.</footnote>
and there may be several ways in which to address this requirement, each leading to a different NLP. Now, given any particular nonconvex NLP formulation, there may then also be more than one way in which to construct its convex relaxation, and selecting the best one is quite another issue.

## 1.6 Outline of this thesis

In view of the discussion in the previous section, it is clear that both the mathematical formulation and the convex relaxation of nonconvex NLPs deserve close attention. Chapter 2 of this thesis presents a review of the related literature, attempting to come up with a rational classification of the quite diverse approaches and techniques that have been reported to date.

Chapter 3 presents a novel automatic reformulation method for nonconvex NLPs that include linear constraints and bilinear terms. The latter are well-known sources of nonconvexity in a wide range of models of engineering systems, for example whenever an extensive variable (e.g. flowrate) multiplies an intensive one (e.g. composition, enthalpy etc.) within a conservation law. The method presented automatically reformulates the nonconvex NLP to one that (a) involves fewer bilinear terms and more linear constraints than the original, and (b) has a tighter convex relaxation when the remaining bilinear terms are convexified using standard techniques.

The objective function and constraints in nonconvex NLPs usually involve a variety of types of sub-expressions (e.g. bilinear, trilinear and linear fractional terms) and other mathematical functions. Relatively tight convex relaxations exist for many of these categories of expressions. One notable exception is that of terms involving monomials of odd degree (i.e. expressions of the form $x^n$ where $n$ is odd) when the range of the variable $x$ includes zero. These occur often (e.g. as cubic or quintic expressions) in practical applications. Chapter 4 presents a novel method for constructing the convex envelope of such monomials, as well as a tight linear relaxation of this envelope.

The automatic reformulation and convex relaxation techniques outlined above fit within the framework of a general sBB algorithm. The form of this algorithm as discussed in the literature is reviewed in Chapter 5. The software implementation of such an algorithm places special demands since, beyond the usual numerical information required by most optimization codes, a substantial amount of symbolic information must also be made available. The architecture of *ooOPS*, a general object-oriented software framework that attempts to fulfill these, is also presented in Chapter 5.
Finally, Chapter 6 summarizes the main points of the work presented and draws some general conclusions.
# Chapter 2

## Overview of reformulation techniques in optimization

A key theme of this thesis is the issue of problem reformulations that can be carried out in an automatic manner via symbolic and numerical procedures. Such reformulations can be used for various tasks, both prior to the actual problem solution and during the solution process. Many types of automatic reformulations have appeared in the literature; however, to the best of our knowledge, no attempts have been made, to date, to review them and classify them in a systematic manner. This chapter, therefore, attempts to present such a literature review.

A reformulation of an optimization problem $P$ is a problem $P'$ which shares some mathematical properties with $P$. A reformulation may be useful as part of an algorithmic procedure to solve $P$, or if it offers qualitative or quantitative insights regarding the properties of $P$.

A reformulation $P'$ is *exact* if the global solution $x \in \mathbb{R}^n$ of $P$ can be directly inferred from the global solution $x'$ of $P$ ("directly inferred" meaning that $x$ can be computed from $x'$ in linear time depending on $n$). An exact reformulation is *convenient* if finding $x'$ requires less computational resources than finding $x$.

A useful reformulation which is not exact is usually called a *relaxation* of the problem. For example, many algorithms for the solution of discrete optimization problems involve a continuous relaxation of the problem (i.e. the discrete variables are reformulated to continuous variables). Most Branch-and-Bound techniques for nonconvex optimization problems use a convex relaxation of the problem to compute the lower bound at each iteration.
Reformulations can be constructed by using algebraic-symbolic manipulation of the equa-
tions in the original problem, or by using numerical computations. In both cases, it is necessary
to prove theoretically, or at least offer strong evidence, that the reformulation is indeed useful.
If a reformulation is sometimes useful but it cannot always be decided a priori whether it will
be useful or not, and to what extent, then it is called heuristic.

In the rest of this chapter we shall give a detailed account of the most useful reformulation
techniques used in conjunction with deterministic methods for global optimization.

## 2.1 Reformulations to standard forms

Closely associated with the idea of a reformulation is the concept of a standard form for an
optimization problem. Solution algorithms often require the optimization problem to be in a
pre-specified form called the *standard form* with respect to that algorithm. Most reformulations
are meant to unearth some structural mathematical property of *P* which was not evident before,
or to transform *P* into a standard form.

In this section we refer to optimization problems in the following form:

$$
\min_{x \in \Omega} f(x) \tag{2.1}
$$

where $\Omega$ is the feasible region.

Notation-wise, by $Ax = b$ we mean a system of $m$ linear equations in $n$ variables, where $A$ is an $m \times n$ matrix, $x \in \mathbb{R}^n$ (unless specified otherwise) and $b \in \mathbb{R}^m$. Likewise, we express a system of linear inequalities as $Ax \le b$. By $x^L \le x \le x^U$ we mean a set of ranges on the variables: $\forall i \le n(x_i^L \le x_i \le x_i^U)$, and $x^L, x^U \in \mathbb{R}^n$.

2.1.1 Box-constrained problems

An optimization problem (2.1) is box-constrained if the only constraints of the problem are the variable ranges (i.e. $\Omega = \{x \mid x^L \le x \le x^U\}$). In other words, no equation or inequality constraints are present (sometimes such problems are called *unconstrained problems*, though in fact a truly unconstrained optimization problem would lack variable bounds as well). This is a very well-studied and interesting class of problems [35, 14, 38], and much effort has gone into reformulations of other types of problems to this type. Moreover, most stochastic optimization
methods are originally devised for box-constrained problems, and are then extended to deal with constraints.

### 2.1.1.1 Penalty and barrier functions

The most widely used reformulation for eliminating equation or inequality constraints from constrained optimization problems is to employ *penalty* or *barrier* functions. This embeds the constraints into the objective function *f* itself, so that *f* attains high values at infeasible points; hence a minimization of the objective function automatically excludes the infeasible points.

Consider the optimization problem (2.1) where:

$$ \Omega = \{ x \in \mathbb{R}^n \mid \forall i \le m (h_i(x) = 0) \wedge \forall i \le m' (g_i(x) \le 0) \wedge x^L \le x \le x^U \}. \quad (2.2) $$

Let $\delta_i, \varepsilon_i: \mathbb{R} \to \{0, \infty\}$ be boolean functions that specify whether a constraint is satisfied: for all $i \le m$ define $\delta_i(h_i(x)) = 0$ if $h_i(x) = 0$ and $\delta_i(h_i(x)) = \infty$ if $h_i(x) \ne 0$. Likewise for inequality constraints, for all $i \le m'$ define $\varepsilon_i(g_i(x)) = 0$ if $g_i(x) \le 0$ and $\varepsilon_i(g_i(x)) = \infty$ otherwise. Now, reformulating the problem to:

$$ \min_{x^L \le x \le x^U} F(x) = f(x) + \sum_{i=1}^{m} \delta_i(h_i(x)) + \sum_{i=1}^{m'} \varepsilon_i(g_i(x)) \quad (2.3) $$

is an exact reformulation of the original problem to a box-constrained form. If *x* is feasible in the original problem, then the objective function *F* of the reformulation reduces to *f*(*x*); if it is infeasible, *F* becomes *∞*.

The main disadvantage of this reformulation is that it is highly non-smooth. Furthermore, the introduction of infinity in the definition of $\delta_i, \varepsilon_i$ is questionable from a numerical point of view. The latter point can be disposed of in case the function *f* is Lipschitz (i.e. there is a real constant $M > 0$ such that for all $x, y \in \Omega$ we have $|f(x) - f(y)| \le M||x - y||$), in which case it is possible to find a global upper bound *L* of *f* on *$\Omega$, and $\infty$ can be replaced by *L*. The non-smoothness is a more delicate question. Unless one devises an entirely symbolic method for solving (2.3), one has to replace the functions $\delta_i, \varepsilon_i$ with smoother functions. Unfortunately, this means that the reformulation may no longer be exact. In [14, 35] a host of different penalty functions are introduced, and their properties are discussed. The simplest “realistic” penalty function involves a reformulated objective function $F(x) = f(x) + \mu (\sum_{i=1}^{m} |h_i(x)| + \sum_{i=1}^{m'} \max(0, g_i(x)))$. In this case the main drawback is the determination of the parameter $\mu$. There is a value of $\mu$ such that this reformulation is exact, but it cannot be known *a priori*. Furthermore, absolute
value and pointwise maximum are not smooth functions everywhere. Other penalty and barrier functions involve logarithms, inverses, powers and so on.

### 2.1.1.2 Lagrangian and Lagrange coefficients

The *Lagrangian* function of an optimization problem in form 2.1, where $\Omega$ is defined as in equation (2.2), is defined as $L(x, \lambda, \mu) = f(x) + \sum_{i=1}^{m} \lambda_i h_i(x) + \sum_{i=1}^{m'} \mu_i g_i(x)$ with the requirement that $\mu_i \ge 0$ for all $i \le m'$. The parameters $\lambda_i, \mu_i$ are called the *Lagrange multipliers*. The Lagrangian is often used in deriving theoretical conditions of optimality in methods of local optimization (Karush-Kuhn-Tucker necessity and sufficiency conditions, see [38, 85]).

In fact, the theory of Lagrange multipliers is also linked to duality theory for general optimization problems. It can be shown that, for a convex problem $P$, its Lagrangian dual:

$$ \max_{\lambda; \mu > 0} \min_{x^L \le x \le x^U} L(x, \lambda, \mu) \qquad (2.4) $$

has exactly the same solution as $P$. For a nonconvex problem it can be guaranteed only that the solution of the dual problem is a lower bound to the solution of the original problem (2.2). The difference between the value of the dual objective function and the original objective function is called the *duality gap*. In view of the fact that the solution to the dual problem (2.4) is often used as a lower bounding technique in Branch-and-Bound algorithms, it is important to find methods to reduce the duality gap.

There are various such methods: in [17] a suitable partition of the variable ranges is employed. The dual problem is then solved over each subset of the partition. In [70, 69] the Lagrangian $L(x, \lambda, \mu)$ is reformulated to:

$$ L_p(x, \lambda, \mu) = f^p(x) + \sum_{i=1}^{m} \lambda_i h_i^p(x) + \sum_{i=1}^{m'} \mu_i g_i^p(x) $$

where $p \ge 1$ is an exponent. The reformulated Lagrangian has local convexity properties in correspondence of local solutions of the original problem. More precisely, under some further regularity conditions, the Hessian of $L_p$ is positive definite at $x^*$ (a locally optimal solution of the original problem (2.2)) for large enough $p$. This is a sufficient condition for $L_p$ to be locally convex in a neighbourhood of $x^*$, and hence for $L_p$ to be a valid lower bounding function for the original problem.
## 2.1.2 Separable problems

Separable programming problems are useful because the objective function can be expressed as a sum of functions of one variable only. Thus, each term in the sum is independent of the other terms. Separable problems have a special structure which offers wide scope for decomposition strategies. If the problem is separable and box-constrained (i.e. $\Omega = \{x | x^L \le x \le x^U\}$), then it can easily be solved to global optimality via interval analysis (see Section 2.1.2.2). If linear or nonlinear constraints are present, interval analysis offers a very fast way to calculate lower bounds on the objective function: it is no surprise that the first Branch-and-Bound approaches to global optimization were restricted to separable problems [34, 117, 15].

Let $x = (x_1, \dots, x_n) \in \Omega$. A function $f: \mathbb{R}^n \to \mathbb{R}$ is separable if and only if there are $f_1, \dots, f_n: \mathbb{R} \to \mathbb{R}$ such that for all $x \in \Omega$ we have $f(x) = \sum_{i=1}^n f_i(x_i)$.

If the objective function $f$ of a problem (2.1) is separable, then the problem is separable. The definition of the feasible region $\Omega$ varies. Usually [14, 35, 79], $\Omega$ is defined by a system of separable inequalities $\sum_{j=1}^n g_{ij}(x_j) \le b_i, \forall i \le m$, where $g_{ij}: \mathbb{R} \to \mathbb{R}$ for each $i, j$. Some authors [27, 58] require $\Omega$ to be a polytope.

### 2.1.2.1 Separation of bilinear forms

In the procedure of Section 2.1.8.1, the basic idea is the separation of a bilinear form $xy$. This idea rests on the relationship $(x+y)^2 = x^2+2xy+y^2 \Rightarrow xy = \frac{1}{2}(w^2-x^2-y^2)$ and $w = x+y$. It was already known in the 18th century that the quadratic form $\sum_{i,j=1}^n a_{ij}x_ix_j$ could always be reduced to a sum or difference of squares:

$$y_1^2 + \dots + y_r^2 - y_{r+1}^2 - \dots - y_{r+t}^2 \quad (2.5)$$

via a real linear transformation $x_i = \sum_j b_{ij}y_j$ (for all $i \le n$) having non-zero determinant. Cauchy, Sylvester and Jacobi all worked on this problem. If $t = 0$, then the form is called *positive definite*; if $t > 0$ then the form is called *semidefinite* (these terms were introduced by Gauss in his *Disquisitiones arithmeticae* [62]). Their work is also relevant to the theory of convex, concave and d.c. functions, since a positive quadratic term is a convex function and a negative quadratic term is a concave function.

The problem of reducing quadratic forms is tightly linked to the problem of diagonalizing a square matrix. Let $x = (x_1, \dots, x_n)$, $s = (s_1, \dots, s_n)$ and $A = (a_{ij})$ be an $n \times n$ matrix. The expression $x^T A s = \sum_{i,j} a_{ij}x_i s_j$ is called a *generalized bilinear form*. We are specially
interested in the case where $x = s$. To reduce a generalized bilinear form to a semidefinite quadratic form, we have to find a non-singular transformation $x = Py$, where $P$ is an $n \times n$ matrix and $y = (y_1, \dots, y_n)$ such that $x^T A x = y^T P^T A P y$ and $P^T A P$ is diagonal. Constructing such a matrix $P$ is a classical problem in matrix theory (see for example [76], p. 322).

Reformulating generalized bilinear forms to semidefinite quadratic forms is usually convenient. It is an exact reformulation and it often gives rise to better convex relaxations (see Section 2.3).

### 2.1.2.2 Global solution of separable box-constrained problems

Finding the global optimum of the separable box-constrained problem

$$ \min_{x^L \le x \le x^U} \sum_{i=1}^{n} f_i(x_i) $$

reduces to finding the global optimum of each of the one-dimensional problems

$$ \min_{x_i^L \le x_i \le x_i^U} f_i(x_i), $$

and can therefore be totally decomposed. Interval arithmetic provides an extremely fast and effective way to find the bounds on a univariate function $f_i(x_i)$ given the range $[x_i^L, x_i^U]$ of $x_i$.

Let $\mathcal{N} = \{N_1, \dots, N_m\}$ be a partition of $\{1, \dots, n\}$, and for each $i \le m$ let $x^{[i]} = \{x_j | j \in N_i\}$. The function $f: \mathbb{R}^n \to \mathbb{R}$ is semi-separable if and only if for each $i \le m$ there exist functions $f_i: \mathbb{R}^{|N_i|} \to \mathbb{R}$ such that $f(x) = \sum_{i=1}^m f_i(x^{[i]})$. The method described in the example below also works when applied to problems with semi-separable objective functions.

#### 2.1.1 Example

Find the global minima of the problem:

$$ \begin{aligned} \min_{x_1, x_2, x_3, x_4} & f(x) = x_1 x_2 + x_3^4 - 10 x_3^3 + \cos(e^{x_4}) \\ & -1 \le x_1, x_2 \le 1 \\ & -5 \le x_3 \le 10 \\ & 0 \le x_4 \le 2. \end{aligned} $$

Notice $f$ is semi-separable, with $f_1(x_1, x_2) = x_1 x_2$, $f_2(x_3) = x_3^4 - 10 x_3^3$ and $f_3(x_4) = \cos(e^{x_4})$.

Interval analysis on the product $x_1 x_2$ shows that the global minima of $f_1(x_1, x_2)$ in the specified
range are attained at $(x_1, x_2) = (-1, 1)$ and $(x_1, x_2) = (1, -1)$. The global minimum of $f_2(x_3)$ can be found either algebraically, as $f'_2(x_3) = 4x_3^3 - 30x_3^2 = 2x_3^2(2x_3 - 15)$ implies a global minimum at $x_3 = \frac{15}{2}$, or numerically, by applying Newton's method in one dimension for finding the roots of the derivative. Finally, the global minimum of $f_3(x_4) = \cos(e^{x_4})$ can be found by setting $\cos(e^{x_4}) = -1$, which implies $e^{x_4} = (2k + 1)\pi$ for all $k \in \mathbb{Z}$, i.e. $x_4 = \log((2k + 1)\pi)$. Given the range of $x_4$, we have to impose $k = 0$ and so $x_4 = \log(\pi) = 1.14473$. Thus the global solutions to the problem are at $(x_1, x_2, x_3, x_4) = (\pm 1, \mp 1, \frac{15}{2}, \log(\pi))$.

### 2.1.3 Linear problems

A linear optimization problem is such that both the objective function and the constraints are linear functions in the problem variables. A very efficient global solution method (called the simplex method) for linear continuous optimization problems was proposed in the 1940s, by Dantzig [27], and has been further refined since then. Nowadays there are many efficient software codes to solve large-scale linear problems.

Unfortunately, linear reformulations are very rarely exact. In fact, most linear reformulations are linear relaxations and are used within more complex methods for the solution of nonlinear problems (like e.g. Branch-and-Bound). We will analyze them in Section 2.4.

#### 2.1.3.1 Reformulating quadratic binary problems to linear binary problems

Although this reformulation only applies to a very special class of optimization problems, it is one of the very few exact linear<fnref n="1" /> reformulations [16, 83]. Any unconstrained quadratic binary problem $\max_{x \in \{0,1\}^n} x^T Q x$ can be reformulated exactly to:

$$ 
\begin{array}{ll}
\max & q^T y \\
\text{s.t.} & \forall i, j \le n \ (y_{ij} \le x_i) \\
& \forall i, j \le n \ (y_{ij} \le x_j) \\
& \forall i, j \le n \ (y_{ij} \ge x_i + x_j - 1) \\
& x \in \{0, 1\}^n, y \in \{0, 1\}^{n^2},
\end{array}
$$

<footnote n="1">¹Since the considered problems are binary, calling this a linear relaxation is a slight abuse of notation. What we mean is that the objective functions and the constraints are linear functions of the decision variables.</footnote>
where $q$ is a vector consisting of the entries of the $n \times n$ matrix $Q$ ordered by column. Such a reformulation belongs to a class of reformulations called *lifting*s because of the fact that they “lift” the geometry of the problem into a higher dimensional Euclidean space (i.e. they add new variables to the problem). Although usually a problem with more variables is more difficult to solve, *lifting*s can nonetheless be useful.

### 2.1.4 Convex problems

An optimization problem is convex if both the objective function $f$ and the feasible region $\Omega$ are convex. This is an interesting class of problems as it is possible to show that any local solution of a convex problem is also a global solution [85, 98]. As in the linear case, convex reformulations are rarely exact. However, convex relaxations (which will be analysed in more detail in Section 2.3) are used within more complex procedures for the solution of nonlinear problems. Their global minimality property makes it possible to compute lower bounds of the objective function of nonconvex problems in given subregions of the feasible region (this notion is used in most spatial Branch-and-Bound algorithms for the solution of NLPs and MINLPs [116, 4]).

One exact convex reformulation can be obtained when the objective function is a positive definite bilinear form (see Section 2.1.2.1). If $t = 0$ in equation (2.5), then the reformulation is exact and it makes the function completely convex.

Sometimes a nonconvex function can be reformulated exactly to a convex function by a nonlinear change of variables [98]. The function $f(x) = ax_1 \cdots x_n$ is nonconvex. For all $i \le n$ let $X_i = \log x_i$; then $f(X) = ae^{X_1+\cdots+X_n}$, which is convex. This, with a suitable adaptation, also applies if the variables are discrete [95]. The following, for example, is a mixed-integer reformulation of the function $ax_1^{r_1} \cdots x_n^{r_n}$ where $x_i$ are discrete variables for all $i \le n$:

$$ \left\{ \begin{array}{l} ax_1^{r_1} \cdots x_n^{r_n} \\ a > 0, \forall i \le n \ (r_i \in \mathbb{R}) \\ \forall i \le n \ (x_i \text{ discrete}) \end{array} \right. \Leftrightarrow \left\{ \begin{array}{l} ae^{r_1 X_1 + \cdots + r_n X_n} \\ X_i = \log d_{i1} + \sum_{j=1}^{n_i-1} \beta_{ij} (\log d_{i,j+1} - \log d_{i1}) \\ \sum_{j=1}^{n_i-1} \beta_{ij} \le 1 \\ X_i \in \mathbb{R}, \beta_{ij} \in \{0, 1\}, \end{array} \right. (*) $$

where the discrete variable $x_i$ can take values in $\{d_{i1}, \dots, d_{in_i}\}$. For an explanation of the equality constraints ($*$), see Section 2.1.5.1.

In [95], another type of reformulation is presented, that deals with terms of the form $-ax^p y^q$ with $a, p, q > 0$. If $p + q \le 1$ the term is convex, therefore we only need to apply the reformulation when $p + q > 1$. If we let $x = X^{\frac{1}{p+q}}$ and $y = Y^{\frac{1}{p+q}}$, the original term is transformed in
$-aX^{\frac{p}{p+q}}Y^{\frac{q}{p+q}}$ which is convex as the sum of the exponents is obviously less than or equal to 1. If $x, y$ were discrete variables taking values in $\{d_1, \dots, d_n\}$ and $\{c_1, \dots, c_m\}$ respectively, then the following integrality constraints would also be necessary:

$$
\begin{align*}
X &= d_1^{p+q} + \sum_{i=1}^{n-1} \beta_i (d_{i+1}^{p+q} - d_1^{p+q}) \\
Y &= c_1^{p+q} + \sum_{i=1}^{n-1} \alpha_i (c_{i+1}^{p+q} - c_1^{p+q}) \\
\sum_{i=1}^{n-1} \beta_i &\le 1 \\
\sum_{i=1}^{m-1} \alpha_i &\le 1 \\
\alpha_i, \beta_i &\in \{0, 1\}.
\end{align*}
$$

2.1.5 Binary problems

A binary optimization problem is such that the feasible region Ω is a subset of Bⁿ, where
B = {0, 1} (i.e. the problem variables can take only values 0 and 1). The solution space of
these problems is well suited to Branch-and-Select type searches. Each node of the search tree
has exactly two subnodes [9]. From a theoretical point of view, it is important because it is
the simplest type of discrete problem, and theorems can be stated and proved more simply if
variables can only take two values (see e.g. [109], p. 1277, Theorem 1).

2.1.5.1 Reformulating discrete problems to binary problems

A discrete optimization problem is such that the problem variables can take only discrete values.
It is possible to reformulate discrete problems to binary problems exactly [27, 144, 110]. Let
$v$ be a discrete problem variables taking values in the set $\{v_1, \dots, v_d\}$. By introducing $d$ new
binary variables $\beta_1 \dots \beta_d$, we can replace $v$ by $\sum_{i=1}^d v_i \beta_i$ and add a linear constraint $\sum_{i=1}^d \beta_i =
1$ to the definition of the feasible region $\Omega$. This is an exact reformulation: the linear constraint
ensures that exactly one binary variable $\beta_i$ takes the value of 1, so that the expression for $v$ takes
exactly one value in $\{v_1, \dots, v_d\}$. This reformulation involves only linear functions, so in this
sense it does not add significant complexity to the problem. This reformulation is a lifting.
A slightly better reformulation (in the sense that it requires one less binary variable) is the following [95]:

$$v = v_1 + \sum_{i=1}^{n-1} \beta_i (d_{i+1} - d_1) \wedge \sum_{i=1}^{n-1} \beta_i \le 1,$$

where $\beta_i \in \{0, 1\}$ for any $i \le n - 1$.

### 2.1.5.2 Reformulating binary problems to continuous problems

Sometimes it is more convenient to approach a discrete problem with continuous methods. To this end, the following quadratic exact reformulation can be used: substitute each binary variable $y \in \{0, 1\}$ with a continuous variable $\bar{y} \in [0, 1]$ and add the equality constraint $\bar{y} = \bar{y}^2$ to the formulation of the problem. This constraint, called an *integrality enforcing constraint*, is equivalent to the equation $\bar{y}(\bar{y} - 1) = 0$ which has solutions 0 and 1. Hence, any feasible solution to the continuous problem satisfying the integrality enforcing constraint will be such that $\bar{y} = 0$ or $\bar{y} = 1$.

This reasoning can be generalized to any discrete variable $y \in S = \{a_i | i \in \mathbb{N}\}$. Thus, we can relax $y$ to a continuous variable $\bar{y} \in [\min S, \max S]$ and add an equality constraint $g(\bar{y}) = 0$ having $S$ as solution set. If $S$ is finite, a possible function is $g(\bar{y}) = \prod_{i=1}^{|S|} (\bar{y} - a_i)$.

It is worth mentioning here that when trying to incorporate this relaxation within a Branch-and-Bound algorithm which requires convex relaxations to find the lower bounds of the objective function in each region, one needs to keep in mind that the function $g(\bar{y})$ will be replaced by its convex and concave relaxations, so $g$ should be chosen so that the “convexity gap” between $g$ and its relaxations is minimal (see also Section 2.3 about convex relaxations). In particular, Smith showed ([114], p. 209-210) that the simplest linear relaxation of the quadratic integrality enforcing constraint $\bar{y} = \bar{y}^2$ just reduces to $\bar{y} \in [0, 1]$ and hence no integrality is enforced on $\bar{y}$ when solving the relaxation.

Numerically, integrality enforcing constraints like $y = y^2$ are usually problematic for most nonlinear local solver codes. Thus, the idea of reformulating a large-scale binary problem to a continuous problem having high numbers of these constraints is not a workable one. However, such constraints are sometimes useful in MINLPs having few integer variables.
### 2.1.6 Concave problems

A concave optimization problem in form 2.1 is such that the objective function $f$ is concave and the feasible region $\Omega$ is convex. Concave optimization is, in a certain sense, the simplest case of a complicated optimization problem: it is multi-extremal (i.e. in general, it has many local minima) so that the techniques of local optimization are not sufficient to find the global optimum. Its complexity is high enough to allow many different optimization problem families to be reformulated exactly to a concave problem. However, its formulation is simple enough that efficient algorithms for its global solution can be designed [58, 18].

#### 2.1.6.1 Reformulating binary problems to concave problems

Let $C$ be a convex set (which may be defined by a set of constraints), $B = \{0, 1\}$, $\bar{B} = [0, 1]$ and consider a binary problem in the form $\min\{f(x) | x \in C \cap B^n\}$, where $f$ is Lipschitz and twice continuously differentiable on $\bar{B}^n$. Then there exists $\mu_0 \in \mathbb{R}$ such that, for all $\mu > \mu_0$, the binary problem above can be reformulated exactly to:

$$ \min_{x \in C \cap \bar{B}^n} F(x) = f(x) + \mu \sum_{i=1}^{n} x_i(1 - x_i), $$

and $F(x)$ is concave on $\bar{B}^n$. For the proof of this theorem see [58], p. 15. One convenient feature of this exact reformulation is that it adds no new variables to the problem. Furthermore, it succeeds in relaxing the binary variables to continuous whilst still keeping the reformulation exact.

#### 2.1.6.2 Reformulating bilinear problems to concave problems

The bilinear programming problem is stated in general terms as:

$$ \min_{x \in X, y \in Y} f(x, y) = px + xQy + qy \qquad (2.6) $$

where $X, Y$ are convex polyhedral sets (defined by sets of linear constraints). Assume $Y$ has at least one vertex and for every $x \in X$ the problem $\min_{y \in Y} f(x, y)$ has a solution. Then problem (2.6) can be reformulated exactly to a concave minimization problem with piecewise linear objective function and linear constraints [58]. In particular, let $V(Y)$ be the set of vertices of $Y$.
Since the solution of a linear problem is attained at least at one point in V(Y), we have

$$
\begin{align*}
\min_{x \in X, y \in Y} f(x, y) &= \min_{x \in X} \min_{y \in V(Y)} f(x, y) = \\
&= \min_{x \in X} f(x),
\end{align*}
$$

where $f(x) = \min_{y \in V(Y)} f(x, y)$. The set $V(Y)$ is finite and for each $y \in V(Y)$, $f(x, y)$ is a linear function of x. Thus $f(x)$ is the pointwise minimum of a finite family of linear functions, and hence is concave and piecewise linear.

2.1.6.3 Reformulating complementarity problems to concave problems

Complementarity problems are feasibility problems rather than optimization problems. The problem is to find $x \in \Omega \subseteq \mathbb{R}^n$ such that for all $i \le m$ we have $g_i(x) \ge 0$, $h_i(x) \ge 0$ and $g_i(x)h_i(x) = 0$, where $g_i, h_i: \mathbb{R}^n \to \mathbb{R}$. Complementarity problems arise in the study of other optimization problems, in the analysis of the computation of equilibria in fixed-point problems (theory of games) and in modelling the logical notion of disjunction (which concerns the modelling of constraints in certain combinatorial problems where a constraint is to be enforced only if a certain binary variable is true) [86].

Assume $\Omega$ is convex and $g_i, h_i$ are concave functions for each $i \le m$, and that the complementarity problem has a solution $x^*$. In this case, the complementarity problem can be reformulated exactly to the following concave optimization problem:

$$
\left.
\begin{array}{l}
\min_{x \in \Omega} f(x) = \sum_{i=1}^{m} \min\{g_i(x), h_i(x)\} \\
\text{s.t.} \quad \forall i \le m (g_i(x) \ge 0, h_i(x) \ge 0)
\end{array}
\right\}
$$

For the proof of this statement, see [58], p. 24.

2.1.6.4 Reformulating max-min problems to concave problems

A max-min optimization problem is stated as:

$$
\left.
\begin{array}{l}
\max_{x \ge 0} \min_{y \ge 0} ax + by \\
\text{s.t.} \quad Ax + By \le c
\end{array}
\right\}
$$

where $x, a \in \mathbb{R}^n$, $y, b \in \mathbb{R}^m$, $A$ is an $s \times n$ matrix, $B$ is an $s \times m$ matrix and $c \in \mathbb{R}^s$. This problem
can be reformulated exactly to a concave optimization problem $- \min_{x \in P} (-f(x) - ax)$, where:

$$
\begin{align*}
P &= \{x \ge 0 \mid \exists y \le 0 (Ax + By \le c)\} \\
f(x) &= \min\{by \mid By \le c - Ax \wedge y \ge 0\}.
\end{align*}
$$
Note that $f(x)$ is a convex piecewise linear function on $P$, so that $-f(x) - ax$ is concave and piecewise linear.

### 2.1.7 D.c. problems

An optimization problem is d.c. if the objective function is d.c. and $\Omega$ is a d.c. set (see Section 1.1) [132, 121, 120]; in fact in [121], the objective function is only required to be upper semicontinuous (i.e. for all $x$ in the domain of $f$ we have $f(x) = \lim_{y \to x} \sup f(y)$; the definition of lower semicontinuity is similar, with sup replaced by inf, and both semicontinuities at one point imply ordinary continuity at that point). In [59], the feasible region is defined as $C \cap M$, where $M = \{x \mid \forall i \le m (g_i(x) \le 0 \land g_i \text{ is d.c.})\}$.

D.c. programming problems have two fundamental properties. The first is that the space of all d.c. functions is dense in the space of all continuous functions. This implies that any continuous optimization problem can be approximated as closely as desired, in the uniform convergence topology, by a d.c. optimization problem [132, 59]. The second property is that it is possible to give explicit necessary and sufficient global optimality conditions for certain types of d.c. problems [132, 121]. Some formulations of these global optimality conditions [120] also exhibit a very useful algorithmic property: if at a feasible point $x$ the optimality conditions do not hold, then the optimality conditions themselves can be used to construct a better point $x'$.

#### 2.1.7.1 Reformulating continuous functions to d.c. functions

Each twice differentiable function $f: \mathbb{R}^n \to \mathbb{R}$ is d.c., since $g(x) = f(x) + \rho x^T x$ is a convex function for all sufficiently large $\rho$. However, since finding a good estimate for $\rho$ is itself a difficult problem [7], there is to date no automatic efficient procedure to reformulate any given continuous problem to a d.c. problem. However, some progress in this field has been made for the case of separable functions [67]. Assuming $f(x) = \sum_{i=1}^n f_i(x_i)$ and each $f_i$ is continuous but neither concave nor convex, then it is possible to reformulate $f$ exactly to a d.c. function. Because $f$ is separable and since the procedure is applied to each $f_i$ individually, we shall dispose of the index $i$ here and just consider $f(x)$ as a function of only one variable $x$. We suppose that $x \in [a, b]$ and aim to find functions $p, q: \mathbb{R} \to \mathbb{R}$ such that $p$ is concave, $q$ is convex and $f(x) = p(x) + q(x)$.

First of all find a partition $\mathcal{N} = \{[a_i, b_i] \mid i = 1, \dots, r\}$ of $[a, b]$ such that:
1. $a_1 = a, b_r = b$;

2. for all $i < r$, we have $b_i = a_{i+1}$;

3. for all $i < r$, either $f$ restricted to $[a_i, b_i]$ is convex and $f$ restricted to $[a_{i+1}, b_{i+1}]$ is concave, or the reverse is true (i.e. $f|_{[a_i, b_i]}$ is concave and $f|_{[a_{i+1}, b_{i+1}]}$ is convex).

Basically $a_2 = b_1, \dots, a_r = b_{r-1}$ are the points where $f$ changes between concavity and convexity<fnref n="2" />. Such a partition can be found by studying the behaviour of the derivatives $f'(x), f''(x)$. We assume without loss of generality that $f$ is concave on the first subinterval $[a_1, b_1]$ (if $f$ is convex on the first subinterval just set $a_1 = b_1$). This makes it possible to infer that, if $i$ is odd, then $f$ is concave on $[a_i, b_i]$; if $i$ is even, then $f$ is convex on $[a_i, b_i]$.

Now let:

$$ \Delta_i = \begin{cases} \min\{f'_{-}(b_i), f'_{+}(b_i)\} & \text{if } i \text{ is odd} \\ \max\{f'_{-}(b_i), f'_{+}(b_i)\} & \text{if } i \text{ is even} \end{cases} $$

where $f'_{-}(x)$ is the left derivative and $f'_{+}(x)$ is the right derivative of $f$ at $x$, and define two affine functions $s_i, t_i : [a, b] \to \mathbb{R}$ as follows:

$$ s_i(x) = f(b_i) + \Delta_i(x - b_i) $$

$$ t_i(x) = \sum_{k=1}^{i} (-1)^{k+i} s_k(x). $$

If we now define:

$$ p(x) = \begin{cases} f(x) - t_{i-1}(x) & \text{if } x \in D_i \text{ and } i \text{ is odd} \\ t_{i-1}(x) & \text{if } x \in D_i \text{ and } i \text{ is even} \end{cases} $$

$$ q(x) = \begin{cases} t_{i-1}(x) & \text{if } x \in D_i \text{ and } i \text{ is odd} \\ f(x) - t_{i-1}(x) & \text{if } x \in D_i \text{ and } i \text{ is even,} \end{cases} $$

then it can be shown [67] that $p$ is concave and $q$ is convex, and that $p(x) + q(x) = f(x)$, so that this reformulation is exact.

### 2.1.8 Factorable problems

Most common functions can be expressed in factorable form [80]; this form is desirable because it makes it relatively easy to construct convex relaxations in a recursive way [126, 141].

<footnote n="2">²If there are any cusp points, the function might fail to alternate in such a way; in this case we simply define an empty interval where $a_i = b_i$.</footnote>
A mathematical programming problem is in factorable form if it is written in the following way:

$$ \left. \begin{array}{l} \min_{x \in \mathbb{R}^n} \quad X^N(x) \\ -\infty < a_i \le X^i(x) \le b_i < \infty, \quad i \in [1, N-1] \end{array} \right\} \qquad (2.7) $$

where $x = (x_1, \dots, x_n)$, $X^i(x) = x_i$ for $1 \le i \le n$ and the other $X^i$ expressions are defined recursively as follows: given $X^p(x)$ for $p = 1, \dots, i-1$ then for $i = n+1, \dots, N$

$$ X^i(x) = \sum_{p=1}^{i-1} T_p^i(X^p(x)) + \sum_{p=1}^{i-1} \sum_{q=1}^{p} V_{q,p}^i(X^p(x)) U_{p,q}^i(X^q(x)) $$

where $T_p^i, U_p^i, V_p^i$ are continuous functions of one variable.

The algebraic properties of the factorable form make it possible to substitute each nonconvex term (univariate functions and bilinear products) with a convex relaxation of it [126]. The $\alpha$BB method [13, 2, 6, 5] is based on a variation of the factorable form problem, with extensions to include fractional terms and general twice-differentiable nonconvex terms (which are underestimated by means of a quadratic convex relaxation).

Sherali [113] applied the RLT (Reformulation-Linearization Technique) relaxation technique (see Section 2.4.1) to factorable problems.

### 2.1.8.1 Reformulation of factorable problems to separable form

In this reformulation [79], each nonseparable term in the objective function and the constraints is recursively replaced by an equivalent separable expression, until no nonseparable terms remain. The two steps to repeat are:

1. Replace any product term of the form $q_1(x)q_2(x)$ by $y_1^2 - y_2^2$ and add the constraints $q_1(x) = y_1 - y_2$ and $q_2(x) = y_1 + y_2$ to the definition of $\Omega$.

2. Replace any term of the form $T(t(x))$ by $T(y)$ and add the constraint $t(x) = y$ to the definition of $\Omega$.

It can be shown that the class of problems that can be reformulated to separability via the procedure above loosely corresponds to the class of factorable problems. In fact, it is possible
to expand this class to include terms like $q_1(x)^{q_2(x)}$: replace this term by a new variable $y$ and add the following constraints to the definition of $\Omega$:

$$ 
\begin{aligned}
y_1 &= q_1(x) \\
y_2 &= q_2(x) \\
y_3 &= \log(y_2) \\
y_4 &= y_1 + y_3 \\
y_5 &= \frac{1}{2}(y_4^2 - y_1^2 - y_3^2) \\
y &= e^{y_5}.
\end{aligned}
$$

### 2.1.9 Smith's standard form

This standard form, a symbolic exact reformulation for NLPs and MINLPs in general form, was first proposed in [114]. It will be discussed in detail in Section 5.2.2.1 as it forms a fundamental part of Smith's Branch-and-Bound algorithm, on which a considerable proportion of this thesis is based. This reformulation is an automatic process that isolates the nonlinear terms of the problem (1.1) in a list of simple constraints which are easy to deal with algorithmically. The following example will suffice, for now, to explain the essence of this reformulation.

#### 2.1.2 Example

In order to reduce the following problem:

$$ \min_{x \in \Omega} f(x_1, x_2) = \log(x_1 x_2) e^{(x_1+x_2)} + 2x_1 $$

to Smith's standard form, we define new problem variables $w_1, \dots, w_5$ via the following list of constraints:

$$ 
\begin{aligned}
w_1 &= x_1 x_2 \\
w_2 &= x_1 + x_2 \\
w_3 &= \log(w_1) \\
w_4 &= e^{w_2} \\
w_5 &= w_3 w_4.
\end{aligned}
$$

We can now express the objective function $f$ in terms of these new variables $f(w_5, x_1) = w_5 + 2x_1$, by adding the above list of constraints to the definition of the feasible region $\Omega$.
Smith's standard form is a lifting (it adds new variables to the problem), so it may not always be convenient; however, the advantage of Smith's standard form is that it is easier for an automatic algorithm to deal with a list of simple constraints like the one above, rather than with a nonlinear problem in the most general formulation (2.1). Such a reformulation makes it easy to implement symbolic algorithms: in short, it is a sort of “starting step” for more complex reformulations.

Smith's standard form can be carried out algorithmically in an extremely efficient way (see section 5.2.2.1). It can be shown that this reformulation is exact: solving a problem in Smith's standard form to global optimality will produce the same solutions as solving the original problem.

## 2.2 Exact reformulations

In this section, we shall present some exact reformulations which are not considered as standard forms.

### 2.2.1 Equality/inequality constrained problems

Sometimes it is convenient, for algorithmic reasons, to express an inequality constraint as an equality constraint, and vice versa. Given an inequality $g(x) \le 0$, it can be reformulated exactly as an equality via the introduction of a slack variable. Write $g(x) \le 0$ as $g(x) + s = 0$ where $s$ is a new problem variable (the slack variable) such that $s \ge 0$. This reformulation is a lifting.

Vice versa, any equation $g(x) = 0$ can be reformulated exactly to a pair of inequalities $g(x) \le 0$ and $g(x) \ge 0$, without adding any new variable to the problem.

### 2.2.2 Dimensionality reduction

Cantor proved that $\mathbb{R}^n$ and $\mathbb{R}$ have the same cardinality $c$ (the cardinality of the continuum) by showing that it was possible to construct a bijection between the two sets. Cantor's bijection, however, is not continuous in the usual topology. Peano and Hilbert suggested bijections (known as Peano-type space-filling curves) that were everywhere continuous and nowhere dif-
ferentiable. These bijections make it possible to devise an exact reformulation of a problem in $\mathbb{R}^n$ to a problem in $\mathbb{R}$, i.e. having a single problem variable.

Given a space-filling curve $g : \mathbb{R} \to [0, 1]^n$ (which exists by transfinite cardinality considerations [65, 24]), the dimensionality of the problem can be reduced from $n$ to 1 by solving the reduced problem $\min_{y \in \mathbb{R}} f(g(y))$. The function $f \circ g$ maps $\mathbb{R}$ into $\mathbb{R}$, and may thus be minimized via efficient 1-dimensional optimization methods. If $y$ is a minimizer of $f \circ g$, then $g(y)$ is a minimizer for $f$.

The difficult part is the choice of an appropriate space-filling function. These functions are usually a set of theoretical tools used to prove existence theorems, and as such are even difficult to describe explicitly. This approach may be worth investigating when the Euclidean space is approximated by a rational or integer lattice. There are some references in the literature which explain how to generate space-filling curves algorithmically [37]. One example is the Hilbert space-filling curve (see Fig. 2.1).

<figure><img src="image_4.png" /><figcaption>Figure 2.1: The Hilbert space-filling curve in 2-dimensional space.</figcaption></figure>

This technique is based on ideas proposed in the 1970s [122]. Interest in the use of space filling curves for global optimization waned during the 80s, but recently there have been renewed
efforts in dimensionality reduction [123, 44, 124].

## 2.3 Convex relaxations

A relaxation cannot be used to solve a difficult problem directly because the solution of the original problem cannot, in general, be directly inferred from the solution of the relaxation. Relaxations are, however, very important in the field of deterministic global optimization. One of the most important tools in this field is the Branch-and-Bound algorithm, which uses a convex (or linear) relaxation at each step to calculate the lower bound in a region.

Convex relaxations for nonconvex problems in form (2.1) are obtained by substituting the (nonconvex) objective function $f(x)$ with a convex relaxation $\underline{f}(x)$ and the (nonconvex) feasible region $\Omega$ with a convex set $\bar{\Omega}$ such that $\Omega \subseteq \bar{\Omega}$. Because every local minimizer of a convex problem is also a global minimizer, solving the convex relaxation:

$$ \min_{x \in \Omega} \underline{f}(x) \qquad (2.8) $$

with a local optimization algorithm will obtain the global solution of problem (2.8) which is guaranteed to be a valid lower bound for the global solution of the original problem (2.1). When the feasible region $\Omega$ is defined by equality and inequality constraints, as in equation (2.2), any convex relaxation $\bar{\Omega} \supseteq \Omega$ must be such that all the equality constraints are linear and all inequality constraints are convex<fnref n="3" />.

In general, there is more than one possible convex relaxation for any given problem; we therefore look for the best one, i.e. generally the one that gives rise to the greatest possible lower bound. Unfortunately, finding convex relaxations and convex envelopes of arbitrary subsets of $\mathbb{R}^n$ is not an easy task. Usually generic problems (1.1) need to be reformulated to some standard form before it is possible to construct a convex relaxation.

### 2.3.1 $\alpha$BB convex relaxation

Floudas and co-workers have proposed a Branch-and-Bound algorithm (called $\alpha$BB [2, 4, 5, 8]) for general nonconvex twice-differentiable problems. The algorithm aims to solve a problem in form (2.1) where the feasible region $\Omega$ is defined as in (2.2). Any nonlinear, twice-differentiable

<footnote n="3">³A more precise characterization is that inequality constraints should be quasi-convex.</footnote>
function $f(x)$ in the problem, be it the objective function or one of the constraints, can be reformulated exactly as:

$$f(x) = c^T x + f_C(x) + \sum_i b_i x_{B_1(i)} x_{B_2(i)} + \sum_i t_i x_{T_1(i)} x_{T_2(i)} x_{T_3(i)} + \\ \sum_i d_i \frac{x_{F_1(i)}}{x_{F_2(i)}} + \sum_i r_i \frac{x_{R_1(i)} x_{R_2(i)}}{x_{R_3(i)}} + \sum_i f_{U(i)}(x_i) + \sum_i f_{N(i)}(x),$$

where:

* $c, x \in \mathbb{R}^n$ and each $b_i, t_i, d_i, r_i$ is a real constant;

* $f_C(x)$ is a general convex function;

* $B_j, T_j, F_j, R_j, U_j, N_j$ are integer functions $\mathbb{N} \to \{1, \dots, n\}$.

* each $f_{U(i)}$ is a concave univariate function term;

* each $f_{N(i)}$ is a general nonconvex function term.

For a bilinear term $xy$, McCormick's underestimators [80] are used. A new variable $w_B$ is added to the problem (it replaces the bilinear term $xy$) and the following inequality constraints are inserted in the relaxed problem:

$$ \begin{align*} w_B &\ge x^L y + y^L x - x^L y^L \\ w_B &\ge x^U y + y^U x - x^U y^U \\ w_B &\le x^U y + y^L x - x^U y^L \\ w_B &\le x^L y + y^U x - x^L y^U. \end{align*} $$

The above linear inequalities have been shown to be the convex envelope of a bilinear term [10]. The maximum separation of the bilinear term $xy$ from its convex envelope $\max(x^L y + y^L x - x^L y^L, x^U y + y^U x - x^U y^U)$ inside the rectangle $[x^L, x^U] \times [y^L, y^U]$ occurs at the middle point $(\frac{x^L+x^U}{2}, \frac{y^L+y^U}{2})$ and is equal to $\frac{(x^U-x^L)(y^U-y^L)}{4}$ [13].
For a trilinear term $xyz$ a new variable $w_T$ is introduced, to replace the trilinear term $xyz$, together with the following constraints [75]:

$$w_T \geq xy^L z^L + x^L y^L z^L + x^L y^L z^L - 2x^L y^L z^L$$

$$w_T \geq xy^U z^U + x^U y^U z^U + x^U y^U z^U - x^U y^U z^U - x^U y^U z^U$$

$$w_T \geq xy^L z^L + x^L y^L z^L + x^L y^L z^L - x^L y^L z^L - x^L y^L z^L$$

$$w_T \geq xy^U z^U + x^U y^U z^U + x^U y^U z^U - x^U y^U z^U - x^U y^U z^U$$

$$w_T \geq xy^L z^L + x^L y^L z^L + x^L y^L z^L - x^L y^L z^L - x^L y^L z^L$$

$$w_T \geq xy^U z^U + x^U y^U z^U + x^U y^U z^U - x^U y^U z^U - x^U y^U z^U$$

$$w_T \geq xy^L z^L + x^L y^L z^L + x^L y^L z^L - x^L y^L z^L - x^L y^L z^L$$

$$w_T \geq xy^U z^U + x^U y^U z^U + x^U y^U z^U - 2x^U y^U z^U$$

Fractional terms $\frac{x}{y}$ are underestimated by replacing them with a new variable $w_F$ and adding two new constraints to the problem:

$$w_F \geq \begin{cases} \frac{x^L}{y} + \frac{x}{y^U} - \frac{x^L}{y^U} & \text{if } x^L \geq 0 \\ \frac{x}{y^U} - \frac{x^L y}{y^L y^U} + \frac{x^L}{y^L} & \text{if } x^L < 0 \end{cases} \quad w_F \geq \begin{cases} \frac{x^U}{y} + \frac{x}{y^L} - \frac{x^U}{y^L} & \text{if } x^U \geq 0 \\ \frac{x}{y^L} - \frac{x^U y}{y^L y^U} + \frac{x^U}{y^U} & \text{if } x^U < 0 \end{cases}$$

Fractional trilinear terms $\frac{xy}{z}$ can be underestimated by replacing them by a new variable $w_{FT}$ and adding the new constraints (for $x^L, y^L, z^L \geq 0$):

$$w_{FT} \geq \frac{xy^L}{z^U} + \frac{x^L y}{z^U} + \frac{x^L y^L}{z} - 2 \frac{x^L y^L}{z^U}$$

$$w_{FT} \geq \frac{xy^L}{z^U} + \frac{x^L y}{z^L} + \frac{x^L y^U}{z} - \frac{x^L y^U}{z^L} - \frac{x^L y^L}{z^U}$$

$$w_{FT} \geq \frac{xy^U}{z^L} + \frac{x^U y}{z^U} + \frac{x^U y^L}{z} - \frac{x^U y^L}{z^U} - \frac{x^U y^L}{z^L}$$

$$w_{FT} \geq \frac{xy^U}{z^U} + \frac{x^U y}{z^L} + \frac{x^L y^U}{z} - \frac{x^L y^U}{z^U} - \frac{x^U y^L}{z^L}$$

$$w_{FT} \geq \frac{xy^L}{z^U} + \frac{x^L y}{z^L} + \frac{x^U y^L}{z} - \frac{x^U y^L}{z^L} - \frac{x^L y^L}{z^U}$$

$$w_{FT} \geq \frac{xy^U}{z^U} + \frac{x^U y}{z^L} + \frac{x^L y^U}{z} - \frac{x^L y^U}{z^U} - \frac{x^U y^L}{z^L}$$

$$w_{FT} \geq \frac{xy^L}{z^U} + \frac{x^L y}{z^L} + \frac{x^U y^L}{z} - \frac{x^U y^L}{z^L} - \frac{x^L y^L}{z^U}$$

$$w_{FT} \geq \frac{xy^U}{z^L} + \frac{x^U y}{z^L} + \frac{x^U y^U}{z} - 2 \frac{x^U y^U}{z^L}$$
To relax a concave univariate function $f_{U(i)}(x)$ over $[x^L, x^U]$, the $\alpha$BB algorithm uses a chord underestimator:

$$f(x^L) + \frac{f(x^U) - f(x^L)}{x^U - x^L}(x - x^L).$$

The main innovation of the $\alpha$BB algorithm is in the underestimation of a general nonconvex function term $f_{N(i)}(x)$. This is underestimated over the entire domain $[x^L, x^U] \subseteq \mathbb{R}^n$ by the function $L(x)$ defined as follows:

$$L(x) = f(x) + \sum_{i=1}^{n} \alpha_i (x_i^L - x_i)(x_i^U - x_i)$$

where the $\alpha_i$ are positive scalars that are sufficiently large to render the underestimating function convex. A good feature of this kind of underestimator is that, unlike other underestimators, it does not introduce any new variable or constraint, so that the size of the relaxed problem is the same as the size of the original problem regardless of how many nonconvex terms it involves.

Since the sum $\sum_{i=1}^n \alpha_i (x_i^L - x_i)(x_i^U - x_i)$ is always negative, $L(x)$ is an underestimator for $f(x)$. Furthermore, since the quadratic term is convex, all nonconvexities in $f(x)$ can be overpowered by using sufficiently large values of the $\alpha_i$ parameters. From basic convexity analysis, it follows that $L(x)$ is convex if and only if its Hessian matrix $H_L(x)$ is positive semi-definite. Notice that:

$$H_L(x) = H_f(x) + 2\Delta$$

where $\Delta \equiv \text{Diag}_{i=1}^n(\alpha_i)$ is the matrix with $\alpha_i$ as diagonal entries and all zeroes elsewhere (diagonal shift matrix). Thus the main focus of the theoretical studies concerning all $\alpha$BB variants is on the determination of the $\alpha_i$ parameters. Some methods are based on the simplifying requirement that the $\alpha_i$ are chosen to be all equal (uniform diagonal shift matrix), others reject this simplification (non-uniform diagonal shift matrix). Under the first condition, the problem is reduced to finding the parameter $\alpha$ that makes $H_L(x)$ positive semi-definite. It has been shown that $H_L(x)$ is positive semi-definite if and only if:

$$\alpha \geq \max\{0, -\frac{1}{2} \min_{i, x^L \leq x \leq x^U} \lambda_i(x)\}$$

where $\lambda_i(x)$ are the eigenvalues of $H_f(x)$. Thus the problem is now of finding a lower bound on the minimum eigenvalue of $H_f(x)$. The most promising method to this end seems to be Interval Matrix Analysis. Various $o(n^2)$ and $o(n^3)$ methods have been proposed to solve both the uniform and the non-uniform diagonal shift matrix problem [39].

Thus, having constructed a convex underestimating function for the reformulated function $f(x)$, the relaxation of the problem is carried out accordingly, bearing in mind that:
* the objective function is replaced by its convex underestimator;

* any nonlinear equality constraint $h_i(x) = 0$ is replaced by the two inequality constraints $h_i(x) \le 0$ and $-h_i(x) \le 0$;

* any nonconvex inequality constraint $g_i(x) \le 0$ is reformulated to $g_i(x) \le 0$ where $g_i$ is the convex underestimator for $g_i(x)$.

### 2.3.2 Smith's convex relaxation

In the work of Smith [114, 116], Smith's standard form (see sections 2.1.9, 5.2.2.1) is used as a starting point for forming the convex relaxation. All the nonconvex terms are isolated in constraints of the forms $z = f(x)$ or $z = g(x, y)$, where $f$ is a univariate nonlinear function, $g$ is a bivariate nonlinear function and $x, y, z \in \mathbb{R}$ are single problem variables. The smallest convex set containing $F = \{(x, z) | z = f(x)\} \subseteq \mathbb{R}^2$ is given by $\bar{F} = \{(x, z) | z \ge \underline{f}(x) \wedge z \le \bar{f}(x)\}$ where $\underline{f}, \bar{f}$ are respectively the convex and the concave envelopes of $f$. Similarly for $g$, the smallest convex set containing $G = \{(x, y, z) | z = g(x, y)\} \subseteq \mathbb{R}^3$ is given by $\bar{G} = \{(x, y, z) | z \ge \underline{g}(x, y) \wedge z \le \bar{g}(x, y)\}$ where $\underline{g}, \bar{g}$ are respectively the convex and the concave envelopes of $g$. When the functions $f$ are convex, concave univariate terms and the functions $g$ are bilinear, trilinear or fractional terms, the convex envelopes are the same as those listed in Section 2.3.1. However, Smith's reformulation makes it possible to extend convex relaxations to nonlinear terms which are not twice-differentiable (as opposed to the $\alpha$BB relaxation).

Smith, in his work, did not construct the convex/concave envelopes of a piecewise convex and concave term (like, for example, a monomial of odd power when the range of the variable includes zero); Chapter 4 addresses this problem.

The disadvantage of Smith's relaxation is that it is based on Smith's standard form, which is a lifting, and therefore adds new variables to the problem (one for each nonlinear term in the problem). This may result in an excessively high number of problem variables.

### 2.3.3 BARON's convex relaxation

BARON stands for “Branch And Reduce Optimization Navigator”. It is a global optimization software written by Sahinidis and co-workers that relies on a Branch-and-Bound algorithm to
solve factorable problems to global optimality [102]. The lower bound to the objective function in each region of the Branch-and-Bound tree is calculated by means of a convex relaxation [126]. The techniques used to form the nonlinear convex relaxation of factorable problems include all the standard convex envelopes for nonconvex factorable terms found in Section 2.3.1 (apart from the α-parameter underestimation for general twice-differentiable nonconvex terms, which is not part of BARON). The most innovative features of BARON in terms of convex relaxations are the following.

* Specific mention of piecewise convex and piecewise concave univariate terms (called concavoconvex by the authors) and the respective convex and concave envelopes (see also the detailed study in Chapter 4) [126]. An alternative to this envelope is suggested which circumvents the issue: by branching on the concavoconvex variable at the point where the curvature changes (i.e. the point where the concavoconvex term changes from concave to convex or vice versa) at a successive Branch-and-Bound iteration, the term becomes completely concave and completely convex in each region [102].

* Convex and concave envelopes are suggested for various types of fractional terms, based on the theory of convex extensions [127, 128]. The proposed convex underestimator for the term $\frac{x}{y}$, where $x \in [x^L, x^U]$ and $y \in [y^L, y^U]$ are strictly positive, is as follows:

$$ 
\left.
\begin{array}{l}
z \ge \frac{x^L}{y_a}(1-\lambda) + \frac{x^U}{y_b}\lambda \\
y^L \le y_a \le y^U \\
y^L \le y_b \le y^U \\
y = (1-\lambda)y_a + \lambda y_b \\
x = x^L + (x^U - x^L)\lambda \\
0 \le \lambda \le 1
\end{array}
\right\} (2.9)
$$

The underestimator is modified slightly when $0 \in [x^L, x^U]$:

$$ 
\left.
\begin{array}{l}
z \ge \frac{x^L(y^L+y^U-y_a)}{y^L y^U}(1-\lambda) + \frac{x^U}{y_b}\lambda \\
y^L \le y_a \le y^U \\
y^L \le y_b \le y^U \\
y = (1-\lambda)y_a + \lambda y_b \\
x = x^L + (x^U - x^L)\lambda \\
0 \le \lambda \le 1
\end{array}
\right\} (2.10)
$$

It is shown that these underestimators are tighter than all previously proposed convex underestimators for fractional terms, in particular:
- the bilinear envelope:

$$ \max \left\{ \frac{xy^U - yx^L + x^L y^U}{(y^U)^2}, \frac{xy^L - yx^U + x^U y^L}{(y^L)^2} \right\} $$

- the nonlinear envelope:

$$ \frac{1}{y} \left( \frac{x + \sqrt{x^L x^U}}{\sqrt{x^L} + \sqrt{x^U}} \right)^2 $$

* Having constructed a reasonably tight nonlinear convex relaxation, the authors discuss various outer approximation techniques to reformulate this to a linear relaxation which can be solved by using a very fast and efficient LP software [126].

BARON is a very efficient global optimization software, paying particular attention to several important implementational aspects, such as: generating valid cuts during pre-processing (optimality bounds tightening [102, 126]) and during execution (feasibility bounds tightening and range reduction techniques [99, 100]); improving the branching scheme [102, 126, 130, 101]; and most importantly of all, targeting particular problem formulations with specialized solvers [126]. These are available for:

* mixed-integer linear programming;

* separable concave quadratic programming;

* indefinite quadratic programming;

* separable concave programming;

* linear multiplicative programming;

* general linear multiplicative programming;

* univariate polynomial programming;

* 0-1 hyperbolic programming;

* integer fractional programming;

* fixed charge programming;

* problems with power economies of scale;

besides the “default” solver for general nonconvex factorable problems.
## 2.4 Linear relaxations

It is possible to use linear over- and underestimators for each nonlinear term in a nonconvex NLP in order to obtain a linear relaxation of the problem. Because a linear problem is always convex, the convexity properties that guarantee the validity of a lower bound remain true. The advantage of a linear relaxation with respect to a convex (possibly nonlinear) relaxation is that linear optimization software can be employed to solve the relaxed underestimating problem. Linear optimization codes are much more efficient than nonlinear optimization software; hence the overall run of a Branch-and-Bound algorithm might be faster. The disadvantage, however, is that a linear relaxation might not be a convex envelope: hence the convexity gap might be increased and the lower bound to the original problem might not be the best possible.

### 2.4.1 Reformulation-linearization technique

The basic idea of the Reformulation-Linearization Technique (RLT), proposed by Sherali and co-workers in a number of papers [109, 111, 107, 113, 112, 110, 108], consists in deriving valid cuts to the problem by multiplying together various factors involving variables and constraints. This technique was initially proposed in conjunction with combinatorial problems with binary variables, and then extended to continuous bilinear problems, polynomial problems and factorable problems.

In this section we shall describe the RLT applied to bilinear problems [111] of the form:

$$ \left. \begin{array}{l} \min_x \quad x^T Q x + c^T x \\ Ax = b \\ x^L \le x \le x^U, \end{array} \right\} \qquad (2.11) $$

where $c, x \in \mathbb{R}^n$, $Q$ is an $n \times n$ matrix, $A$ is an $m \times n$ matrix and $b \in \mathbb{R}^m$.

In order to form a linear (convex) relaxation of problem (2.11), the RLT applied to bilinear problems considers the following sets of algebraic expressions:

* the bound factor set $B_F = \{x_i - x_i^L \mid i \le n\} \cup \{x_i^U - x_i \mid i \le n\}$;

* the constraint factor set $C_F = \{\sum_{j=1}^n a_{ij}x_j - b_i \mid i \le m\}$.

Note that for each $\beta \in B_F$ the constraint $\beta \ge 0$ is a valid problem constraint, and so is $\gamma = 0$ for all $\gamma \in C_F$.
The RLT procedure for forming the convex relaxation consists in creating new linear valid constraints (reformulation step) by multiplying together bound factors and constraint factors as follows:

1. for all $\beta_1, \beta_2 \in B_F$, $\beta_1\beta_2 \ge 0$ is a valid constraint (generation via bound factors);

2. for all $\beta \in B_F$ and for all $\gamma \in C_F$, $\beta\gamma = 0$ is a valid constraint (mixed generation);

3. for all $\gamma_1, \gamma_2 \in C_F$, $\gamma_1\gamma_2 = 0$ is a valid constraint (generation via constraint factors).

Having created all these new constraints, we define new variables $w_j^i = x_i x_j$ for all $i, j$ between 1 and $n$, and use them to replace the corresponding bilinear products appearing in problem (2.11) or in the newly generated constraints (linearization step). Assuming there are $t$ distinct bilinear terms$^4$, we end up with a linear relaxation whose variable vector $(x, w)$ is in $\mathbb{R}^{n+t}$.

Let $S_F$ be the region defined by the newly generated constraints. The linear relaxation of problem (2.11) is as follows:

$$ \left. \begin{array}{llll} \min_{x} & c^T x & + & p^T w \\ & Ax & = & b \\ & (x, w) & \in & S_F \\ & x^L \le & x & \le x^U \\ & w^L \le & w & \le w^U \end{array} \right\} \qquad (2.12) $$

where $w^L, w^U \in \mathbb{R}^t$ are the variable bounds on the $w$ variables (obtained by simple interval arithmetic on the bounds of the $x$ variables via the defining relations $w_j^i = x_i x_j$).

It is also possible to derive RLT constraints from a system of inequalities $Hx \le d$, where $H$ is an $m' \times n$ matrix and $d \in \mathbb{R}^{m'}$. Define the inequality constraint factor set $C'_F = \{\sum_{j=1}^n h_{ij}x_j - d_i \mid i \le m'\}$ and notice that for all $\gamma' \in C'_F$, we have $\gamma' \le 0$. Since all the bound factors $\beta \in B_F$ are nonnegative, any constraint of the form $\beta\gamma' \le 0$ is a valid problem constraint.

When the RLT is applied to general polynomial programming problems, arbitrary multilinear products:

$$ \prod_j \beta_{ij} \prod_j \gamma_{ij} $$

may arise between the bound factors $\beta_{ij}$ and constraint factors $\gamma_{ij}$.

<footnote n="4">⁴Note that $t \le \frac{n(n+1)}{2}$.</footnote>
The RLT gives rise to very tight convex relaxations but it has a disadvantage: it is a heuristic reformulation. If implemented to its full extent (i.e. adding all possible constraints deriving from all possible bound and constraint factor products), it gives rise to relaxations which have too many constraints and therefore take inordinate amounts of time to solve. More precisely, with 2*n* bound factors and *m* constraint factors, we obtain:

* $4n^2$ possible bound factor products;

* $m^2$ possible constraint factor products;

* $2mn$ possible mixed factor products,

for a total of $4n^2 + 2mn + m^2$ new constraints. When the RLT is applied to polynomial problems, more added constraints can be obtained recursively from the newly formed factor product constraints and the original factors, so this number increases even more.

## 2.5 Other reformulations

There are several other reformulation studies in the literature. These are not analysed in detail here either because they are specific to a certain type of optimization problem (and, hence, difficult to generalize) or because they address software implementation issues. Here follows a short list of some of the most notable such works.

A qualitative knowledge-based software implementation, REFORM, for the generic reformulation of MINLPs, is proposed in [12]; it is aimed at increasing the general robustness of an optimization problem before submitting it to an optimization software. Frangioni [41] showed how to reformulate a bilevel programming problem (a hierarchical two-stage optimization problem) as a MILP. Fischer [36] reformulated the complementarity problem as a minimization problem with nonnegativity constraints. Some regularity conditions ensure that a stationary point of the reformulation is a solution of the original problem. Bomze and co-workers [20, 21] reformulated the maximum clique problem as a quadratic problem over a standard simplex, and used copositivity-based procedures to solve the standardized quadratic problem. In [19], a general method for handling disjunctive constraints in MINLPs is proposed, which relies on reformulating disjunctive constraints to a mathematical form that can be solved using standard MILP software. Sahinidis and co-workers [1, 129] have studied various reformulations of the pooling problem extensively.
## 2.6 Conclusion

This chapter has provided a literature review of various techniques for the reformulation of optimization problems, with emphasis on the most useful and generally applicable reformulation techniques for deterministic global optimization. This includes exact reformulations leading to standard forms (box-constrained problems, separable problems, binary problems, concave problems, d.c. problems, factorable problems, generic nonlinear problems), exact reformulations which address a specific feature of a problem, and relaxations (convex and linear).

It is clear from the review presented in this chapter that much progress has been achieved in both exact reformulation and convex relaxation of nonconvex NLPs, especially over the past decade. Equally clearly, there are several areas for which no satisfactory solutions exist to date. Two of these will receive special attention in this thesis; these are respectively:

1. *The tight linear relaxation of bilinear terms.* Bilinear terms are almost ubiquitous in practical applications. The RLT technique described in Section 2.4.1 represents a good attempt to tighten the convex relaxation of NLPs involving such terms. However, the large numbers of constraints and variables that it introduces represents a significant obstacle to its use in large problems. An algorithm which can be applied to large sparse NLPs and which leads to convex relaxations that are both tighter and smaller is presented in Chapter 3.

2. *The tight convex relaxation of monomials of odd power.* Although such terms appear frequently in applications (e.g. as cubics or quintics), no convex envelope has been proposed for them to date. As a result, they are currently treated either as general non-convexities (cf. Section 2.3.1), or via specialized branching schemes (cf. Section 2.3.3). Neither of these two approaches is satisfactory. A convex envelope that addresses this issue in a more efficient manner is presented in Chapter 4 of this thesis.
# Chapter 3

## Reduction constraints for sparse bilinear programs

We consider the solution of Nonlinear Programs (NLPs) of the following standard form:

$$ \min_{z} z_{l} \qquad (3.1) $$

$$ Az = b \qquad (3.2) $$

$$ z_i = z_j z_k \quad \forall i, j, k \in B \qquad (3.3) $$

$$ z_i = \frac{z_j}{z_k} \quad \forall i, j, k \in F \qquad (3.4) $$

$$ z_i = f_i(z_j) \quad \forall i, j \in N \qquad (3.5) $$

$$ z^L \leq z \leq z^U \qquad (3.6) $$

where $z = (z_1, \dots, z_p) \in \mathbb{R}^p$ are the problem variables, $l$ is an index in the set $\{1, \dots, p\}$, $A = (a_{ij})$ is an $m \times p$ matrix of rank $m$, $b \in \mathbb{R}^m$, $B, F$ are sets of index triplets $\{(i, j, k) \mid 1 \le i, j, k \le p\}$, $N$ is a set of index pairs $\{(i, j) \mid 1 \le i, j \le p\}$, $f_i : \mathbb{R} \to \mathbb{R}$ are nonlinear univariate functions and $z^L, z^U \in \mathbb{R}^p$ are variable bounds. The above standard form is practically important as it can be shown that all NLPs can automatically be reformulated to it using symbolic manipulations [116, 114] (see Chapter 5). Therefore, any theoretical results, manipulations and solution algorithms derived on the basis of this standard form are generally applicable.

Spatial Branch-and-Bound (sBB) algorithms [133] are among the most effective methods
currently available for the global solution of nonconvex NLPs (see Chapter 1). An important requirement for any sBB algorithm is to be able to construct a tight convex underestimator of the NLP within any given region of the space of the variables. For the standard form [P], the lower bound to the objective function can be generated by replacing the nonconvex constraints (3.3), (3.4), (3.5) with their convex relaxations; in this manner, a convex relaxation of the whole problem can be obtained in a simple and completely automatic fashion.

Tight convex relaxations for the most common nonconvex terms are available in the literature. One of the best known is that proposed by McCormick (1976) for the relaxation of bilinear terms. Linear fractional terms like those appearing in constraint (3.4) are reformulated to $z_i z_k = z_j$ and replaced by the McCormick convex relaxation. For the nonlinear terms in constraint (3.5), the convex relaxation depends on the function $f_i$. When $f_i$ is wholly concave or wholly convex, the function itself and its secant provide the tightest convex relaxation. For functions $f_i$ which are partially concave and partially convex, like e.g. $z_i = z_j^{2k+1}$, where $k \in \mathbb{N}$, building the convex relaxation may not be so straightforward (see, for example, [72]).

### 3.0.1 Example

*Consider the problem*

$$ 
\begin{array}{ll}
\min_{x,y,z,w} & z \\
z - x + w & = \ 1 \\
w & = \ xy \\
(x^L, y^L, z^L, w^L) & \le (x, y, z, w) \le (x^U, y^U, z^U, w^U)
\end{array}
$$

in standard form. To obtain the convex relaxation, we replace the nonconvex constraint $w = xy$ with its McCormick convex relaxation:

$$ 
\left. 
\begin{array}{l}
w \ge x^L y + y^L x - x^U y^U \\
w \ge x^U y + y^U x - x^U y^L \\
w \le x^U y + y^L x - x^U y^L \\
w \le x^L y + y^U x - x^L y^U
\end{array}
\right\} . \qquad (3.7)
$$

This chapter presents a technique to automatically reformulate optimization problems in standard form [P] in such a way that some of the nonlinear constraints (3.3) are replaced by linear constraints. This is possible because, in certain instances, feasible regions described by
nonlinear constraints are, in fact, (linear) hyperplanes. We propose an automatic algorithm to identify such instances in large scale systems. The solution of the reformulated problems by sBB algorithms often requires reduced computational effort, sometimes even by several orders of magnitude. This occurs because replacing nonlinear constraints with linear ones makes the feasible region of the convex relaxation is much tighter.

The creation of new linear constraints via multiplication of existing linear constraints by problem variables was proposed in [111, 108] under the name “reformulation-linearization technique” (RLT). The RLT uses linear constraints built in this way to provide a lower bound to bilinear programming problems. The maximum possible number of new linear constraints is created by multiplying all linear constraints by all problem variables; thus, the method may lead to excessive computational complexity in large problems. In contrast, the algorithm presented in this chapter identifies precisely the set of multiplications of linear constraints by variables that are beneficial.

The rest of this chapter is organized as follows. In Section 3.1, we introduce the basic concepts and ideas behind reduction constraints. Section 3.2 considers the existence of reduction constraints in a more formal manner. This provides the basis of the fast algorithm for the identification of reduction constraints presented in Section 3.3. An example of the application of the algorithm is presented in detail in Section 3.4. The effects of reduction constraints on the global solution of an important class of problems with bilinear constraints, namely pooling and blending problems, are considered in Section 3.5. Finally, Section 3.6 proposes an extension of the algorithm of Section 3.3 that may result in better reformulations at the expense of higher computational time and memory.

## 3.1 Basic concepts

This section introduces in an informal manner the motivation and general ideas behind the work presented in this chapter. Let $y$ be a single problem variable, and let $w = (w_1, \dots, w_n)$ and $x = (x_1, \dots, x_n)$ be problem variables (with $n < p$) such that the following constraints exist in the problem:

$$ \forall j \le n (w_j = x_j y). \qquad (3.8) $$
Now suppose that the problem also involves the linear constraint (cf. equation (3.2)):

$$ \sum_{j=1}^{n} a_{ij}x_j = b_i \qquad (3.9) $$

Multiplying this constraint by $y$ and making use of (3.8) leads to a new linear constraint:

$$ \sum_{j=1}^{n} a_{ij}w_j - b_iy = 0. \qquad (3.10) $$

The linear constraint (3.10) is redundant with respect to the original constraints. Indeed, it can be used to replace one of the bilinear constraints $w_j = x_jy$ in problem [P] without affecting the feasible region of the problem. To see this, assume $a_{ik} \neq 0$ for some $k \in \{1, \dots, n\}$ and discard the nonlinear constraint $w_k = x_ky$ from the set (3.8) above. We now replace $b_i$ in (3.10) with the left hand side of (3.9) to obtain:

$$ \sum_{j=1}^{n} a_{ij}w_j - \sum_{j=1}^{n} a_{ij}x_jy = \sum_{j=1}^{n} a_{ij}(w_j - x_jy) = 0. $$

Since $w_j = x_jy$ for all $j \neq k$, the above reduces to $a_{ik}(w_k - x_ky) = 0$, which implies $w_k = x_ky$. We have thus recovered the discarded bilinear constraint from the other $n-1$ bilinear constraints and the new linear constraint (3.10). We call the latter a *reduction constraint* as it reduces the number of bilinear terms in the original problem.

The geometric significance of reduction constraints is that, given any $i \leq m$, the set:

$$ \{(w, x, y) \mid \forall j \leq n (w_j = x_jy) \land \sum_{j=1}^{n} a_{ij}x_j = b_i\} \qquad (3.11) $$

is less nonlinear than might appear to be the case: in fact, provided that $a_{ik} \neq 0$ for some $k \leq n$, it is equal to the set:

$$ \{(w, x, y) \mid \forall j \neq k (w_j = x_jy) \land \sum_{j=1}^{n} a_{ij}x_j = b_i \land \sum_{j=1}^{n} a_{ij}w_j - b_iy = 0\}, \qquad (3.12) $$

where the $k$th bilinear constraint has been replaced by a linear constraint. Consequently, the convex relaxation of the set (3.12) is tighter than that of set (3.11), which may be important in the context of sBB algorithms.

Below, we consider an extreme example of such a reformulation.
<figure><img src="image_5.png" /><figcaption>Figure 3.1: Linearity embedded in a bilinear constraint.</figcaption></figure>

### 3.1.1 Example

Let $(w, x, y) \in \mathbb{R}^3$ and consider the set:

$$A_1 = \{(w, x, y) | w = xy \land x = 1\}.$$

This set is defined as the intersection of the bilinear surface $xy$ and the vertical plane $x = 1$ (see Fig. 3.1). By multiplying $x = 1$ by $y$, we get the linear constraint $w - y = 0$ (the skew plane in Fig. 3.1). It is evident that the set

$$A_2 = \{(w, x, y) | x = 1 \land w - y = 0\}$$

is the same as the set $A_1$. However, $A_2$ is defined only by linear relationships whereas $A_1$ is not. The convexification of $A_1$ in the context of an sBB algorithm (e.g. by introducing the McCormick relaxation of $xy$) would be both unnecessary and counter-productive.

A limitation of the reformulation technique presented above is that all constraints of the type $w_j = x_j y$ have to exist in the problem (3.1) before the linear constraint (3.10) can be created. Consider, however, a problem that already includes bilinear constraints $w_j = x_j y$ for
all $j \le n - 1$, and two linear constraints:

$$ \sum_{j=1}^{n-1} a_{ij}x_j + a_{in}x_n = b_i \qquad (3.13) $$

$$ \sum_{j=1}^{n-1} a_{hj}x_j + a_{hn}x_n = b_h \qquad (3.14) $$

where $i, h \le m$ and $a_{in}, a_{hn}$ are nonzero constants. On multiplying (3.13) and (3.14) by $y$, we obtain the two linear constraints:

$$ \sum_{j=1}^{n-1} a_{ij}w_j + a_{in}w_n - b_iy = 0 \qquad (3.15) $$

$$ \sum_{j=1}^{n-1} a_{hj}w_j + a_{hn}w_n - b_hy = 0. \qquad (3.16) $$

where $w_n$ is a new variable defined via a new bilinear constraint $w_n = x_n y$. By forming a suitable linear combination of the two constraints, the new variable $w_n$ can be eliminated, and we now have a linear constraint similar to (3.10) which can be used to eliminate one of the original bilinear constraints in the manner indicated earlier.

## 3.2 Fundamental properties

Section 3.1 has shown that it is possible, in some cases, to replace nonconvex nonlinear constraints by linear ones. We now proceed to provide a more rigorous foundation for these ideas.

Consider an optimization problem in standard form [P] and subsets of problem variables $w, x \in \mathbb{R}^n$ and $y \in \mathbb{R}$. Suppose the problem includes the set of $m$ linear equality constraints $Ax = b$, where the matrix $A$ is of full row rank $m$.

Multiplying $Ax = b$ by $y$, we create $m$ reduction constraints of the form $Aw - yb = 0$, where the variable vector $w \in \mathbb{R}^n$ is defined as in equation (3.8). The following theorem shows that the reduction constraints can be used to replace $m$ of the $n$ bilinear constraints in (3.8) without changing the feasible region of the problem.
### 3.2.1 Theorem

Let $J \subseteq \{1, \dots, n\}$ be an index set of cardinality $|J| = n - m$ and consider the sets:

$$C = \{(x, w, y) \mid Ax = b \wedge \forall j \le n (w_j = x_j y)\} \qquad (3.17)$$

$$R_J = \{(x, w, y) \mid Ax = b \wedge Aw - yb = 0 \wedge \forall j \in J (w_j = x_j y)\} \qquad (3.18)$$

Then, there is at least one set $J$ such that $C = R_J$.

*Proof.* The fact that $C \subseteq R_J$ for any set $J$ is straightforward: if $(x, w, y) \in C$, then it also satisfies the constraints $Aw - yb = 0$, hence $(x, w, y) \in R_J$.

We shall now prove the converse inclusion. Since $Ax = b$, the reduction constraint system $Aw - yb = 0$ implies $Aw - yAx = 0$. If we now define $u_j = w_j - x_j y, \forall j \le n$, this can be written as the linear homogeneous system $Au = 0$. Since rank($A$) = $m$, there exists a permutation $\pi$ of the columns of $A$ such that:

$$(A_0 \ A_1) \begin{pmatrix} u^{(0)} \\ u^{(1)} \end{pmatrix} = 0,$$

where $A_0$ is a non-singular $m \times m$ matrix, and $(u^{(0)})^T, (u^{(1)})^T$ is the corresponding permutation of $u$ with $u^{(0)} \in \mathbb{R}^m$ and $u^{(1)} \in \mathbb{R}^{n-m}$. Let $J$ be the image of $\{n - m + 1, \dots, n\}$ under the permutation $\pi$, and let $(x, w, y) \in R_J$. Since $(x, w, y)$ satisfies $w_j = x_j y$ for all $j \in J$, then $u^{(1)} = 0$; since $(x, w, y)$ also satisfies the reduction constraint system, we have $Au = 0$, which implies $u^{(0)} = 0$ as well. Hence $w_j = x_j y$ for all $j \le n$. $\square$

The above theorem concerns the multiplication of a set of linear constraints by a single variable $y$. The latter could, in fact, be any one of the system variables $x$, so considering multiplications by variable $x_k$, we can modify the definitions (3.17) and (3.18) to:

$$C^k = \{(x, w^k) \mid Ax = b \wedge \forall j \le n (w_j^k = x_j x_k)\}, \quad k \le n \qquad (3.19)$$

$$R_J^k = \{(x, w^k) \mid Ax = b \wedge Aw^k - x_k b = 0 \wedge \forall j \in J (w_j^k = x_j x_k)\}, \quad k \le n \quad (3.20)$$

where $w^k = (w_1^k, \dots, w_n^k)$. Now, by virtue of the theorem, we will have $C^k = R_J^k, \forall k \le n$. We note that the index set $J$ appearing in (3.20) is the same for all $k$: the proof of the theorem indicates that $J$ is fully determined by the nature of the matrix $A$ and does not depend on the specific variable $y$ being used to multiply the constraints. Consequently, the union of the sets $C^k$ over all $k \le n$ must be equal to the corresponding union of sets $R_J^k$. This proves that the
two sets:

$$ \bar{C} = \{(x, w) \mid Ax = b \wedge \forall k \le n \ j \le n \ (w_j^k = x_j x_k)\} \quad (3.21) $$

$$ \bar{R}_J = \{(x, w) \mid Ax = b \wedge Aw - x_k b = 0 \wedge \forall k \le n, \ \forall j \in J (w_j^k = x_j x_k)\} \quad (3.22) $$

are equal to each other.

The geometrical implications of the above results are that the intersection in $\mathbb{R}^n$ of a set of bilinear terms, and the linear form $Ax = b$ is a hypersurface containing a degree of linearity. By exploiting this linearity, we are able to replace some of the bilinear terms with linear constraints.

## 3.3 An algorithm for the identification of valid reduction constraints

Section 3.2 has established that, in principle, reduction constraints exist for any NLP involving bilinear terms and linear constraints. Here we are concerned with two more practical questions. Section 3.3.1 introduces a criterion for determining whether the introduction of reduction constraints is likely to be beneficial. Then sections 3.3.2 and 3.3.3 present an efficient algorithm for the identification of such reduction constraints based on the above criterion; this is particularly important in the case of large NLPs for which neither identification by inspection nor indiscriminate multiplication of linear constraints by problem variables are practical propositions.

### 3.3.1 Valid reduction constraint sets

Given an NLP in the standard form [P], consider a subset $\mathcal{L}$ of its linear constraints (3.2). Let the set of variables that occur in these constraints be denoted by $\mathcal{V}(\mathcal{L})$, i.e. $\mathcal{V}(\mathcal{L}) = \{k \mid \exists l \in \mathcal{L} (a_{lk} \neq 0)\}$.

Now consider multiplying the linear constraint set $\mathcal{L}$ by a problem variable $z_l$. This will create bilinear terms of the form $z_l z_k, \forall k \in \mathcal{V}(\mathcal{L})$. Some of these terms will already occur in the NLP, i.e. $\exists i : (i, l, k) \in B$. Let $K(l, \mathcal{L})$ be the subset of $\mathcal{V}(\mathcal{L})$ that leads to new bilinear terms, i.e.:

$$ K(l, \mathcal{L}) = \{k \mid k \in \mathcal{V}(\mathcal{L}) \wedge \exists i : (i, l, k) \in B\} $$
The theorem of Section 3.2 indicates that we can now construct a problem that has the same feasible region as [P] by replacing $|\mathcal{L}|$ bilinear constraints by linear constraints. The latter are said to form a valid reduction constraint system if the substitution leads to a reduction in the overall number of bilinear constraints in comparison with the original problem [P], i.e. if the number of new bilinear constraints introduced by the multiplication by $z_l$ is smaller than the number of bilinear constraints eliminated by the reduction constraints:

$$|K(l, \mathcal{L})| < |\mathcal{L}| \qquad (3.23)$$

Despite the apparent simplicity of this criterion, applying it in a brute force manner as a means of identifying reduction constraint systems is impractical: for an NLP containing $M$ linear constraints, one would have to examine $2^M - 1$ subsets $\mathcal{L}$ for each candidate multiplier variable $z_l$.

### 3.3.2 A graph-theoretical representation of linear constraints

Here we consider a fast graph-theoretical algorithm for the identification of linear constraint subsets $\mathcal{L}$ that satisfy (3.23) for a given multiplier variable $z_l$. We start by constructing a bipartite graph [64, 45, 51] $G_l$ where the set of nodes is partitioned into two disjoint subsets $\mathcal{N}^L$ and $\mathcal{N}^V$. We call these the “constraint” and “variable” nodes respectively. The former correspond to the set of linear constraints in [P] while the latter correspond to those problem variables which do not appear in any bilinear term (3.3) multiplied by $z_l$, i.e.:

$$\mathcal{N}_l^V \equiv \{k | \nexists i : (i, l, k) \in B\}$$

The set of edges $\mathcal{E}_l$ in graph $G_l$ is defined as:

$$\mathcal{E}_l \equiv \{(\lambda, k) | \lambda \in \mathcal{N}^L \wedge k \in \mathcal{N}_l^V \wedge a_{\lambda k} \neq 0\}$$

Thus, edge $(\lambda, k)$ exists if variable $z_k$ occurs in linear constraint $\lambda$.

Suppose, now, that $\mathcal{L}$ is a subset of the nodes $\mathcal{N}^L$. Let $\bar{K}(l, \mathcal{L})$ denote the subset of nodes $\mathcal{N}_l^V$ that are connected to nodes in $\mathcal{L}$, i.e.:

$$\bar{K}(l, \mathcal{L}) \equiv \{k | k \in \mathcal{N}_l^V \wedge \exists \lambda \in \mathcal{L} : (\lambda, k) \in \mathcal{E}_l\}$$

It can be verified that, in view of the definitions of sets $\mathcal{V}$, $\mathcal{N}_l^V$ and $\mathcal{E}_l$, the two sets $K(l, \mathcal{L})$ and $\bar{K}(l, \mathcal{L})$ are, in fact, identical. Consequently, a valid reduction constraint set (i.e. one that
satisfies criterion (3.23)) will correspond to a *dilation* in graph $G_l$, i.e. a subset $\mathcal{L}$ of, say, $m$ nodes $\mathcal{N}^L$ that are connected to fewer than $m$ distinct nodes $\mathcal{N}_l^V$. This provides a practical basis for the identification of valid reduction constraint sets using efficient existing algorithms for the determination of dilations in bipartite graphs.

### 3.3.3 Efficient identification of dilations

Dilations in bipartite graphs are closely related to the existence of *output set assignments* (OSA) in such graphs [87]. A subset $\mathcal{A}_l$ of the edges $\mathcal{E}_l$ in graph $G_l$ is an OSA if no node is incident to more than one edge in $\mathcal{A}_l$. A *complete* OSA is one in which each and every node $\lambda \in \mathcal{N}^L$ is incident to exactly one edge in it; in this case, $|\mathcal{A}_l| = |\mathcal{N}^L| = M$.

Given a non-complete OSA $\mathcal{A}_l$ of cardinality $m < M$, it may be possible to obtain one of cardinality $m + 1$ by identifying an unassigned node $\lambda \in \mathcal{N}^L$ (i.e. one that is not incident to any edge in $\mathcal{A}_l$) and tracing an *augmenting path* emanating from this node [29]. An augmenting path is a sequence of $2\alpha - 1$ edges of the form:

$$ (\lambda_1, k_1), (\lambda_2, k_1), (\lambda_2, k_2), \dots, (\lambda_\alpha, k_{\alpha-1}), (\lambda_\alpha, k_\alpha) $$

such that $(\lambda_\beta, k_\beta) \notin \mathcal{A}_l, \forall \beta = 1, \dots, \alpha$ and $(\lambda_\beta, k_{\beta-1}) \in \mathcal{A}_l, \forall \beta = 2, \dots, \alpha$. If such a path can be found, then an OSA of cardinality $m+1$ can be obtained from $\mathcal{A}_l$ simply by replacing the $\alpha-1$ edges $(\lambda_\beta, k_{\beta-1}), \beta = 2, \dots, \alpha$ in it by the $\alpha$ edges $(\lambda_\beta, k_\beta), \beta = 1, \dots, \alpha$.

Since the empty set is a valid (albeit trivial) OSA, the repeated tracing of augmenting paths can be used for the construction of OSAs of maximal cardinality.

A complete output set assignment for the graph $G_l$ exists if and only if the graph does not contain any dilations. Consequently, dilations can be identified by the application of efficient algorithms designed for the identification of output set assignments such as that proposed by [29]. A basic algorithm [29] for constructing an augmenting path emanating from a node $\lambda$ of graph $G_l$ is shown in pseudo-code notation in Figure 3.2. The procedure `AugmentPath` returns a boolean flag `PathFound` indicating whether or not an augmenting path has been found. `PathFound` is initially set to false.

The variable `Assigned[k]` contains the index of the constraint node to which variable node $k$ is currently assigned; a value of 0 indicates a variable node that has not yet been assigned. The boolean variable `ConstraintVisited[\lambda]` takes a value of true if constraint $\lambda$ has been visited as part of the attempted construction of the augmenting tree. Similarly, `VariableVisited[k]` marks
whether variable *k* has been visited. Initially, all elements of the array *Assigned* are set to 0, and all elements of *ConstraintVisited* and *VariableVisited* to false.

Line 2 of the algorithm simply marks the current constraint node $\lambda$ as “visited”. Then lines 3-7 aim to determine whether node $\lambda$ is adjacent to any node $k$ that is not yet assigned. If so, an augmenting path has been found (line 4); $k$ is immediately assigned to $\lambda$ (line 5), and the algorithm terminates.

If no immediate assignment can be found for constraint node $\lambda$ in lines 3-7, then this means that all variable nodes $k$ adjacent to $\lambda$ are already assigned to other constraint nodes. Lines 8-14 consider all such nodes $k$ provided they have not already been visited as part of the current search (i.e. *VariableVisited*[k] is still false, cf. line 8). If such a variable node $k$ exists, then it is immediately marked as “visited” (line 9). Node $k$ must be already assigned to a node *Assigned*[k] which is different to $\lambda$. At line 10, the algorithm attempts to determine an augmenting path starting from this other node via a recursive invocation to procedure *AugmentPath*. If such a path is found (line 11), then node $k$ is re-assigned to $\lambda$ (line 12), and no further search is necessary (line 13).

If neither of the searches in lines 3-7 and 8-15 manage to find an augmenting path, then the procedure *AugmentPath* terminates at line 16, having left the flag *PathFound* at its original value of false.

On termination of the algorithm, if *PathFound* is true, then the above recursive algorithm must have first executed an assignment of a yet-unassigned variable node *k* at line 5 at the deepest level of the recursion; this will have been followed by exactly one re-assignment (line 12) of an already assigned variable node *k* at each higher level of the recursion. This will have increased the total number of assignments by one, which is, of course, the intention of tracing an augmenting path.

On the other hand, if, on termination of the algorithm, *PathFound* is false, then it can be verified that (a) the number of constraint nodes marked as “visited” will exceed the number of variable nodes by exactly one, and (b) the algorithm will have visited each and every variable node adjacent to each visited constraint node. Consequently, the sets of constraint and variable nodes that have been visited form a dilation.

We are now in a position to formulate an algorithm that identifies the set of linear constraints $L_i$ that lead to valid reduction constraints when multiplied by a given variable $z_i$. The algorithm is shown, again in pseudo-code notation, in Figure 3.3. The initialization part (lines 2-6) sets the set $L_i$ to empty, constructs the bipartite graph $G_i$, and declares all variable nodes in it
1 PROCEDURE AugmentPath (λ, PathFound)
2 Set ConstraintVisited [λ] := TRUE
3 IF ∃ a node k : (λ, k) ∈ E AND Assigned [k] = 0 THEN
4 Set PathFound := TRUE
5 Set Assigned [k] := λ
6 RETURN
7 END IF
8 FOR every k : (λ, k) ∈ E AND NOT VariableVisited [k] DO
9 Set VariableVisited [k] := TRUE
10 AugmentPath (Assigned [k], PathFound)
11 IF PathFound THEN
12 Set Assigned [k] := λ
13 RETURN
14 END IF
15 END FOR
16 END AugmentPath

<caption>Figure 3.2: Algorithm for the construction of an augmenting path emanating from a linear constraint node λ in graph G_l.</caption>
1 PROCEDURE ValidReductionConstraints(l, Ll)
2 Set Ll := 
3 Construct bipartite graph Gl as described in 3.3.2
4 FOR each variable node k DO
5 Set Assigned[k] := 0
6 END FOR
7 FOR each constraint node λ DO
8 FOR each constraint node λ' DO
9 Set ConstraintVisited[λ'] := FALSE
10 END FOR
11 FOR each variable node k DO
12 Set VariableVisited[k] := FALSE
13 END FOR
14 Set PathFound := FALSE
15 AugmentPath(λ, PathFound)
16 IF NOT PathFound THEN
17 Set Ll := Ll ∪ {λ'|ConstraintVisited[λ'] = TRUE}
18 END IF
19 END FOR
20 END ValidReductionConstraints

<caption>Figure 3.3: Algorithm for identification of set of valid reduction constraint set L for variable zl.</caption>
as unassigned. The algorithm then considers each constraint node $\lambda$ in sequence (lines 7-19), attempting to construct an augmenting path emanating from it. Lines 8-14 perform the initializations necessary for the correct operation of procedure AugmentPath invoked on line 15 (cf. Figure 3.2). If an augmenting path is not found, then all constraints visited during the search for one are added to the set of constraints $\mathcal{L}_l$ (lines 16-18).

Procedure AugmentPath considers each edge in the bipartite graph at most once. Since procedure ValidReductionConstraints invokes AugmentPath once for each constraint node, the theoretical complexity of the algorithm is at most of the order (number of constraint nodes $\times$ number of edges). However, as observed by Duff (1981), in practice, the computational complexity of such OSA algorithms is usually nearer to (number of constraint nodes + number of edges). This results in an efficient algorithm that is practically applicable to large problems.

## 3.4 A detailed example

In order to illustrate the operation of the algorithm proposed in Section 3.3, we consider the following bilinear optimization problem in standard form:

$$ \left. \begin{array}{ll} \min_{z} & z_{24} \\ & z_{24} = c^T z \\ & Az' = b \\ & z_i = z_j z_k \quad \forall (i, j, k) \in B \\ & 0 \le z \le 10, \end{array} \right\} \qquad (3.24) $$

where:

$$ c^T = (0, 0, 0, 0, 0, 0, 1, 1, 3, 1, 2, 2, 1, 2, 1, -1, 2, -1, 4, 3, 6, 9, 1), $$

$$ z = (z_1, \dots, z_{23})^T, \quad z' = (z_1, \dots, z_6)^T, \quad b = (1, 2, -1, 1)^T, $$

$$ A = \begin{pmatrix} 1 & 2 & 1 & 1 \\ 2 & -1 & 1 & 3 \\ 1 & 6 & 2 & -3 \\ 2 & 1 & 3 & -1 \end{pmatrix} $$

and:

$$ B = \{(7, 1, 1), (8, 2, 2), (9, 4, 4), (10, 5, 5), (11, 6, 6), (12, 1, 4), \\ (13, 1, 5), (14, 1, 6), (15, 2, 4), (16, 2, 5), (17, 2, 6), (18, 3, 4), \\ (19, 3, 5), (20, 3, 6), (21, 4, 5), (22, 4, 6), (23, 5, 6)\} \qquad (3.25) $$
<figure><img src="image_6.png" /><figcaption>Figure 3.4: Generation of valid reduction constraints for multiplier variable $z_1$.</figcaption></figure>

The above is, effectively, a bilinear programming problem involving 6 variables $z_1, \dots, z_6$ which are, in turn, related via 4 linear constraints. The problem involves all bilinear combinations (including simple squares) of these 6 variables apart from $z_1z_2, z_1z_3, z_2z_3, z_3^2$. New variables $z_7, \dots, z_{23}$ have been introduced to represent these bilinear terms, and these are linearly combined to form the objective function represented by variable $z_{24}$.

In this problem, it is quite easy to see that only variables $z_1, \dots, z_6$ can possibly give rise to reduction constraints when multiplied by the linear constraints $Az' = b$. We consider each of these 6 variables in turn, applying to it procedure ValidReductionConstraints of Figure 3.3.

### 3.4.1 Valid reduction constraints by multiplication with variable $z_1$

The bipartite graph $G_1$ constructed as described in Section 3.3.2 for variable $z_1$ is shown in Figure 3.4a. Nodes $c_1, \dots, c_4$ correspond to the four linear constraints. Edges correspond to the creation of new bilinear terms if a particular constraint were to be multiplied by $z_1$. For example, edges $(c_1, z_1)$ and $(c_1, z_3)$ exist because, if constraint $c_1$ were to be multiplied by $z_1$, then this would give rise to new terms $z_1^2$ and $z_1z_3$ respectively.

As there are initially no assignments of variable nodes to constraint nodes, all edges are shown as dotted lines in Figure 3.4a. We now consider each constraint node in turn (lines 7-19 of Figure 3.3), attempting to construct an augmenting path emanating from it (line 15).

Starting from constraint node $c_1$, an augmenting path is found immediately to its adjacent variable node $z_3$. This is assigned to $c_1$ (cf. line 5 of Figure 3.2); this assignment is indicated as a solid line in Figure 3.4b.
Considering constraint $c_2$ also results in an augmenting path being found immediately, this time assigning $z_2$ to $c_2$. The resulting output set assignment is also shown as a solid line in Figure 3.4b.

So far, we have not identified any dilations. However, if we now attempt to construct an augmenting path emanating from constraint node $c_3$, we note that $c_3$ is adjacent to only one variable node, namely $z_2$ which is already assigned to node $c_2$. Consequently, the first search loop (lines 3-7) of the AugmentPath procedure fails to detect an augmenting path, and we therefore have to enter the second loop (lines 8-15). We mark $z_2$ as visited (line 9), and then recursively invoke AugmentPath (line 10) asking it to search for an augmenting path emanating from the constraint node currently assigned to $z_2$, namely $c_2$. However, $c_2$ is not adjacent to any variables other than $z_2$ which has already been visited; consequently, the recursive invocation of AugmentPath returns PathFound false. As there are no more variable nodes adjacent to $c_3$, the first invocation of AugmentPath also returns PathFound false. Consequently, line 17 of Figure 3.3 adds the two constraint nodes visited during this unsuccessful search for an augmenting path, i.e. $c_2$ and $c_3$ to the set of linear constraints $\mathcal{L}_1$ to be multiplied by $z_1$.

Finally, constraint node $c_4$ is completely isolated in $\mathcal{G}_1$, i.e. there are no edges incident to it. Thus, the corresponding invocation of procedure AugmentPath at line 15 of Figure 3.3 returns immediately with PathFound false. Consequently, $c_4$ is also added to set $\mathcal{L}_1$.

As there are no more constraint nodes to be considered, we conclude that $\mathcal{L}_1 = \{c_2, c_3, c_4\}$. Essentially, the algorithm has identified automatically that multiplying constraints $c_2$:

$$2z_1 - z_2 + z_4 + 3z_6 = 2$$

and $c_3$:

$$z_2 + 6z_4 + 2z_5 - 3z_6 = -1$$

by $z_1$ creates only one new bilinear term, namely $z_1z_2$, beyond those that already exist in the problem. Consequently, if we define a new variable $z_{25} \equiv z_1z_2$, we end up with two linear constraints:

$$2z_7 - z_{25} + z_{12} + 3z_{14} = 2z_1$$

and

$$z_{25} + 6z_{12} + 2z_{13} - 3z_{14} = -z_1$$

respectively that can be used to eliminate two bilinear terms from the problem.

Moreover, the algorithm has also detected that multiplying constraint $c_4$:

$$2z_1 + z_4 + 3z_5 = 1$$
<figure><img src="image_7.png" /><figcaption>Figure 3.5: Generation of valid reduction constraints for multiplier variable $z_2$.</figcaption></figure>

by $z_1$ does not generate any new bilinear terms. The new linear constraint:

$$2z_7 + z_{12} + 3z_{13} = z_1$$

obtained by this multiplication can be used to eliminate one more of the problem's bilinear terms. Overall, the number of bilinear terms in the problem can be reduced by two via the use of valid reduction constraints.

### 3.4.2 Valid reduction constraints by multiplication with variable $z_2$

The bipartite graph $G_2$ is shown in Figure 3.5a. An augmenting path emanating from constraint node $c_1$ can be found immediately, resulting in the assignment of $z_1$ to it (see solid line in Figure 3.5b.

We now proceed to consider constraint node $c_2$. The only variable adjacent to it is $z_1$ but this is already assigned. Thus, no augmenting path is found in lines 3-7 of Figure 3.2 and we are forced to enter the second search loop (lines 8-15). This considers variable $z_1$ again and recursively invokes procedure AugmentPath (line 10) asking it to search for an augmenting path emanating from the constraint node currently assigned to $z_1$, namely $c_1$. This recursive invocation immediately identifies an unassigned variable $z_3$ as shown in Figure 3.5c. Variable $z_3$ is assigned to $c_1$ while $z_1$ is re-assigned to $c_2$, resulting in the output set assignment shown in Figure 3.5d.

Constraint node $c_3$ is completely isolated. Thus, the corresponding invocation of procedure AugmentPath at line 15 of Figure 3.3 returns immediately with *PathFound* false. Consequently, $c_3$ is added to the set of linear constraints $\mathcal{L}_2$ (line 17).

Finally, attempting to construct an augmenting path emanating from constraint node $c_4$ visits first variable node $z_1$ and from there, the constraint node $c_2$ to which $z_1$ is currently assigned.
<figure><img src="image_8.png" /><figcaption>Figure 3.6: Generation of valid reduction constraints for multiplier variable $z_3$.</figcaption></figure>

No further progress is possible, and therefore `AugmentPath` returns with `PathFound` false, having visited both $c_4$ and $c_2$. Consequently, these are also added to the set $\mathcal{L}$.

We conclude that $\mathcal{L}_2 = \{c_2, c_3, c_4\}$. In this case, the algorithm has detected that $c_3$ can be multiplied by $z_2$ without generating any new bilinear terms. Also, the multiplication of $c_2$ and $c_4$ by $z_2$ leads to the creation of only one new bilinear term, $z_1z_2$ and two linear constraints.

### 3.4.3 Valid reduction constraints by multiplication with variable $z_3$

The bipartite graph $G_3$ is shown in Figure 3.6a. In this case, augmenting paths emanating from both nodes $c_1$ and $c_2$ can be found immediately, resulting in the assignment of $z_1$ to $c_1$, and $z_2$ to $c_2$ respectively (see solid lines in Figure 3.6b).

The application of procedure `AugmentPath` to constraint node $c_3$ results in tracing an augmenting path $c_3 \to z_2 \to c_2 \to z_1 \to c_1 \to z_3$, as shown in Figure 3.6c. We therefore assign $z_3$ to $c_1$, and re-assign $z_1$ to $c_2$ and $z_2$ to $c_3$. The output set assignment up to this point is shown in Figure 3.6d.

Finally, attempting to locate an augmenting path emanating from $c_4$ visits nodes $z_1$, $c_2$, $z_2$, $c_3$. However, no further progress is possible in this case. Procedure `AugmentPath` returns `PathFound` false. Consequently, all three constraint nodes visited are added to set $\mathcal{L}_3$.

We conclude that $\mathcal{L}_3 = \{c_2, c_3, c_4\}$. Here the algorithm has detected that multiplying these three constraints by $z_3$ results in the creation of only two new bilinear terms, namely $z_1z_3$ and $z_2z_3$.
### 3.4.4 Valid reduction constraints by multiplication with variables $z_4, \dots, z_6$

The bipartite graphs $G_l$ for $l = 4, 5, 6$ are very simple as they contain no variable nodes (i.e. $N_l^V = \emptyset$). Consequently, the invocation of procedure AugmentPath for each and every constraint node $c_\lambda$, $\lambda = 1, \dots, 4$ immediately returns with *PathFound* false. This results in the constraint sets:

$$ \mathcal{L}_l = \{c_1, c_2, c_3, c_4\}, \quad l = 4, 5, 6 $$

This, of course, is a simple consequence of the fact that this problem already contains all possible bilinear terms involving $z_4$ or $z_5$ or $z_6$. Consequently, multiplying any linear constraint by one of these variables does not result in any new bilinear terms and can, therefore, be used to generate a valid reduction constraint.

### 3.4.5 The reformulated NLP

In summary, the operations described in sections 3.4.1-3.4.4 above result in the following set of 21 reduction constraints:
From multiplication by $z_1$ :

$$z_1 \times c_2 \Rightarrow 2z_7 - z_{25} + z_{12} + 3z_{14} - 2z_1 = 0$$

$$z_1 \times c_3 \Rightarrow z_{25} + 6z_{12} + 2z_{13} - 3z_{14} + z_1 = 0$$

$$z_1 \times c_4 \Rightarrow 2z_7 + z_{12} + 3z_{13} - z_1 = 0$$

From multiplication by $z_2$ :

$$z_2 \times c_2 \Rightarrow 2z_{25} - z_8 + z_{15} + 3z_{17} - 2z_2 = 0$$

$$z_2 \times c_3 \Rightarrow z_8 + 6z_{15} + 2z_{16} - 3z_{17} + z_2 = 0$$

$$z_2 \times c_4 \Rightarrow 2z_{25} + z_{15} + 3z_{16} - z_2 = 0$$

From multiplication by $z_3$ :

$$z_3 \times c_2 \Rightarrow 2z_{26} - z_{27} + z_{18} + 3z_{20} - 2z_3 = 0$$

$$z_3 \times c_3 \Rightarrow z_{27} + 6z_{18} + 2z_{19} - 3z_{20} + z_3 = 0$$

$$z_3 \times c_4 \Rightarrow 2z_{26} + z_{18} + 3z_{19} - z_3 = 0$$

From multiplication by $z_4$ :

$$z_4 \times c_1 \Rightarrow z_{12} + 2z_{18} + z_{21} + z_{22} - z_4 = 0$$

$$z_4 \times c_2 \Rightarrow 2z_{12} - z_{15} + z_9 + 3z_{22} - 2z_4 = 0$$

$$z_4 \times c_3 \Rightarrow z_{15} + 6z_9 + 2z_{21} - 3z_{22} + z_4 = 0$$

$$z_4 \times c_4 \Rightarrow 2z_{12} + z_9 + 3z_{21} - z_4 = 0$$

From multiplication by $z_5$ :

$$z_5 \times c_1 \Rightarrow z_{13} + 2z_{19} + z_{10} + z_{23} - z_5 = 0$$

$$z_5 \times c_2 \Rightarrow 2z_{13} - z_{16} + z_{21} + 3z_{23} - 2z_5 = 0$$

$$z_5 \times c_3 \Rightarrow z_{16} + 6z_{21} + 2z_{10} - 3z_{23} + z_5 = 0$$

$$z_5 \times c_4 \Rightarrow 2z_{13} + z_{21} + 3z_{10} - z_5 = 0$$

From multiplication by $z_6$ :

$$z_6 \times c_1 \Rightarrow z_{14} + 2z_{20} + z_{23} + z_{11} - z_6 = 0$$

$$z_6 \times c_2 \Rightarrow 2z_{14} - z_{17} + z_{22} + 3z_{11} - 2z_6 = 0$$

$$z_6 \times c_3 \Rightarrow z_{17} + 6z_{22} + 2z_{23} - 3z_{11} + z_6 = 0$$

$$z_6 \times c_4 \Rightarrow 2z_{14} + z_{22} + 3z_{13} - z_6 = 0$$
Three new bilinear terms were introduced:

$$z_{25} = z_1 z_2$$

$$z_{26} = z_1 z_3$$

$$z_{27} = z_2 z_3$$

augmenting the original set B (cf. equation (3.3)) with the triplets (25, 1, 2), (26, 1, 3) and (27, 2, 3).

In order to determine which bilinear terms may be replaced by these linear constraints, we write the latter in matrix form:

$$Bz'' + Cz' = 0. \qquad (3.26)$$

where $z' = (z_1, \dots, z_6)^T \in \mathbb{R}^6$ and $z'' = (z_7, \dots, z_{23}, z_{25}, \dots, z_{27})^T \in \mathbb{R}^{20}$, $B: \mathbb{R}^{20} \to \mathbb{R}^{21}$ and $C: \mathbb{R}^6 \to \mathbb{R}^{21}$. By performing Gaussian elimination with row pivoting on B (and replicating the same row operations on C) we obtain a system of the form:

$$\begin{pmatrix} U & \tilde{B} \\ 0 & 0 \end{pmatrix} \begin{pmatrix} \xi \\ \zeta \end{pmatrix} = \begin{pmatrix} \tilde{C} \\ 0 \end{pmatrix} z'$$

where $\xi = (z_7, \dots, z_{23}, z_{25})^T \in \mathbb{R}^{18}$, $\zeta = (z_{26}, z_{27})^T \in \mathbb{R}^2$, $U: \mathbb{R}^{18} \to \mathbb{R}^{18}$ is a nonsingular, upper-triangular matrix, $\tilde{B}: \mathbb{R}^2 \to \mathbb{R}^{18}$ and $\tilde{C}: \mathbb{R}^6 \to \mathbb{R}^{18}$. Thus, the reduction constraints determine the variables $\xi$ in terms of $\zeta$ and $z'$ via the solution of $U\xi = \tilde{C}z' - \tilde{B}\zeta$, and consequently, the corresponding triplets $(i, j, k)$ can be deleted from the set B in equation (3.25). Overall, the reformulation of the original problem (3.24) involves the introduction of the reduction constraints (3.26) and a much smaller set of just two triplets, namely $B = \{(26, 1, 3), (27, 2, 3)\}$.

## 3.5 Computational results

We have chosen a selection of bilinear test problems relative to pooling and blending problems to test the efficiency of our algorithm. Pooling and blending involves the determination of optimal amounts of different raw materials that need to be mixed to produce required amounts of end-products with desired properties. Such problems occur frequently in the petrochemical industry and are well known to exhibit multiple local minima.
Here, we use the general blending problem formulation found in [1]:

$$ \min_{f,q,x} \sum_{j=1}^{p} \sum_{i=1}^{n_j} c_{ij}f_{ij} - \sum_{k=1}^{r} d_k \sum_{j=1}^{p} x_{jk} \quad (3.27) $$

$$ \sum_{i=1}^{n_j} f_{ij} - \sum_{k=1}^{r} x_{jk} = 0, \quad \forall j \le p \quad (3.28) $$

$$ q_{jw} \sum_{k=1}^{r} x_{jk} - \sum_{i=1}^{n_j} \lambda_{ijw} f_{ij} = 0, \quad \forall j \le p, \forall w \le l \quad (3.29) $$

$$ \sum_{j=1}^{p} x_{jk} \le S_k, \quad \forall k \le r \quad (3.30) $$

$$ \sum_{j=1}^{p} q_{jw}x_{jk} - Z_{kw} \sum_{j=1}^{p} x_{jk} \le 0, \quad \forall k \le r, \forall w \le l \quad (3.31) $$

$$ f^L \le f \le f^U, \quad q^L \le q \le q^U, \quad x^L \le x \le x^U, \quad (3.32) $$

where $f_{ij}$ is the flow of input stream $i$ into pool $j$, $x_{jk}$ is the total flow from pool $j$ to product $k$ and $q_{jw}$ is the $w$-th quality of pool $j$; $p$ is the number of pools, $r$ the number of products, $l$ the number of qualities, $n_j$ the number of streams; $c_{ij}$, $d_k$, $S_k$, $Z_{kw}$, $\lambda_{ijw}$ are given parameters.

When the blending problem (3.27)-(3.32) is reformulated to the standard form [P], new variables $u_j$ are created to replace the term $\sum_{k=1}^{r} x_{jk}$ in constraint set (3.29). Thus, we introduce the linear constraints:

$$ \sum_{k=1}^{r} x_{jk} - u_j = 0, \quad \forall j \le p \quad (3.33) $$

and re-write (3.29) as

$$ t_{jw} - \sum_{i=1}^{n_j} \lambda_{ijw} f_{ij} = 0, \quad \forall j \le p, \forall w \le l $$

where we have introduced new variables $t_{jw}$ derived via the bilinear constraints:

$$ t_{jw} = q_{jw}u_j, \quad \forall j \le p, \forall w \le l $$

More new variables $z_{jwk}$ are created in the standard form to replace the bilinear terms $q_{jw}x_{jk}$ in constraint set (3.31):

$$ z_{jwk} = q_{jw}x_{jk}, \quad \forall j \le p, \forall w \le l, \forall k \le r $$

which allows (3.31) to be re-written in linear form as:

$$ \sum_{j=1}^{p} z_{jwk} - Z_{kw} \sum_{j=1}^{p} x_{jk} = 0, \quad \forall k \le r, \forall w \le l. $$
The standard form reformulation of (3.27)-(3.32) is shown below:

$$ \min_{f,q,x} \sum_{j=1}^{p} \sum_{i=1}^{n_j} c_{ij} f_{ij} - \sum_{k=1}^{r} d_k \sum_{j=1}^{p} x_{jk} $$

subject to linear constraints:

$$ \sum_{i=1}^{n_j} f_{ij} - \sum_{k=1}^{r} x_{jk} = 0 \quad \forall j \le p $$

$$ t_{jw} - \sum_{i=1}^{n_j} \lambda_{ijw} f_{ij} = 0 \quad \forall j \le p, \forall w \le l $$

$$ \sum_{k=1}^{r} x_{jk} - u_j = 0 \quad \forall j \le p $$

$$ \sum_{j=1}^{p} x_{jk} \le S_k \quad \forall k \le r $$

$$ \sum_{j=1}^{p} z_{jwk} - Z_{kw} \sum_{j=1}^{p} x_{jk} = 0 \quad \forall k \le r, \forall w \le l $$

with bilinear terms:

$$ B = \{ (t_{jw}, q_{jw}, u_j), (z_{jwk}, q_{jw}, x_{jk}) \mid j \le p, w \le l, k \le r \} $$

and bounds:

$$ f^L \le f \le f^U, \quad q^L \le q \le q^U, \quad x^L \le x \le x^U. $$

We now apply the algorithm of this chapter to the above standard form. This results in a set of reduction constraints derived by multiplying constraints (3.33) by the quality variables $q_{jw}$:

$$ \begin{align*} \forall j \le p \forall w \le l \quad & q_{jw} \left( \sum_{k=1}^{r} x_{jk} - u_j \right) = \\ & \sum_{k=1}^{r} q_{jw} x_{jk} - q_{jw} u_j = \\ & \sum_{j=1}^{r} z_{jwk} - t_{jw} = 0. \end{align*} $$

These constraints can be used to eliminate the bilinear constraints $t_{jw} = q_{jw}u_j$.

Another way to interpret what is happening in this particular case is that the constraint sets (3.29) and (3.31) (the bilinear constraints in the general blending problem formulation above) define more bilinear products than is really necessary: if we were to re-write the first term on the left hand side of (3.29) as $\sum_{k=1}^{r} q_{jw}x_{jk}$, we would not need to create all the variables $t_{jw}$.
<table><caption>Table 3.1: Test problem statistics.</caption><thead><tr><th></th><th>n<sub>org</sub></th><th>m<sub>org</sub></th><th>n<sub>std</sub></th><th>m<sup>l</sup><sub>std</sub></th><th>|BT|</th><th>|RC|</th><th>|NB|</th></tr></thead><tbody><tr><td>Haverly 1</td><td>9</td><td>8</td><td>18</td><td>11</td><td>6</td><td>2</td><td>0</td></tr><tr><td>Haverly 2</td><td>9</td><td>8</td><td>18</td><td>11</td><td>6</td><td>2</td><td>0</td></tr><tr><td>Haverly 3</td><td>9</td><td>8</td><td>18</td><td>11</td><td>6</td><td>2</td><td>0</td></tr><tr><td>Foulds 2</td><td>26</td><td>16</td><td>51</td><td>21</td><td>20</td><td>4</td><td>0</td></tr><tr><td>Foulds 3</td><td>168</td><td>48</td><td>313</td><td>57</td><td>136</td><td>8</td><td>0</td></tr><tr><td>Foulds 4</td><td>168</td><td>48</td><td>313</td><td>57</td><td>136</td><td>8</td><td>0</td></tr><tr><td>Foulds 5</td><td>100</td><td>40</td><td>173</td><td>45</td><td>68</td><td>4</td><td>0</td></tr><tr><td>Ben-Tal 4</td><td>10</td><td>8</td><td>19</td><td>11</td><td>6</td><td>2</td><td>0</td></tr><tr><td>Ben-Tal 5</td><td>41</td><td>27</td><td>94</td><td>32</td><td>48</td><td>8</td><td>0</td></tr><tr><td>example 1</td><td>21</td><td>30</td><td>64</td><td>33</td><td>40</td><td>8</td><td>0</td></tr><tr><td>example 2</td><td>25</td><td>42</td><td>88</td><td>45</td><td>60</td><td>12</td><td>0</td></tr><tr><td>example 3</td><td>38</td><td>49</td><td>132</td><td>53</td><td>90</td><td>18</td><td>0</td></tr><tr><td>example 4</td><td>26</td><td>35</td><td>77</td><td>38</td><td>48</td><td>8</td><td>0</td></tr></tbody></table>

Yet another way of saying this is that distributing products over sums is advantageous. This is in accordance with the considerations found in [130], p. 73. However, it is important to note that, here, this reformulation is determined automatically by a generic algorithm.

Table 3.1 summarizes the main characteristics of the test problems. Here, $n_{org}$, $m_{org}$ are respectively the number of variables and constraints in the original problem formulation, $n_{std}$ is the number of variables in the standard form $[P]$, $m_{std}^l$ is the number of linear constraints in the standard form, $|RC|$ is the number of reduction constraints created, $|BT|$ is the number of bilinear constraints in the standard form before reduction constraint creation and $|NB|$ is the number of new bilinear terms created during the reduction constraint creation procedure. We note that $|NB| = 0$ in all test problems considered.

We now proceed to consider how the addition of reduction constraints affects the tightness of the convex relaxation of the NLP in the context of its solution using deterministic global optimization algorithms. Table 3.2 compares the number of nodes needed to solve the problems of Table 3.1 using eight different codes. The first six are codes described in the literature:
1 = [40] (Foulds, 1992)

2 = [138] (Visweswaran, 1993)

3 = [17] (Ben-Tal, 1994)

4 = [139] (Visweswaran, 1996)

5 = [1] (Adhya, 1999)

6 = [126] (Tawarmalani, 1999)

and their performance is shown in Table 3.2, taken from the corresponding papers. Codes 7 and 8 both correspond to a rather basic implementation of sBB based on the algorithm proposed in [116]. The only difference between them is that code 7 does not incorporate the reduction constraints while code 8 does.

<table><caption>Table 3.2: Numerical results.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th></tr></thead><tbody><tr><td>Haverly 1</td><td>5</td><td>7</td><td>3</td><td>12</td><td>3</td><td>3</td><td>31</td><td>1</td></tr><tr><td>Haverly 2</td><td></td><td>19</td><td>3</td><td>12</td><td>9</td><td>9</td><td>43</td><td>7</td></tr><tr><td>Haverly 3</td><td></td><td></td><td>3</td><td>14</td><td>5</td><td>3</td><td>39</td><td>7</td></tr><tr><td>Foulds 2</td><td>9</td><td></td><td></td><td></td><td>1</td><td>1</td><td>131</td><td>7</td></tr><tr><td>Foulds 3</td><td>1</td><td></td><td></td><td></td><td>1</td><td>1</td><td>&gt; 20,000</td><td>1</td></tr><tr><td>Foulds 4</td><td>25</td><td></td><td></td><td></td><td>1</td><td>1</td><td>&gt; 20,000</td><td>1</td></tr><tr><td>Foulds 5</td><td>125</td><td></td><td></td><td></td><td>1</td><td>1</td><td>&gt; 20,000</td><td>1</td></tr><tr><td>Ben-Tal 4</td><td></td><td>47</td><td>25</td><td>7</td><td>3</td><td>3</td><td>101</td><td>1</td></tr><tr><td>Ben-Tal 5</td><td></td><td>42</td><td>283</td><td>41</td><td>1</td><td>1</td><td>&gt; 200,000</td><td>1</td></tr><tr><td>example 1</td><td></td><td></td><td></td><td></td><td>6174</td><td>1869</td><td>11245</td><td>5445</td></tr><tr><td>example 2</td><td></td><td></td><td></td><td></td><td>10743</td><td>2087</td><td>83051</td><td>11049</td></tr><tr><td>example 3</td><td></td><td></td><td></td><td></td><td>79944</td><td>7369</td><td>&gt; 200,000</td><td>7565</td></tr><tr><td>example 4</td><td></td><td></td><td></td><td></td><td>1980</td><td>157</td><td>2887</td><td>1467</td></tr></tbody></table>

As can be seen from the results in column 7 of Table 3.2, the lack of sophisticated algorithmic features regarding, for example, the choice of branching strategy, makes the performance of our basic implementation of the sBB algorithm significantly worse than that of earlier codes presented in the literature. However, the mere addition of the reduction constraints to this basic
code improves its performance dramatically, to the extent that it outperforms all but one of the other codes (cf. column 8 of Table 3.2).

It is worth noting that the BARON software (columns 5 and 6 of Table 3.2) implements a Branch-and-Bound method that employs various methods to generate valid cuts; however, none of these is currently equivalent to our method. Our results are in the same order of magnitude as those of BARON (the version published in 1999, fifth column of Table 3.2). However, at the end of 1999 a new version of BARON was described [126] that used an improved branching scheme that cut the number of iterations by an order of magnitude (sixth column). With the notable exception of Example 3, our code cannot generally attain these levels. Nevertheless, our reformulation is complementary to the techniques implemented in BARON and their combination could lead to a further improvement in performance.

## 3.6 Generalization of the graph-theoretical algorithm

The graph-theoretical algorithm described in Section 3.3 for the identification of useful reduction constraints has some limitations. These arise from the fact that potential multiplications are considered separately for each variable, which may result in some useful multiplications being missed. To understand this point, consider the following very simple example:

$$ \min \quad \begin{array}{l} x_1^2 + x_2^2 \\ x_1 + x_2 \\ 0 \le x_1, x_2 \le 10. \end{array} \quad \left. \begin{array}{l} = 1 \\ \end{array} \right\} \qquad (3.34) $$

On multiplying the linear constraint $x_1 + x_2 = 1$ by $x_1$ we obtain $x_1^2 + x_1x_2 = x_1$, which introduces one new bilinear term $x_1x_2$; thus, this multiplication would not appear to bring any benefit. Similarly, on multiplying the same linear constraint by $x_2$, we would get $x_1x_2 + x_2^2 = x_2$ and thus, again, one new bilinear term $x_1x_2$; hence this multiplication would not be considered to be beneficial either. Consequently, the algorithm of Section 3.3 applied to this system will not create any reduction constraint.

However, considering the combined effect of the two multiplications, we note that they produce two new linearly independent reduction constraints while introducing only one new bilinear term (namely, $x_1x_2$). These two reduction constraints can be used to eliminate two bilinear
terms from the problem, e.g. $x_1^2$ and $x_2^2$, leading to the problem reformulation:

$$
\begin{align*}
\min \quad & w_1 + w_3 \\
& x_1 + x_2 = 1 \\
& w_1 + w_2 = x_1 \\
& w_2 + w_3 = x_2 \\
& w_2 = x_1 x_2 \\
& 0 \le x_1, x_2 \le 10 \\
& 0 \le w_i \le 100, \quad i = 1, 2, 3.
\end{align*}
$$

Essentially, the algorithm of Section 3.3 correctly identifies that multiplying the linear constraint by either $x_1$ or $x_2$ results in a bilinear term $x_1x_2$ that did not occur in the original problem. What it misses is the fact that it is the *same* bilinear term in both cases. This is unfortunate as such reformulation may be very beneficial. For instance, we have found that the numerical solution of the example above in its original form with a simple sBB algorithm requires the examination of 255 nodes, while the reformulated one can be solved in a single node.

This motivates an extension to the algorithm of Section 3.3 to consider simultaneously mul-
tiplication of linear constraints by all system variables. Instead of creating one bipartite graph
for each multiplier variable, we consider one unified bipartite graph comprising two disjoint
sets of nodes:

* The $\rho$-nodes which comprise $m \times n$ nodes $\rho_{ij}$ representing the potential multiplication of constraint $i$ by variable $x_j$.

* The $\sigma$-nodes which comprise $n^2$ nodes $\sigma_{hk}$ representing the bilinear term $x_h x_k$.

An edge connecting node $\rho_{ij}$ and $\sigma_{hk}$ exists if the multiplication of constraint $i$ by variable $x_j$ would give rise to a new bilinear term $x_h x_k$; obviously either $h = j$ or $k = j$ holds.

Having created the bipartite graph, we attempt to trace an augmenting path algorithm similar
to that described in Figure 3.2 emanating from each node $\rho_{ij}$. If no such path is found, then
we must have identified a dilation involving $\nu$ nodes of type $\rho$ and $\nu - 1$ nodes of type $\sigma$. This
implies that the variable-constraint multiplications corresponding to these $\rho$ nodes will result in
$\nu$ new constraints but only $\nu - 1$ new bilinear terms – which is exactly what we are trying to
establish.

We apply this generalized algorithm to the simple example problem (3.34). Assuming that
the constraint $x_1 + x_2 = 1$ is labelled with index $i = 1$, the unified bipartite graph is shown in
Fig. 3.7. Tracing an augmenting path from node $\rho_{11}$ is, resulting in node $\sigma_{12}$ being assigned to it. However, no augmenting path emanating from node $\rho_{12}$ can be found. Instead, we identify a dilation comprising nodes $\rho_{11}, \sigma_{12}$ and $\rho_{12}$, i.e. in this case, $\nu = 2$. This simply implies that the reformulated problem should involve multiples of the linear constraint by both variables, which introduces a single bilinear term.

<figure><img src="image_9.png" /><figcaption>Figure 3.7: Unified bipartite graph for problem (3.34) in the generalized algorithm.</figcaption></figure>

As has been mentioned, the worst-case computational complexity of the dilation-finding algorithm of Section 3.3.3 is proportional to the product of the number of nodes from which the augmenting paths emanate and the number of edges; on the other hand, the average-case complexity is nearer to the sum of these two numbers. The original procedure described in Section 3.3 was applied to a bipartite graph with $m$ linear constraint nodes; in principle, each of these nodes could be connected with each and every one of the $n$ variable nodes; consequently, the number of edges is bounded from above by $mn$. Therefore, the worst-case and average complexities of this procedure are $O(m^2n)$ and $O(m + mn)$ respectively. Of course, the procedure has to be applied separately for each candidate multiplier variable; therefore, the corresponding total complexities are $O(m^2n^2)$ and $O(mn + mn^2)$ respectively. On the other hand, the procedure described in this section is applied to a bipartite graph with $mn$ $\rho$-nodes and $n^2$ $\sigma$-nodes. Each $\rho$-node is potentially connected with up to $n$ $\sigma$-nodes, and therefore the number of edges is bounded from above by $mn^2$ edges. Consequently, the worst-case and average complexities are $O(m^2n^3)$ and $O(mn + mn^2)$ respectively. In conclusion, the procedure of this section has worse worst-case complexity than that of Section 3.3, but similar average complexity.

On the other hand, the memory requirements of the two procedures are very different. With the original algorithm, a graph consisting of $m + n$ nodes and up to $mn$ vertices needs to be stored in memory at any one time, whereas the unified bipartite graph of this section will have
$mn + n^2$ nodes and up to $mn^2$ edges. For extremely large problems, these requirements may be excessive and special attention may have to be paid to the implementation of the algorithm (one could create the graph “on-the-fly” as the algorithm progresses, deriving the nodes and edges from the problem data at each step).

## 3.7 Concluding remarks

The work presented in this chapter is based on the fact that geometrical intersections of hyperplanes and nonlinear hypersurfaces corresponding to bilinearities may embed a higher degree of linearity than what is apparent by mere inspection of the defining equations. We have shown that it is possible to exploit this fact so as to reformulate an NLP involving such equality constraints to a form with fewer bilinearities and more linear constraints.

The basic idea of the reformulation is to multiply subsets of the NLP’s linear constraints by one of the system variables $x$. This creates new linear “reduction” constraints expressed in terms of variables $w$, each one of which corresponds to a bilinear product, viz. $w_{jk} \equiv x_j x_k$. In general, some of these $w$ variables will already exist in the original NLP while others are new variables introduced by the multiplication.

In principle, any of the original linear constraints in the NLP can be multiplied by any variable; this is the basis of the RLT procedure proposed in [111, 108]. However, here we have focussed on identifying multiplications which result in “valid” sets of reduction constraints, i.e. sets in which the number of constraints exceeds the number of new variables $w$ introduced by the multiplication.

Valid reduction constraints do not affect the feasible region of the original NLP, but they do reduce the feasible region of its convex relaxation, thereby rendering it tighter. To see this, consider the multiplication of a subset of $\nu$ of the original linear constraints $Ax = b$ by a single variable $x_k$. This subset can be written in the form:

$$A'x' + A''x'' = b' \qquad (3.35)$$

where the partitioning is such that the bilinear terms $w' \equiv x_k x'$ already occur in the NLP, while $w'' \equiv x_k x''$ do not. Now, if the above is to form a valid reduction constraint set, we must have $x'' \in R^{\nu'}$ where $\nu' < \nu$. We multiply the above constraint by $x_k$ to obtain:

$$A'w' + A''w'' = x_k b' \qquad (3.36)$$
If we now apply Gaussian elimination with row pivoting on $A''$, we can bring the above system to the form:

$$ \begin{pmatrix} \tilde{A}' \\ \bar{A}' \end{pmatrix} w' + \begin{pmatrix} \tilde{A}'' \\ 0 \end{pmatrix} w'' = x_k \begin{pmatrix} \tilde{b}' \\ \bar{b}' \end{pmatrix} \qquad (3.37) $$

We now note that the bottom block row of the above equation, $\tilde{A}'w' = \bar{b}'$ constrains the admissible values of the variables $w'$ which already existed in the original NLP. Although these constraints are redundant with respect to the NLP itself, they are non-redundant with respect to any relaxation which does not enforce the equality $w' = x_kx'$ exactly.

The above argument helps explain the beneficial effects of the proposed reformulation. On one hand, the (linear) valid reduction constraints always tighten the convex relaxation of the NLP. On the other hand, the elimination of some of the bilinear terms reduces the size of the convex relaxation (e.g. by obviating the need for McCormick relaxations for these terms). Overall, we have a relaxation that is both tighter and smaller, something that is not always the case with earlier reformulation methods such as RLT.

Finally, it is worth noting that our reformulation relies on the use of graph theoretical algorithms for the identification of valid reduction constraint sets. These algorithms have the advantage of being quite fast even when applied to relatively large systems. On the other hand, they may fail to identify some valid reduction constraint sets. This may occur in cases in which the matrix $A''$ in equation (3.35) is numerically singular but structurally non-singular.
# Chapter 4

## A convex relaxation for monomials of odd degree

One of the most effective techniques for the solution of nonlinear programming problems (NLPs) to global optimality is the spatial Branch-and-Bound (sBB) method. This requires the computation of a lower bound to the solution, usually obtained by solving a convex relaxation of the original NLP. The formation and tightness of such a convex relaxation are critical issues in any sBB implementation.

As will be shown in sections 5.2.2.1 and 5.2.2.2, it is possible to form a convex relaxation of any NLP by isolating the nonconvex terms and replacing them with their convex relaxation. Tight convex underestimators are already available for many types of nonconvex term, including bilinear and trilinear products, linear fractional terms, and concave and convex univariate functions. However, terms which are piecewise concave and convex are not explicitly catered for. A frequently occurring example of such a term is $x^{2k+1}$, where $k \in \mathbb{N}$ and the range of $x$ includes zero. A detailed analysis of the conditions required for concavity and convexity of polynomial functions has been given in [75]; however, the results obtained therein only apply to the convex underestimation of multivariate polynomials with positive variable values. For monomials of odd degree, where the variable ranges over both negative and positive values, no special convex envelopes have been proposed in the literature, and one therefore has to rely either on generic convex relaxations such as those given by Floudas and co-workers (see [13, 7]) or on reformulation in terms of other types of terms for which convex relaxations are available.
In this chapter, we propose convex/concave nonlinear envelopes for odd power terms of the form $x^{2k+1}$ ($k \in \mathbb{N}$), where $x \in [a, b]$ and $a < 0 < b$. These envelopes are continuous and differentiable everywhere in $[a, b]$. We also derive tight linear relaxations. We compare both of these relaxations with relaxations for the same terms derived using other methods. We shall, with a slight abuse of notation, speak about “envelope” to mean the region enclosed between the convex and the concave envelopes.

## 4.1 Statement of the problem

In [75], the generation of convex envelopes for general univariate functions was discussed. Here we consider the monomial $x^{2k+1}$ in the range $x \in [a, b]$ where $a < 0 < b$. Let $c, d$ be the $x$-coordinates of the points $C, D$ where the tangents from points $A$ and $B$ respectively meet the curve (see Figure 4.1 below). The shape of the convex underestimator of $x^{2k+1}$ depends on the relative magnitude of $b$ and $c$. In particular, if $c < b$ (as is the case in Figure 4.1), a convex underestimator can be formed from the tangent from $x = a$ to $x = c$ followed by the curve $x^{2k+1}$ from $x = c$ to $x = b$. On the other hand, if $c > b$ (cf. Figure 4.2), a convex underestimator is simply the straight line passing through $A$ and $B$.

The situation is similar for the concave overestimator of $x^{2k+1}$ in the range $x \in [a, b]$. If $d > a$, the overestimator is given by the upper tangent from $B$ to $D$ followed by the curve $x^{2k+1}$ from $D$ to $A$, as shown in Figure 4.1. On the other hand, if $d < a$, the overestimator is just the straight line from $A$ and $B$. It should be noted that the conditions $c > b$ and $d < a$ cannot both hold simultaneously.

<figure><img src="image_10.png" /><figcaption>Figure 4.1: Convex envelope of $x^{2k+1}$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 4.2: The case when $c > b$.</figcaption></figure>

## 4.2 The tangent equations

The discussion in Section 4.1 indicates that forming the envelope of $x^{2k+1}$ requires the determination of the tangents that pass through points A, C and B, D. Considering the first of these two tangents and equating the slope of the line $\overline{AC}$ to the gradient of $x^{2k+1}$ at $x = c$, we derive the tangency condition:

$$ \frac{c^{2k+1} - a^{2k+1}}{c - a} = (2k + 1)c^{2k} \quad (4.1) $$

Hence c is a root of the polynomial:

$$ P^k(x, a) \equiv (2k)x^{2k+1} - a(2k+1)x^2 + a^{2k+1} \quad (4.2) $$
It can be shown by induction on $k$ that:

$$P^k(x, a) = a^{2k-1}(x-a)^2 Q^k\left(\frac{x}{a}\right) \qquad (4.3)$$

where the polynomial $Q^k(x)$ is defined as:

$$Q^k(x) = 1 + \sum_{i=2}^{2k} ix^{i-1}. \qquad (4.4)$$

Thus, the roots of $P^k(x, a)$ can be obtained from the roots<fnref n="1" /> of $Q^k(x)$. Unfortunately, polynomials of degree greater than 4 cannot generally be solved by radicals (what is usually called an “analytic solution”). This is the case for $Q^k(x)$ for $k > 2$. For example, the Galois group of $Q^3(x) = 6x^5 + 5x^4 + 4x^3 + 3x^2 + 2x + 1$ is isomorphic to $S_5$ (i.e. the symmetric group of order 5) which is not soluble since its biggest proper normal subgroup is $A_5$, the smallest non-soluble group. For details on Galois theory and the solvability of polynomials, see [118].

## 4.3 The roots of $Q^k(x)$ and their uniqueness

Unlike $P^k(x, a)$, the polynomial $Q^k(x)$ does not depend on the range of $x$ being considered. Moreover, as shown formally in Section 4.3.1 below, $Q^k(x)$ has exactly one real root, $r_k$, for any $k \ge 1$, and this lies in $[-1 + 1/2k, -0.5]$. Hence, the roots of $Q^k(x)$ for different $k$ can be computed *a priori* to arbitrary precision using simple numerical schemes (e.g. bisection). A table of these roots is presented in Table 4.1 for $k \le 10$.

<table><caption>Table 4.1: Numerical values of the roots of $Q^k(x)$ for $k = 1, \dots, 10$ (to 10 significant digits).</caption><thead><tr><th>k</th><th>r<sub>k</sub></th><th>k</th><th>r<sub>k</sub></th></tr></thead><tbody><tr><td>1</td><td>-0.5000000000</td><td>6</td><td>-0.7721416355</td></tr><tr><td>2</td><td>-0.6058295862</td><td>7</td><td>-0.7921778546</td></tr><tr><td>3</td><td>-0.6703320476</td><td>8</td><td>-0.8086048979</td></tr><tr><td>4</td><td>-0.7145377272</td><td>9</td><td>-0.8223534102</td></tr><tr><td>5</td><td>-0.7470540749</td><td>10</td><td>-0.8340533676</td></tr></tbody></table>

<footnote n="1">¹Although $P^k(x, a)$ has the additional root $x = a$, this is not of practical interest.</footnote>
### 4.3.1 Bounding the roots of $Q^k(x)$

In this section, we show that $Q^k(x)$ has exactly one real root, which lies in the interval $[-1 + \frac{1}{2k}, -\frac{1}{2}]$.

#### 4.3.1 Proposition

For all $k \in \mathbb{N}$, the following properties hold:

$$ \left. \begin{array}{l} Q^k(0) = 1 \\ Q^k(-1) = -k \end{array} \right\} \qquad (4.5) $$

$$ \forall x > 0 \left( \frac{dQ^k(x)}{dx} > 0 \right) \qquad (4.6) $$

$$ \forall x \le -1 \ (Q^k(x) < 0) \qquad (4.7) $$

**Proof.** (4.5): $Q^k(0) = 1$ by direct substitution in (4.4). Also $Q^k(-1) = 1 + \sum_{i=2}^{2k} i(-1)^{i-1} = \sum_{i=1}^{k} (2i - 1) - \sum_{i=1}^{k} 2i = -k$.

(4.6): $\frac{dQ^k(x)}{dx} = \sum_{i=1}^{2k-1} i(i+1)x^{i-1}$, hence it is greater than zero whenever $x > 0$.

(4.7): For $x \ne 0$, we can rewrite $Q^k(x)$ as $\sum_{i=1}^{k} x^{2i-2}[2i(x+1) - 1]$. For $x \le -1$, we have $x^{2i-2} > 0$ and $[2i(x+1) - 1] < 0$, thus each term of the sum is negative. $\square$

From the above proposition and the continuity of $x^{2k+1}$, we can conclude that:

1. there is at least one root between -1 and 0 (property (4.5));

2. there are no roots for $x \ge 0$ (property (4.6) and the fact that $Q^k(0) > 0$);

3. there are no roots for $x \le -1$ (property (4.7)).

#### 4.3.2 Lemma

For all $k \in \mathbb{N}$, the real roots of $Q^k(x)$ lie in the interval $[-1 + \frac{1}{2k}, -\frac{1}{2}]$.

**Proof.** This is proved by induction on $k$. For $k = 1$, $Q^1(x) = 1 + 2x$ has one real root at $x = -\frac{1}{2}$ which lies in the set $[-1 + \frac{1}{2}, -\frac{1}{2}]$. In particular, $Q^1(x) < 0$ for all $x < -1 + \frac{1}{2}$ and $Q^1(x) > 0$ for all $x > -\frac{1}{2}$.

We now make the inductive hypothesis that, for all $j < k$, $Q^j(x) > 0$ for all $x > -\frac{1}{2}$ and $Q^j(x) < 0$ for all $x < -1 + \frac{1}{2j}$ and prove that the same holds for $j = k$. Using (4.4), we can
write $Q^k(x) = Q^{k-1}(x) + x^{2k-2}(2kx + 2k - 1)$ for all $k > 1$. Since $x^{2k-2}$ is always positive, we have that:

$$Q^k(x) > Q^{k-1}(x) \quad \text{if} \quad x > -1 + \frac{1}{2k}$$

$$Q^k(x) < Q^{k-1}(x) \quad \text{if} \quad x < -1 + \frac{1}{2k}$$

for all $k > 1$. Now, since $-\frac{1}{2} > -1 + \frac{1}{2k}$ for all $k > 1$, $Q^k(x) > Q^{k-1}(x) > 0$ for all $x > -\frac{1}{2}$ by inductive hypothesis.

Furthermore, by the inductive hypothesis, $Q^{k-1}(x) < 0$ for all $x < -1 + \frac{1}{2(k-1)}$; since $\frac{1}{2k} < \frac{1}{2(k-1)}$, it is also true that $Q^{k-1}(x) < 0$ for all $x < -1 + \frac{1}{2k}$. But since, as shown above, $Q^k(x) < Q^{k-1}(x)$ for all $x < -1 + \frac{1}{2k}$, we can deduce that $Q^k(x) < 0$ for all $x < -1 + \frac{1}{2k}$.

We have thus proved that, for all $k > 0$,

$$Q^k(x) > 0 \quad \text{if} \quad x > -\frac{1}{2} \qquad (4.8)$$

$$Q^k(x) < 0 \quad \text{if} \quad x < -1 + \frac{1}{2k}. \qquad (4.9)$$

The proof of the lemma follows from (4.8), (4.9) and the continuity of $Q^k(x)$. $\square$

### 4.3.3 Theorem

For all $k \in \mathbb{N}$, $Q^k(x)$ has exactly one real root, which lies in the interval $[-1 + \frac{1}{2k}, -\frac{1}{2}]$.

*Proof.* Consider the polynomial $P^k(x, 1) = 2kx^{2k+1} - (2k+1)x^{2k} + 1$ defined by (4.2). By virtue of (4.3), we have the relation $P^k(x, 1) = (x - 1)^2 Q^k(x)$. Consequently, $P^k(x, 1)$ and $Q^k(x)$ have exactly the same roots for $x < 1$. Therefore (Lemma 4.3.2), all negative real roots of $P^k(x, 1)$ lie in the interval $[-1 + 1/2k, -1/2]$, and there is at least one such root.

Now, $P^k(x, 1)$ can be written as $P^k(x, 1) = q_1^k(x) + q_2^k(x) + 1$, where $q_1^k(x) = 2kx^{2k+1}$ and $q_2^k(x) = -(2k+1)x^{2k}$. Since $q_1$ is a monomial of odd degree, it is monotonically increasing in $[-1, 0]$. Since $q_2$ is a monomial of even degree with a negative coefficient, it is also monotonically increasing in $[-1, 0]$.

Overall, then, $P^k(x, 1)$ is monotonically increasing in $[-1, 0]$, and consequently in the interval $[-1 + 1/2k, -1/2]$ where all its negative real roots lie. Therefore, there can be only one such root, which proves that $Q^k(x)$ also has a unique root in this interval. $\square$
## 4.4 Nonlinear convex envelopes

If the roots shown in the second column of Table 4.1 are denoted by $r_k$, then the tangent points $c$ and $d$ in Figure 4.1 are simply $c = r_k a$ and $d = r_k b$. The lower and upper tangent lines are given respectively by:

$$a^{2k+1} + \frac{c^{2k+1} - a^{2k+1}}{c-a} (x-a) \qquad (4.10)$$

$$b^{2k+1} + \frac{d^{2k+1} - b^{2k+1}}{d-b} (x-b) \qquad (4.11)$$

Hence, the envelope for $z = x^{2k+1}$ when $x \in [a, b]$ and $a < 0 < b$:

$$l_k(x) \le z \le u_k(x) \qquad (4.12)$$

is as follows:

* If $c < b$, then:

$$l_k(x) = \begin{cases} a^{2k+1} (1 + R_k (\frac{x}{a} - 1)) & \text{if } x < c \\ x^{2k+1} & \text{if } x \ge c \end{cases} \qquad (4.13)$$

otherwise:

$$l_k(x) = a^{2k+1} + \frac{b^{2k+1} - a^{2k+1}}{b-a} (x-a) \qquad (4.14)$$

* If $d > a$, then:

$$u_k(x) = \begin{cases} x^{2k+1} & \text{if } x \le d \\ b^{2k+1} (1 + R_k (\frac{x}{b} - 1)) & \text{if } x > d \end{cases} \qquad (4.15)$$

otherwise:

$$u_k(x) = a^{2k+1} + \frac{b^{2k+1} - a^{2k+1}}{b-a} (x-a) \qquad (4.16)$$

where we have used the constant $R_k \equiv \frac{r_k^{2k+1}-1}{r_k-1}$. If the range of $x$ is unbounded either below or above we take the limits of $l_k(x)$ and $u_k(x)$ as $a \to \infty$ or $b \to \infty$.

By construction, the above convex underestimators and concave overestimators are continuous and differentiable everywhere. Moreover, they form the convex envelope of of $x^{2k+1}$, as the following theorem shows.

### 4.4.1 Theorem

The convex underestimator and concave overestimator of $x^{2k+1}$ for $x \in [a, b]$ where $a < 0 < b$ and $k \in \mathbb{N}$, given in equations (4.12)-(4.16), are as tight as possible.
*Proof.* First, consider the case where $a < d < 0 < c < b$. As the convex underestimator between $c$ and $b$ is the curve itself, no tighter one can be found in that range. Furthermore, the convex underestimator between $a$ and $c$ is a straight line connecting two points on the original curve, so again it is the tightest possible.

It only remains to show that $l_k(x)$ is convex for any small neighbourhood of $c$. Consider the open interval $(c - \varepsilon, c + \varepsilon)$, and the straight line segment $\Gamma(c, c + \varepsilon)$ with endpoints $(c, l_k(c))$, $(c + \varepsilon, l_k(c + \varepsilon))$. Because for all $x \ge c$, $l_k(x) \equiv x^{2k+1}$ is convex, all points in $\Gamma(c, c + \varepsilon)$ lie above the underestimator. If we now consider $\Gamma(c - \varepsilon, c + \varepsilon)$, its slope is smaller than the slope of $\Gamma(c, c + \varepsilon)$ (because the point with coordinate $c - \varepsilon$ moves on the tangent of the curve at $c$), yet the right endpoint $c + \varepsilon$ of the segments is common. Thus all points in $\Gamma(c - \varepsilon, c + \varepsilon)$ also lie above the underestimator $l_k(x)$. Since $\varepsilon$ is arbitrary, the claim holds. A similar argument holds for the overestimator between $a$ and $d$.

The cases where $a < d < 0 < b \le c$ and $d \le a < 0 < c < b$ are simpler as the underestimator is a straight line in the whole range of $x$. $\square$

## 4.5 Tight linear relaxation

The convex envelope presented in Section 4.4 is nonlinear. As convex relaxations are used to solve a local optimization problem at each node of the search tree examined by the sBB algorithm, using a linear relaxation instead of a nonlinear one may have a significant impact on computational cost. We can relax the nonlinear envelope to a linear relaxation by dropping the “follow the curve” requirements on either side of the tangency points, and using the lower and upper tangent as convex underestimator and concave overestimator respectively, as follows:

$$a^{2k+1} \left(1 + R_k \left(\frac{x}{a} - 1\right)\right) \le z \le b^{2k+1} \left(1 + R_k \left(\frac{x}{b} - 1\right)\right) \quad (4.17)$$

We can tighten this relaxation further by drawing the tangents to the curve at the endpoints $A, B$, as shown in Figure 4.3. This is equivalent to employing the following constraints:

$$(2k + 1)b^{2k}x - 2kb^{2k+1} \le z \le (2k + 1)a^{2k}x - 2ka^{2k+1} \quad (4.18)$$

in addition to those in (4.17).

As has been noted in Section 4.1, when $c > b$, the underestimators on the left hand sides of (4.17) and (4.18) should be replaced by the line $a^{2k+1} + \frac{b^{2k+1}-a^{2k+1}}{b-a}(x-a)$ through points $A$ and
<figure><img src="image_12.png" /><figcaption>Figure 4.3: Tight linear relaxation of $x^{2k+1}$.</figcaption></figure>

B (see Figure 4.2). On the other hand, if $d < a$, this line should be used to replace the concave overestimators on the right hand sides of (4.17) and (4.18). The linear relaxation constraints are summarized in Table 4.2.

<table><caption>Table 4.2: Summary of linear relaxations for $z = x^{2k+1}$, $x \in [a, b]$, $a < 0 < b$.</caption><thead><tr><th>c &lt; b and d &gt; a</th><th>c &gt; b and d &gt; a</th><th>c &lt; b and d &lt; a</th></tr></thead><tbody><tr><td>z &ge; a<sup>2k+1</sup>(1 + R<sub>k</sub>(x/a - 1))</td><td>z &ge; a<sup>2k+1</sup> + <sup>b<sup>2k+1</sup>-a</sup>&frasl;<sub>b-a</sub>(x - a)</td><td>z &ge; a<sup>2k+1</sup>(1 + R<sub>k</sub>(x/a - 1))</td></tr><tr><td>z &le; b<sup>2k+1</sup>(1 + R<sub>k</sub>(x/b - 1))</td><td>z &le; b<sup>2k+1</sup>(1 + R<sub>k</sub>(x/b - 1))</td><td>z &le; a<sup>2k+1</sup> + <sup>b<sup>2k+1</sup>-a</sup>&frasl;<sub>b-a</sub>(x - a)</td></tr><tr><td>z &ge; (2k + 1)b<sup>2k</sup>x - 2kb<sup>2k+1</sup></td><td>z &le; (2k + 1)a<sup>2k</sup>x - 2ka<sup>2k+1</sup></td><td>-</td></tr><tr><td>z &le; (2k + 1)a<sup>2k</sup>x - 2ka<sup>2k+1</sup></td><td>-</td><td>z &ge; (2k + 1)b<sup>2k</sup>x - 2kb<sup>2k+1</sup></td></tr></tbody></table>

## 4.6 Comparison to other relaxations

This section considers two alternative convex/concave relaxations of the monomial $x^{2k+1}$ where the range of $x$ includes 0, and compares them with both the nonlinear envelope and linear relaxation proposed in this paper.
4.6.1 Reformulation in terms of bilinear products

One possible way of determining a convex/concave relaxation for $z = x^{2k+1}$, where $a \le x \le b$ and $a < 0 < b$, is via its exact reformulation in terms of a bilinear product of $x$ and the convex monomial $x^{2k}$:

$$
\begin{align*}
z &= wx \\
w &= x^{2k} \\
a &\le x \le b \\
0 &\le w \le w^U = \max\{a^{2k}, b^{2k}\}
\end{align*}
$$

By replacing the bilinear term $wx$ with the standard linear convex envelope proposed by [80]
(see equations (5.10)-(5.13)), and the convex univariate term $x^{2k}$ with the envelope given by the
function itself as the underestimator and the secant as the overestimator, we obtain the following
constraints:

$$
\begin{align*}
aw &\le z \le bw \\
w^U x + bw - w^U b &\le z \le w^U x + aw - w^U a \\
x^{2k} &\le w \le a^{2k} + \frac{b^{2k} - a^{2k}}{b-a}(x-a) \\
a &\le x \le b \\
0 &\le w \le w^U
\end{align*}
$$

After some algebraic manipulation, we can eliminate $w$ to obtain the following nonlinear
convex relaxation for $z$:

$$
\frac{w^U a}{a-b} (x-b) \le z \le \frac{w^U b}{b-a} (x-a) \quad (4.19)
$$

$$
bx^{2k} + w^U (x-b) \leq z \leq ax^{2k} + w^U (x-a) \quad (4.20)
$$

$$
a \left( a^{2k} + \frac{b^{2k} - a^{2k}}{b-a} (x-a) \right) \le z \le b \left( a^{2k} + \frac{b^{2k} - a^{2k}}{b-a} (x-a) \right) \quad (4.21)
$$

Figure 4.4 shows the convex/concave relaxation for $x^3$ for $x \in [-1, 1]$ obtained using (4.19)-(4.21). It also compares it with the nonlinear envelope of Section 4.4 (dashed lines in Figure 4.4a) and the linear relaxation of Section 4.5 (dashed lines in Figure 4.4b).