Mixed Integer Convex Optimization

Paritosh Mokhasi, Kernel Developer

Join the Conversation #WolframTechConf

ABSTRACT

This talk will provide an overview of solving mixed-integer convex optimization problems using dedicated Mathematica solver functions. Examples will be shown that demonstrate how a variety of problems can be specified using the mixed-integer formulations and how the results can be interpreted. An overview of the underlying algorithm that solves mixed-integer convex problems will also be  given in this talk.

Convex Optimization

A convex optimization problem is one of the form
			
				minimize    (Subscript[f, 0])(x)
				subject to    (Subscript[f, i])(x)≤0, i=1,2,...,m,
				
	where the functions (Subscript[f, 0])(x),...,(Subscript[f, m])(x): (ℝ^n)->ℝ are convex, i.e., satisfy
	
				(Subscript[f, i])(α x+β y)≤ α (Subscript[f, i])(x)+β (Subscript[f, i])(y)
				
	for all x,y∈(ℝ^n) and all α,β∈ℝ with α+β=1,α≥0,β≥0.

Optimization Models

Minimize 	(Subscript[f, 0])(x),
		     subject to	         (Subscript[f, i])(x)≤0, i=1,2,...,m,
         x∈(ℝ^n)

Optimization Models

Minimize     	(Subscript[f, 0])(x),
		 subject to	 (Subscript[f, i])(x), i=1,2,...,m,
 x∈(ℝ^n)

Linear Optimization 			(Subscript[f, 0])(x)=c.x,                                   (Subscript[f, i])(x)=((Subscript[a, i]).x+(Subscript[b, i]))≥0.

Linear Fractional Optimization		(Subscript[f, 0])(x)=(α.x+β)/(γ.x+δ),   (Subscript[f, i])(x)=((Subscript[a, i]).x+(Subscript[b, i]))≥0.

Quadratic Optimization 		(Subscript[f, 0])(x)=(1/2)x.Q.x+c.x,              (Subscript[f, i])(x)=((Subscript[a, i]).x+(Subscript[b, i]))≥0.

Second Order Cone Optimization 	(Subscript[f, 0])(x)=c.x,                                   (Subscript[f, i])(x)=(Subscript[α, i])x+(Subscript[β, i])-||(Subscript[A, i]).x+(Subscript[b, i])(Subscript[||, 2])≥0.

Semidefinite Optimization 		(Subscript[f, 0])(x)=c.x,                                   f(x)=(Subscript[A, 0])+(Subscript[A, 1])(Subscript[x, 1])+…+(Subscript[A, n])(Subscript[x, n]) (Subscript[⪰, (SemidefiniteConeList[n])])0.

Mixed Integer Convex Optimization Models

Minimize 	(Subscript[f, 0])(x),
		     subject to	        (Subscript[f, i])(x), i=1,2,...,m,
(Subscript[x, j])∈ℝ,j=1,2,…,k,
(Subscript[x, s]) ∈ℤ, s=k+1,…,n.

Linear Optimization 			(Subscript[f, 0])(x)=c.x,                                   (Subscript[f, i])(x)=((Subscript[a, i]).x+(Subscript[b, i]))≥0.

Linear Fractional Optimization		(Subscript[f, 0])(x)=(α.x+β)/(γ.x+δ),   (Subscript[f, i])(x)=((Subscript[a, i]).x+(Subscript[b, i]))≥0.

Quadratic Optimization 		(Subscript[f, 0])(x)=(1/2)x.Q.x+c.x,              (Subscript[f, i])(x)=((Subscript[a, i]).x+(Subscript[b, i]))≥0.

Second Order Cone Optimization 	(Subscript[f, 0])(x)=c.x,                                   (Subscript[f, i])(x)=(Subscript[α, i])x+(Subscript[β, i])-||(Subscript[A, i]).x+(Subscript[b, i])(Subscript[||, 2])≥0.

Semidefinite Optimization 		(Subscript[f, 0])(x)=c.x,                                   f(x)=(Subscript[A, 0])+(Subscript[A, 1])(Subscript[x, 1])+…+(Subscript[A, n])(Subscript[x, n]) (Subscript[⪰, (SemidefiniteConeList[n])])0.

Mixed Integer Linear Optimization

Minimize    	c.x
								subject to	a.x+b\[VectorGreaterEqual]0,
(Subscript[a, eq]).x+(Subscript[b, eq])=0,
(Subscript[x, i])∈ℝ, i=1,2,…,m, (Subscript[x, s])∈ℤ, s=m+1,…,n

LinearOptimization[objective, constraints, vars]
LinearOptimization[c,{a,b}, {(Subscript[dom, 1]),(Subscript[dom, 2]),...}]
LinearOptimization[c,{a,b},{aeq,beq},{(Subscript[dom, 1]),(Subscript[dom, 2]),...}]

• Also called Mixed Integer LinearProgramming (MILP).

• Branch and Cut method used to solve mixed integer problems.

• Mixed integer problems can only handle machine precision computations.

Mixed Integer Linear Optimization

Minimize    	c.x
								subject to	a.x+b\[VectorGreaterEqual]0,
(Subscript[a, eq]).x+(Subscript[b, eq])=0,
(Subscript[x, i])∈ℝ, i=1,2,…,m, (Subscript[x, s])∈ℤ, s=m+1,…,n

LinearOptimization[objective, constraints, vars]
LinearOptimization[c,{a,b}, {(Subscript[dom, 1]),(Subscript[dom, 2]),...}]
LinearOptimization[c,{a,b},{aeq,beq},{(Subscript[dom, 1]),(Subscript[dom, 2]),...}]

Minimize x+y subject to the constraints x+2y==3,x≥0.5,y≥0, x∈ℤ, y∈ℝ:

LinearOptimization[x+y,{x+2y==3,x≥0.5,y≥0},{x∈Integers,y}]

Use the equivalent matrix representation:

LinearOptimization[{1,1},{{{1,0},{0,1}},{-0.5,0}},{{{1,2}},{-3}},{Integers,Reals}]

Conic Optimization Models

Mixed Integer Conic Optimization

Minimize                   c.x,
		     subject to	        (Subscript[a, 0]).x+(Subscript[b, 0])=0,
(Subscript[a, j]).x+(Subscript[b, j])(Subscript[\[VectorGreaterEqual], (Subscript[κ, j])])0, j=1... k,
(Subscript[x, j])∈ℝ,j=1,2,…,k,
(Subscript[x, s]) ∈ℤ, s=k+1,…,n.

Mixed Integer Conic Optimization

Minimize 	         c.x,
		     subject to	        (Subscript[a, 0]).x+(Subscript[b, 0])=0,
(Subscript[a, j]).x+(Subscript[b, j])(Subscript[\[VectorGreaterEqual], (Subscript[κ, j])])0, j=1... k,
(Subscript[x, j])∈ℝ,j=1,2,…,k,
(Subscript[x, s]) ∈ℤ, s=k+1,…,n.

Mixed Integer Conic Optimization

Minimize 	         c.x,
		     subject to	        (Subscript[a, 0]).x+(Subscript[b, 0])=0,
(Subscript[a, j]).x+(Subscript[b, j])(Subscript[\[VectorGreaterEqual], (Subscript[κ, j])])0, j=1... k,
(Subscript[x, j])∈ℝ,j=1,2,…,k,
(Subscript[x, s]) ∈ℤ, s=k+1,…,n.

• A quadratic optimization problem can transformed into a second order cone optimization:
								
		Minimize             xᵀ.Q.x+c.x
subject to           (Subscript[a, i]).x+(Subscript[b, i])≥0		Minimize             y+c.x
subject to             xᵀ(Subscript[Q, f])ᵀ.(Subscript[Q, f]).x≤y,
		(Subscript[a, i]).x+(Subscript[b, i])≥0	Minimize             y+c.x
subject to             ||(({{(Subscript[Q, f]).x}, {(((1-y))/2)}}))||≤(((1+y))/2),
		(Subscript[a, i]).x+(Subscript[b, i])≥0

Outer Approximation Algorithm

The original conic problem is split into two sub-problems:

(1) A linear problem with integer constraints:

Min. c.x + (OverHat[c]).z,
s.t: A.x+(OverHat[A]).z=b,
 z∈(ℤ^m),
β.x≥ 0 ∀ β∈T.

* Lemma: Let 𝒦 be a closed convex cone, then x∈𝒦 iff β.x≥0, ∀β∈(SuperStar[𝒦])

(2) A conic problem with the integer values (SuperStar[z])obtained from the linear problem substituted in.

Min(Subscript[., x]) c.x + (OverHat[c]).(SuperStar[z]),
s.t: A.x=b-(OverHat[A]).(SuperStar[z]),
x∈𝒦.

(2a) After solving this problem, β∈(SuperStar[𝒦]) can be obtained from the dual of the primal problem. Then 		       T←T⋃{(Subscript[β, x])}.

Outer Approximation Algorithm: An example

Lets solve the following mixed-integer problem:

ConicOptimization[-x-y,{x≤-0.1,y≥0.2,Norm[{x,y}]≤1.5},{x∈Integers,y}]

1. Solve the relaxed conic problem:

{primal,dual,affList}=ConicOptimization[-x-y,{x≤-0.1,y≥0.2,Norm[{x,y}]≤1.5},{x,y},{"PrimalMinimizerRules","DualMaximizer", "ConicConstraintAffineLists"}, Method->"CSDP"];

2. Find the hyper-plane formed by the dual cone:

{a,b} = dual[[-1]].#&/@affList[[-1]];
hplane1=a.{x,y}+b≥0;
plot=RegionPlot[{Norm[{x,y}]≤1.5 &&x≤-0.1&&y≥0.2, hplane},{x,-2,0},{y,0,2},Epilog->{PointSize[0.02],Point[{x,y}/.primal]}]

3, Solve the linear sub-problem and get the integer values:

xIntVal = x/.LinearOptimization[-x-y,{x≤-0.1,y≥0.2,hplane1},{x∈Integers,y}]

4. Substitute the integer values back into the relaxed conic problem and solve:

{primal,dual,affList}=ConicOptimization[-x-y,{x≤-0.1,y≥0.2,Norm[{x,y}]≤1.5,x==xIntVal},{x,y},{"PrimalMinimizerRules","DualMaximizer", "ConicConstraintAffineLists"},Method->"CSDP"];

5. Find the hyper-plane formed by the dual cone:

{a,b} = dual[[-1]].#&/@affList[[-1]];
hplane2=a.{x,y}+b≥0;
Show[plot,RegionPlot[hplane2,{x,-2,0},{y,0,2}],Graphics[{PointSize[0.02],Point[{x,y}/.primal]}]]

6. Repeat step 3-5 adding more hyper-planes:

x/.LinearOptimization[-x-y,{x≤-0.1,y≥0.2,hplane1,hplane2},{x∈Integers,y}]

Tips and Tricks to modelling with integer variables.

•  Semi-continuous variable x∈0 ⋃{a,b}:
						z a≤x≤z b
z∈{0, 1}

•  Given the linear constraints (Subscript[a, i])x≤b, i=1,2,…,n, at-least k of those constraints should be satisfied:
						(Subscript[a, i]).x≤b+M(1-(Subscript[z, i])),
(Subscript[z, 1])+(Subscript[z, 2])+…+(Subscript[z, n])≥k,
(Subscript[z, i])∈{0, 1}

•  Optimization problem should either satisfy (Subscript[a, 1]).x≤(Subscript[b, 1]) OR (Subscript[a, 2]).x≤(Subscript[b, 2]):
						(Subscript[a, 1]).x≤(Subscript[b, 1])+M z,
(Subscript[a, 2]).x≤(Subscript[b, 2])+M(1-z),
z∈{0, 1}

•  If (Subscript[a, 1]).x≤(Subscript[b, 1]) is satisfied then so must (Subscript[a, 2]).x≤(Subscript[b, 2]):
						(Subscript[a, 1]).x≥ (Subscript[b, 1])-L z,
(Subscript[a, 2]).x≤(Subscript[b, 2])+M(1-z),
z∈{0, 1}

•  Allowing x to only take values from a specified set {(Subscript[a, 1]),(Subscript[a, 2]),…,(Subscript[a, n])} can be modeled as 
						x=(Subscript[z, 1])(Subscript[a, 1])+(Subscript[z, 2])(Subscript[a, 2])+…+(Subscript[z, n])(Subscript[a, n]),
(Subscript[z, 1])+(Subscript[z, 2])+…+(Subscript[z, n])=1,
(Subscript[z, i])∈{0, 1}

•  t=Max.{(Subscript[x, 1]),(Subscript[x, 2]),…,(Subscript[x, n])} can be modeled as 
						(Subscript[x, i])≤t≤(Subscript[x, i])+M(1-(Subscript[z, i])), i=1,2,…,n,
(Subscript[z, 1])+(Subscript[z, 2])+…+(Subscript[z, n])=1,
(Subscript[z, i])∈{0, 1}

•  Exact 1-norm (∑_(i=1)^(i=n))|(Subscript[x, i])|=c:
						x = (SuperPlus[x])-(SuperMinus[x]),
(∑_(i=1)^(i=n))(Subscript[x, i]^(+))+(∑_(i=1)^(i=n))(Subscript[x, i]^(-))=c,
(SuperPlus[x])≤c z,
(SuperMinus[x])≤c(1-z),
(SuperPlus[x])≥0,(SuperMinus[x])≥0,
z∈(({0, 1})^n),(SuperPlus[x]),(SuperMinus[x])∈(ℝ^n)

Solve Sudoku game using binary variables

sudokuDisplay[board_?MatrixQ]:=;
sudokuDisplay[known = SparseArray[,{9,9},_]]

Let (Subscript[z, i,j])∈(ℤ^9) be the variable for sub-square (i,j). Let (Subscript[z, i,j])(k), k=1,2,…,9, be the element k of vector (Subscript[z, i,j]). 
	When (Subscript[z, i,j])(k)=1, then (i,j) holds the number k. Each (i,j) contains only one number, so (Subscript[z, i,j]) can contain only one non-zero element, i.e. 		(∑_(k=1)^(k=9))(Subscript[z, i,j])(k)=1:

squareConstraints=Table[{Total[z[i,j]]==1,0\[VectorLessEqual]z[i,j]\[VectorLessEqual]1,z[i,j]∈Vectors[9,Integers]},{i,9},{j,9}];

Each row must contain all the numbers, i.e. (∑_(j=1)^(j=9))(Subscript[z, i,j])= (Overscript[1, ->]),i=1,2,…,9, where (Overscript[1, ->]) is a nine-dimensional vector of ones:

onesVector=ConstantArray[1,9];
rowConstraints=Table[Sum[z[i,j],{j,9}]==onesVector,{i,9}];

Each column must contain all the numbers, i.e. (∑_(i=1)^(i=9))(Subscript[z, i,j])= (Overscript[1, ->]),j=1,2,…,9:

columnConstraints = Table[Sum[z[i,j],{i,9}]==onesVector,{j,9}];

Each 3x3 block must contain all the numbers, i.e. (∑_(m=1)^(m=3))(∑_(n=1)^(n=3))(Subscript[z, i+m,j+n])= (Overscript[1, ->]),i,j∈{0,3,6}:

blockConstraints=Table[Sum[z[i+m,j+n],{m,3},{n,3}]==onesVector,{i,{0,3,6}},{j,{0,3,6}}];

Collect all the variables:

vars=Flatten[Table[z[i,j],{i,9},{j,9}]];

Convert the known values into constraints. If sub-square (i,j) holds number k, then (Subscript[z, i,j])(k)=1:

knownConstraints = MapThread[Indexed[z@@#1,#2]==1&,{known["NonzeroPositions"],known["NonzeroValues"]}];

Solve the puzzle:

res=LinearOptimization[0,{squareConstraints, rowConstraints, columnConstraints, blockConstraints,knownConstraints},vars];

Visualize the result:

sudokuDisplay[SparseArray@MapThread[List@@#1->Range[9].#2&,{vars,vars/.res}]]

Path planning around obstacles:

Find the shortest path between two points while avoiding obstacles. Specify the obstacle:

SeedRandom[1234];mesh = ConvexHullMesh[RandomReal[{.2,.8},{10,2}]];

Extract the half-spaces that form the convex obstacle:

{a,b} =  LinearOptimization[0,{},Element[x,mesh],"LinearInequalityConstraints"];

Specify the start and end points of the path:

start={0,0};
end={1,1};
Show[mesh, Graphics[{PointSize[0.02],Point[{start, end}]}]]

The path can be discretized using n points. Let (Subscript[p, i])∈(ℝ^2) represents the position vector:

n=10;
pvars=Table[Element[p[i],Vectors[2,Reals]],{i,1,n}];

The objective is to minimize (∑_(i=1)^(n-1))||(Subscript[p, i+1])-(Subscript[p, i])||. Let (Subscript[t, i])=(Subscript[p, i+1])-(Subscript[p, i]). The objective is transformed to  (∑_(i=1)^(n-1))((Subscript[t, i])^2):

tvars=Table[t[i],{i,n-1}];
transformationConstraint=Table[t[i]==p[i+1]-p[i],{i,n-1}];
objective=Sum[t[i].t[i],{i,n-1}];

Specify the end point constraints:

endPointConstraints={p[1]==start,p[n]==end};

The distance between any two subsequent points should not be too large:

maxDist =2 Norm[start - end]/n;
separationConstraints =Table[-maxDist \[VectorLessEqual]t[i]\[VectorLessEqual]maxDist ,{i,n-1}];

A point (Subscript[p, i]) is outside the object if at-least one element of (Subscript[a, k]).(Subscript[p, i])+(Subscript[b, k]),k=1,2,…,m is less than zero. 
	To enforce this constraint, let (Subscript[z, k])∈(ℤ^n) be a decision vector and (Subscript[z, k])(i) is the (i^th) element of z such that (∑_(i=1)^n)(Subscript[z, k])(i)≤m-1, then 				a.(Subscript[p, i])+b\[VectorLessEqual]M (Subscript[z, i]),i=1,2,…,n and M is large enough such that M≥Max.|a.(Subscript[p, i])+b|:

m = 10;
zvars = Table[Element[z[i], Vectors[Length[a],Integers]],{i,n}];
binaryConstraints = Table[0 \[VectorLessEqual]z[i] \[VectorLessEqual]1,{i,n}];
avoidanceConstraints = Table[{Inactive[Plus][a.p[i], b]\[VectorLessEqual] m z[i],Total[z[i]]≤ Length[a] - 1},{i,n}];

Find the minimum distance path around the obstacle:

res = QuadraticOptimization[objective, {endPointConstraints, transformationConstraint,separationConstraints, binaryConstraints, avoidanceConstraints},Join[pvars,tvars, zvars]];

Extract and display the path:

path = Table[p[i],{i,n}] /. res;
Show[mesh, Graphics[{Point[path],Line[path]}]]

Graph Coloring Problem

Find the minimum number of colors needed to color the nodes of a graph such that no two adjacent nodes have the same color. 
	Start with more than a sufficient number of colors for any graph:

nColors=6;

Generate a random graph:

SeedRandom[1234];
nNodes=30;
graph=RandomGraph[DegreeGraphDistribution[Table[3,nNodes]],VertexSize->{"Scaled",0.02}]

Extract the adjacency matrix associated with the graph. Since the matrix is symmetric only use the upper-triangular part of the matrix:

graphMat=UpperTriangularize[AdjacencyMatrix[graph]];

Let x be a matrix such that (Subscript[x, i,j])=1 if node i is assigned color j. Let y be a vector such that (Subscript[y, k])=1 if color k is used in the final set. 
	Each node must be given a color:

nodeConstraint=Total[x,{2}]==1;

If two nodes share an edge, then they both cannot be assigned the same color:

sameColorConstraint=Table[{i,j}=ci;Table[Indexed[x,{i,k}]+Indexed[x,{j,k}]≤1,{k,nColors}],{ci,graphMat["NonzeroPositions"]}];

If node i is assigned color k, then color k is used in the final set of colors:

colorConstraint=Table[Indexed[x,{i,k}]≤Indexed[y,k],{i,nNodes},{k,nColors}];

The matrix x and the vector y are binary variables:

binaryConstraint={0\[VectorLessEqual]x\[VectorLessEqual]1,0\[VectorLessEqual]y\[VectorLessEqual]1};

Solve the graph coloring problem:

res=LinearOptimization[Total[y],{nodeConstraint,sameColorConstraint,
colorConstraint,binaryConstraint},{x∈Matrices[{nNodes,nColors},Integers],y∈Vectors[nColors,Integers]}]

The minimum number of colors used are:

Total[y]/.res

Color the graph:

colors={Red,Green,Blue,Brown,Yellow,Black};
HighlightGraph[graph,SparseArray[x/.res]["NonzeroPositions"]/.{i_,j_}:>Style[i,colors[[j]]]]

Thank you!