Linear And Quadratic Optimization Models

Paritosh Mokhasi, Kernel Developer

Join the Conversation #WolframTechConf

ABSTRACT

This talk will provide an overview of new capabilities for solving linear and quadratic convex optimization problems. Examples will be shown that demonstrate how a large class of optimization problems can be transformed into linear or quadratic forms. Different applications will be used to demonstrate the use of new convex optimization solver functions and to provide explanation on interpreting the results through the use of dual formulations and sensitivity analysis.

Convex Optimization

A convex optimization problem is one of the form
			
				minimize    (Subscript[f, 0])(x)
				subject to    (Subscript[f, i])(x), 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.

Convex Optimization Models

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

Convex Optimization Models

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

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

Convex Optimization Models

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

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.

Convex Optimization Models

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

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.

Convex Optimization Models

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

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[A, i]).x+(Subscript[b, i])(Subscript[||, 2])≤(Subscript[α, i])x+(Subscript[β, i]).

Some Popular Transformations.

•  Convert ||a.x-b(Subscript[||, 1]) using ||a.x-b(Subscript[||, 1])=(∑_i)(Subscript[s, i]),(Subscript[s, i])=|(Subscript[a, i]).x-(Subscript[b, i])| and -(Subscript[s, i])≤ (Subscript[a, i]).x-(Subscript[b, i])≤ (Subscript[s, i]), i=1,2,..,n.

•  Convert ||a.x-b(Subscript[||, ∞]) using ||a.x-b(Subscript[||, ∞])=Max(Subscript[., i])|(Subscript[a, i]).x-(Subscript[b, i])|=s and -s≤ (Subscript[a, i]).x-(Subscript[b, i])≤ s, i=1,2,..,n.

•  Convert |x| into using |x|=(SuperPlus[x])+(SuperMinus[x]), x=(SuperPlus[x])-(SuperMinus[x]), (SuperPlus[x])≥0,(SuperMinus[x])≥0. Alternatively t=|x|,-t≤x≤t.

•  Convert Max(Subscript[., i]) ((Subscript[a, i]).x-(Subscript[b, i])) into a linear function using s=Max(Subscript[., i]) ((Subscript[a, i]).x-(Subscript[b, i])), (Subscript[a, i])x-(Subscript[b, i])≤s, i=1,2,..,n.

•  Convert Min(Subscript[., i]) ((Subscript[a, i]).x-(Subscript[b, i])) into a linear function using -s=Min(Subscript[., i]) ((Subscript[a, i]).x-(Subscript[b, i])), (Subscript[a, i])x-(Subscript[b, i])≥s, i=1,2,..,n.

Linear Optimization

Minimize    	c.x
										subject to	a.x+b≥0,
(Subscript[a, eq]).x+(Subscript[b, eq])=0

LinearOptimization[objective, constraints, vars]
LinearOptimization[c,{a,b}]
LinearOptimization[c,{a,b},{aeq,beq}]

• Also called LinearProgramming (LP).

• Dedicated solvers capable of handling machine and arbitrary precision computations.

Ways to provide input

LinearOptimization[objective, constraints, vars]
LinearOptimization[c,{a,b}]
LinearOptimization[c,{a,b},{aeq,beq}]

Minimize x+y subject to the constraints x+2y≥3,x≥-1,y≥-1:

LinearOptimization[x+y,{x+2y≥3,x≥-1,y≥-1},{x,y}]

LinearOptimization[{1,1},{{{1,2},{1,0},{0,1}},{-3,1,1}}]

Ways to provide input

LinearOptimization[objective, constraints, vars]

Minimize x+y subject to the constraints x+2y≥3,x≥-1,y≥-1:

Use a vector variable:

LinearOptimization[{1,1}.v,{{1,2}.v≥3,{1,0}.v≥-1,{0,1}.v≥-1},v]

Specify multiple constraints at once using VectorGreaterEqual (or VectorLessEqual):

LinearOptimization[{1,1}.v,VectorGreaterEqual[{{{1,2},{1,0},{0,1}}.v,{3,-1,-1}}],v]

Use Esc v>= Esc for compact VectorGreaterEqual notation:

LinearOptimization[{1,1}.v,{{1,2},{1,0},{0,1}}.v\[VectorGreaterEqual]{3,-1,-1},v]

Use parameter equations:

LinearOptimization[c.v,{a.v\[VectorGreaterEqual]b,c=={1,1},a=={{1,2},{1,0},{0,1}},b=={3,-1,-1}},v]

Linear Optimization: Robust data fitting

Given data {{(Subscript[x, 1]),(Subscript[y, 1])},{(Subscript[x, 2]),(Subscript[y, 2])},...,{(Subscript[x, n]),(Subscript[y, n])}},  Find a robust fit to nonlinear discrete data by minimizing ||a.λ-b(Subscript[||, 1]):

SeedRandom[123];data=Table[{x,Sin[2 x]+Cos[(x^2)]+RandomVariate[NormalDistribution[0,0.4]]},{x,RandomReal[{-3,3},500]}];
dataPlot = ListPlot[data, PlotRange->All]

fit the data using the bases (Subscript[ϕ, i])(x)=(((1+(((x-(Subscript[c, i])))^2)))^(1/2)). The interpolating function will be s(x)=(∑_(i=1))(Subscript[λ, i])(Subscript[ϕ, i])(x). Construct the input data matrix:

bases = Sqrt[1+(x-#)^2]&/@Subdivide[-3,3,10];
input=DesignMatrix[data,bases,x,IncludeConstantBasis->False];

Find the coefficients (Subscript[λ, i]) of the interpolating function. Use transformation ||a.λ-b(Subscript[||, 1])=(∑_i)(Subscript[s, i]),(Subscript[s, i])=|(Subscript[a, i]).λ-(Subscript[b, i])| and -(Subscript[s, i])≤ (Subscript[a, i]).λ-(Subscript[b, i])≤ (Subscript[s, i]), i=1,2,..,n.:

res=LinearOptimization[Total[s],{a==input,b==data[[All,2]],-s\[VectorLessEqual]a.λ-b\[VectorLessEqual]s},{λ,s}];
Short[res,6]

Visualize the fit:

Show[dataPlot,Plot[Evaluate[bases.λ/.res],{x,-3,3}]]

Linear Optimization: Classification Problem

SeedRandom[1234];
set1=RandomReal[{0,1/2},{40,2}];
set2=RandomReal[{1/2,1},{60,2}];
dataPlot=ListPlot[{set1,set2},PlotMarkers->{Automatic, PointSize[0.025]}]

For separation, set 1 must satisfy a.(Subscript[p, i])+b \[VectorGreaterEqual]1and set 2 must satisfy a.(Subscript[q, j])+b\[VectorLessEqual]-1:

res=LinearOptimization[0,{set1.a+b\[VectorGreaterEqual]1,set2.a+b\[VectorLessEqual]-1},{a,b}]

Show[RegionPlot[Evaluate[a.{x,y}+b ≤0/.res],{x,0,1},{y,0,1}],dataPlot]

Linear Fractional Optimization

Minimize 	(α.x+β)/(γ.x+δ) 
									subject to	a.x+b≥0,
(Subscript[a, eq]).x+(Subscript[b, eq])=0

LinearFractionalOptimization[objective, constraints, vars]
LinearFractionalOptimization[{α,β,γ,δ},{a,b}]
LinearFractionalOptimization[{α,β,γ,δ},{a,b},{aeq,beq}]

• If γ=0 and δ=1, problem reverts to linear optimization.

• Linear fractional optimization is solved by transforming the problem into a linear optimization: z = (x/(γ x + δ)), t=(1/(γ x + δ))=>x=z/t

Minimize 	α.z + β t 
					subject to 	(Subscript[a, i]).z - (Subscript[b, i])t ≥ 0, 
γ.z+δ t=1.

• Additional constraint: If γ.v+δ>0 then t≥0 else t≤0.

• Typically used when you want to maximize and minimize an objective simultaneously.

• There could be cases in which the objective may be singular:

Plot3D[(3+(Subscript[x, 1])+2(Subscript[x, 2]))/(1-5(Subscript[x, 1])-5(Subscript[x, 2])),{(Subscript[x, 1]),(Subscript[x, 2])}∈ImplicitRegion[3 (Subscript[x, 1])+4 (Subscript[x, 2])≤7&&5 (Subscript[x, 1])-6 (Subscript[x, 2])≤20&&-7 (Subscript[x, 1])+8 (Subscript[x, 2])≥27&&(Subscript[x, 1])≥-3,{(Subscript[x, 1]),(Subscript[x, 2])}]]

Linear Fractional Optimization: Transportation Problem

Choose a number of units to transport for five goods by minimizing cost and maximizing profit. The cost and profit of transporting one unit of each good are:

costPerUnit = {100,100,100,100,100};
profitPerUnit={115,108,130,120,110};

The weight of each unit of goods and the maximum units available for transport are:

weightPerUnit={100,100,110,200,115};
maxUnits={30,20,20,30,20};

The minimum operating cost of running the truck is $2000 for the trip. The total transportation cost is:

cost=2000+costPerUnit.x;

The truck makes a minimum profit of $1000 if it rides empty. The total profit is:

profit=1000+profitPerUnit.x;

The truck can carry a maximum of 8000 lbs. The constraints on the truck are:

wtConstraint=weightPerUnit.x\[VectorLessEqual]8000;
goodsConstraints=0\[VectorLessEqual]x\[VectorLessEqual]maxUnits;

The optimal mix of goods to transport if found by minimizing the ratio of cost to profit:

res = LinearFractionalOptimization[cost/profit,{wtConstraint,goodsConstraints},x∈Vectors[5]]

TableForm[{x/.res}, TableHeadings->{{"Units"},Array[(Subscript[x, #])&,5]}]

Quadratic Optimization

Minimize 	(1/2)xᵀ.Q.x + c.x 
									subject to	a.x+b≥0,
(Subscript[a, eq]).x+(Subscript[b, eq])=0

QuadraticOptimization[objective, constraints, vars]
QuadraticOptimization[{Q,c},{a,b}]
QuadraticOptimization[{Q,c},{a,b},{aeq,beq}]]
QuadraticOptimization[{{Qf,0},c},...]
QuadraticOptimization[{{Qf,cf},c},...]

• The matrix Q must be a symmetric positive definite matrix.

• The quadratic part can be specified in the factored form xᵀ.(Subscript[Q, f])ᵀ(Subscript[Q, f]).x  or  (xᵀ.(Subscript[Q, f])ᵀ+(Subscript[c, f])).((Subscript[Q, f]).x+(Subscript[c, f])).

• If Q=0, problem reverts to linear optimization problem.

• Transform the problem into a conic optimization problem and use conic solver:
						
						Minimize 	t
						subject to	(1/2)xᵀ.Q.x+c.x≤ t,
a.x+b≥0.

Quadratic Optimization: Data fitting

Given data {{(Subscript[x, 1]),(Subscript[y, 1])},{(Subscript[x, 2]),(Subscript[y, 2])},...,{(Subscript[x, n]),(Subscript[y, n])}},  Find a robust fit to nonlinear discrete data by minimizing ||a.λ-b||:

SeedRandom[123];data=Table[{x,Sin[2 x]+Cos[(x^2)]+RandomVariate[NormalDistribution[0,0.4]]},{x,RandomReal[{-3,3},500]}];
dataPlot = ListPlot[data, PlotRange->All]

fit the data using the bases (Subscript[ϕ, i])(x)=(((1+(((x-(Subscript[c, i])))^2)))^(1/2)). The interpolating function will be s(x)=(∑_(i=1))(Subscript[λ, i])(Subscript[ϕ, i])(x). Construct the input data matrix:

bases = Sqrt[1+(x-#)^2]&/@Subdivide[-3,3,10];
input=DesignMatrix[data,bases,x,IncludeConstantBasis->False];

Find the coefficients (Subscript[λ, i]) of the interpolating function:

λ=QuadraticOptimization[{{input,-data[[All,2]]},0},{}]

Visualize the fit:

Show[dataPlot,Plot[{Evaluate[bases.λ],Sin[2 x]+Cos[(x^2)]},{x,-3,3}, PlotLegends->{"Fit","Reference"}]]

Quadratic Optimization: Geometry

Find the distance between two convex polyhedra:

SeedRandom[1234];
region1=ConvexHullMesh[RandomReal[1,{50,3}]];
region2=ConvexHullMesh[RandomReal[{1,2},{50,3}]];

{dist,minPts}=QuadraticOptimization[Norm[x1-x2]^2,{x1∈region1,x2∈region2},{x1,x2},{"PrimalMinimumValue", "PrimalMinimizerRules"}]

Show the nearest points with the line connecting them:

Show[region1,region2,Graphics3D[{Line[{x1,x2}],Blue,PointSize[0.025],Point[{x1,x2}]} /. minPts]]

Quadratic Optimization: Using the dual

Find the radius r and center (Subscript[c, 0]) of a minimal enclosing ball that encompasses a given region, i.e., Minimize r subject to ||(Subscript[p, i])-(Subscript[c, 0])||≤r, i=1,2,.n:

region=;
pts=MeshCoordinates[region];

The dual of this problem is to maximize -(∑_i)((((Subscript[w, i])(Subscript[p, i])))^2)+(∑_i)(Subscript[w, i])(∑_(j=1))(((Subscript[p, i])((Subscript[x, j])))^2) subject to (∑_(i-1))(Subscript[w, i])=1,(Subscript[w, i])≥0:

obj=-((Norm[w.pts])^2)+Total[(pts^2),{2}].w;
constraints={w\[VectorGreaterEqual]0,Total[w]==1};

{weights,minval}=QuadraticOptimization[-obj,constraints,w∈Vectors[Length[pts]], {"PrimalMinimizer","PrimalMinimumValue"}];
Short[weights,3]

The center of the minimal enclosing ball is (Subscript[c, 0])=(∑_(i=1))(Subscript[w, i])(Subscript[p, i]). The radius of the minimal enclosing ball is Sqrt of the maximum value:

center=weights.pts

radius=Sqrt[-minval]

Visualize the enclosing ball:

Show[Graphics3D[{Opacity[0.2],Red,Ball[center,radius]}], region]

Second Order Cone Optimization

Minimize 	c.x 
	 								subject to	||(Subscript[A, i]).x+(Subscript[b, i])(Subscript[||, 2])≤(Subscript[α, i])x+(Subscript[β, i]),i=1,2,...,m,

SecondOrderConeOptimization[objective, constraints, vars]
SecondOrderConeOptimization[c,{{(Subscript[A, 1]),(Subscript[b, 1]),(Subscript[α, 1]),(Subscript[β, 1])},…,{(Subscript[A, k]),(Subscript[b, k]),(Subscript[α, k]),(Subscript[β, k])}}]

• If (Subscript[A, i])=0 and (Subscript[b, i])=0, the problem reverts to a linear optimization problem.

• If (Subscript[α, i])=0 and (Subscript[β, i])=0, the constraints represent equality constraints.

• 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

Second Order Cone Optimization: Geometry

Find the minimum distance between two disks of radius 1 centered at (0,0) and (2,1). Using auxiliary variable t, minimize t subject to ||(Subscript[p, 1])-(Subscript[p, 2])||≤t, ||(Subscript[p, 1])||≤1,||(Subscript[p, 2])-(2,1)||≤1:

res=SecondOrderConeOptimization[t,{Norm[p1-p2]≤t,Norm[p1]≤1,Norm[Inactive[Plus][p2,-{2,1}]]≤1},{t,p1∈Vectors[2],p2∈Vectors[2]}]

Visualize the position of the two points:

Graphics[{{Lighter[Blue],Disk[]},{Lighter[Red],Disk[{2,1}]},  {Green,PointSize[0.03],Point[{p1,p2}/.res]}}]

The auxiliary variable t gives the distance between the points:

{t,EuclideanDistance[p1,p2]}/.res

Second Order Cone Optimization: Using the dual

Find the plane that separates two non-intersecting convex polygons:

SeedRandom[12345];
(Subscript[Ω, 1])=ConvexHullMesh[RandomReal[{-1,1},{20,2}]];
(Subscript[Ω, 2])=ConvexHullMesh[({2,1}+#1&)/@RandomReal[{-1,1},{20,2}]];

Let (Subscript[p, 1]) be a point on (Subscript[Ω, 1]). Let (Subscript[p, 2]) be a point on (Subscript[Ω, 2]). Minimize ||(Subscript[p, 1])-(Subscript[p, 2])||. Using auxiliary variable t, the transformed objective is to minimize t subject to ||(Subscript[p, 1])-(Subscript[p, 2])||≤t, (Subscript[p, 1])∈(Subscript[Ω, 1]),(Subscript[p, 2])∈(Subscript[Ω, 2]):

res=SecondOrderConeOptimization[t,{p1∈(Subscript[Ω, 1]),p2∈(Subscript[Ω, 2]),Norm[p1-p2]≤t},{p1,p2,t}]

According to the separating hyperplane theorem, the dual associated with the constraint ||(Subscript[p, 1])-(Subscript[p, 2])||≤t will give the normal of the hyperplane:

{dual,constraints}=SecondOrderConeOptimization[t,{p1∈(Subscript[Ω, 1]),p2∈(Subscript[Ω, 2]),Norm[p1-p2]≤t},{p1,p2,t},{"DualMaximizer","SecondOrderConeConstraints"}]

The constraint i internally has the form {(Subscript[a, i]),(Subscript[b, i]),(Subscript[α, i]),(Subscript[β, i])}. ||(Subscript[p, 1])-(Subscript[p, 2])||≤t is the only constraint containing a matrix a:

i=Position[constraints[[All,1]],x_?MatrixQ][[1,1]]

The hyperplane is constructed as:

plane=Hyperplane[dual[[i,1]],(p1+p2)/2/.res];
Show[(Subscript[Ω, 1]),(Subscript[Ω, 2]), Graphics[{{Point[{p1,p2}/.res]},{Red,plane}}]]

Sensitivity Analysis

Find the mix of products to manufacture that maximizes the profits of a company. The company makes three products. Each product is made using a combination of four machines. The time a machine spends on each product is. Each machine can run a maximum of 2400 minutes per week:

machineTime ={{20,10,10},{12,28,16},{15,6,5},{10,15,0}};
Row[BarChart[machineTime[[#]],,ImageSize->200]&/@Range[4]]

Company can make a maximum of 100 units of (Subscript[x, 1]), 40 units of (Subscript[x, 2]) and 60 units of (Subscript[x, 3]). Profit made by selling (Subscript[x, 1])is $45, (Subscript[x, 2]) is $60 and (Subscript[x, 3]) is $50:

timeConstraint=machineTime.x\[VectorLessEqual]2400;
capacity=x\[VectorLessEqual]{100,40,60};
profit = {45.,60.,50.}.x;

The profit based on the current capacity is:

max0 = -LinearOptimization[-profit,{timeConstraint,capacity},x,"PrimalMinimumValue"]

Sensitivity Analysis

Extract all the data needed to analyze the impact of the capacity change on profit:

{sens,{a,b}}=LinearOptimization[-profit,{timeConstraint,capacity} ,x,
{"ConstraintSensitivity","LinearInequalityConstraints"}]

Extract the sensitive constraints from a.x+b⪰0:

spos=Flatten@Position[Chop[sens[[1]]],_?(Abs[#]>0&)];
Table[a[[i]].{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}+b[[i]]≥0,{i,spos}]//TableForm

If the production capacity for product 3 is increased, there should be an increase in profit:

max1=-LinearOptimization[-profit,{timeConstraint,x\[VectorLessEqual]{100,400,60+10}} ,x,"PrimalMinimumValue"]

This corresponds exactly to the predicted increase from the sensitivity:

sens[[1,spos[[1]]]]×10

Thank you!