Wolfram Optimization

Paritosh Mokhasi, 
Kernel Developer, Algorithms R&D,
Wolfram Research Inc.

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

ABSTRACT

This talk will give an overview of the various optimization functions that can be used to solve a wide variety of convex, non-convex and multi-domain problems. The Wolfram optimization functionality will be
demonstrated using a diverse set of examples. Visiting this talk will
enable you to solve optimization type problems efficiently.

Core areas of optimization

Optimization Modeling

Linear and Convex Optimization

Local Optimization

Global Optimization

Mixed Integer Optimization

Parametric Optimization

Symbolic Optimization

Optimal Fitting

Special Purpose Optimization

Linking to commercial solvers.

Optimization functions

WordCloud[{NMinimize, Minimize,FindMinimum, FindFit,FindRoot,FindArgMin,FindMinValue,ConvexOptimization, ParametricConvexOptimization,ConicOptimization,SemidefiniteOptimization, LeastSquares,SecondOrderConeOptimization,QuadraticOptimization,LinearOptimization,LinearFractionalOptimization,LinearModelFit, RobustConvexOptimization, NonlinearModelFit, FindFormula, GeometricOptimization, BayesianMinimization}]

A typical optimization example:

Design a minimum-weight cantilever beam subject to a vertical load F on the free end. The beam consists of n segments of unit length. Each segment has a width (Subscript[w, i]) and height (Subscript[h, i]):

n=16;
weight=w.h;

The load causes the beam to deflect and induces stress (Subscript[σ, i])=(6i F)/((Subscript[w, i])((Subscript[h, i])^2)), i=1,2,…,n on each segment of the beam. The stress must be less than the maximum allowable stress (Subscript[σ, max]):

stressConstraints=Table[6*i*F/(Indexed[w,i]*Indexed[h,i]^2)≤(Subscript[σ, max]),{i,n}];

A typical optimization example:

Let (Subscript[y, i]) be the deflection at the right end of segment i. The deflection (Subscript[y, 1]) at the free end can be found recursively as (Subscript[y, i])=(Subscript[y, i+1])+6(i-1/3) F/(Y (Subscript[w, i])((Subscript[h, i])^3))+(Subscript[v, i+1]), where Y is the Young's modulus and the slope (Subscript[v, i])=12(i-1/2)F/(Y (Subscript[w, i])((Subscript[h, i])^3))+(Subscript[v, i+1]). Since segment n is fixed, (Subscript[y, n])=0,(Subscript[v, n])=0:

v[n+1]=y=0;
Do[v[i]=(12(i-1/2)*(F/(Y*Indexed[w,i]*Indexed[h,i]^3)))+v[i+1],{i,n,1,-1}];
Do[y+=6*(i-1/3)*(F/(Y*Indexed[w,i]*Indexed[h,i]^3))+v[i+1],{i,n,1,-1}];

The deflection on the free end of the beam must be less than the maximum allowable deflection (Subscript[y, max]):

deflectionConstraint=y≤(Subscript[y, max]);

The width, height and aspect ratios are constrained to be in a certain range:

dimensionConstraints={0.1\[VectorLessEqual]w\[VectorLessEqual]20,0.1\[VectorLessEqual]h\[VectorLessEqual]6,0.1\[VectorLessEqual]h/w\[VectorLessEqual]5};

A typical optimization example:

Specify the parameters:

params={Y->1,F->0.1,(Subscript[σ, max])->1,(Subscript[y, max])->10};

Find the dimensions of the individual segments by minimizing the weight:

AbsoluteTiming[res=NMinimize[{weight,{stressConstraints,
deflectionConstraint,dimensionConstraints}}/.params,{w∈Vectors[n,Reals],h∈Vectors[n,Reals]}] ]

Visualize the cantilever beam:

beamSegments=MapThread[Translate[Cuboid[{-1/2,-#1/2,-#2/2},
{1/2,#1/2,#2/2}],{-#3,0,0}]&,{w,h,Range[0,n-1]}/.res[[2]]];
Graphics3D[{Opacity[0.5],beamSegments},BoxRatios->{3,1,1}]

A fitting example: Basis Pursuit

Use LASSO (L1) regularization to find a sparse fit (basis pursuit):

{σ,nτ,times,nω}={0.05,500,Subdivide[0.,1.,nτ],30};

Here is a signal:

signal = (1+.5 Sin[11 times])Sin[30 Sin[5 times]];
sdata = Transpose[{times, signal}];
ListLinePlot[sdata]

The goal is to approximate the signal with just of few of the thousands of Gabor basis functions:

Length[basis = Flatten[Table[Exp[-((((((t-τ))/σ)))^2)]
{1, Table[{Cos[ω t],Sin[ω t]},{ω,Drop[Subdivide[0.,150.,nω],1]}]},
{τ,Subdivide[0.,1.,nτ]}]]]

Plot[basis[[10000]],{t,0,1},PlotRange->All]

A fitting example: Basis Pursuit

With λ=1.`, a fit is found using a small number of basis elements:

AbsoluteTiming[sparseFit = Fit[sdata, basis, t,"BestFitParameters",FitRegularization->{"LASSO", 1}]]

The error is quite small:

fitt = basis.sparseFit /. t->times;
ListLinePlot[Transpose[{times, signal - fitt}]]

Once the important elements of the basis have been found, error can be reduced by finding the least-squares fit to these elements:

bpos = Flatten[sparseFit["NonzeroPositions"]];
sbasis = basis[[bpos]];
sbfit = Fit[sdata, sbasis, t];
ListLinePlot[Transpose[{times, signal-(sbfit /. t->times)}],PlotRange->All]

Exploring beyond the obvious

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];

We can solve this directly:

AbsoluteTiming[res=NMinimize[{r,Table[Norm[Inactive[Plus][p,-c]]<=r,{p,pts}]},{r,Element[c,Vectors[3,Reals]]}]]

Visualize the enclosing ball:

Show[Graphics3D[{Opacity[0.2],Red,Ball[c,r]}/.res[[2]]], region]

But....What points touch the ball??

{rr,cc}={r,c}/.res[[2]];
Flatten[Position[Table[Abs[Norm[p-cc]-rr]<=10^-5,{p,pts}],True]]

Exploring beyond the obvious: Reformulation

The same problem can be reformulated as 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};
AbsoluteTiming[Short[{minval,weights}=NMaximize[{obj,constraints},w∈Vectors[Length[pts]]],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:

{radius,center}={Sqrt[minval],(w/.weights).pts}

Only the extreme points have non-trivial weights:

Flatten[Position[Abs[w/.weights],x_/;x>10^-5]]

Exploring beyond the obvious: Use hidden information

Sometimes very useful answers and queries are hidden in plane sight:

{primal, dual}=ConicOptimization[r,Table[Norm[Inactive[Plus][p,-c]]<=r,{p,pts}],{r,Element[c,Vectors[3,Reals]]},{"PrimalMinimizerRules","DualMaximizer"}];
Short[dual,3]

index=Flatten[Position[Thread[Map[Total,dual]==0.],False]]

Visualize the enclosing ball with the extreme points:

Show[Graphics3D[{{Opacity[0.2],Red,Ball[c,r]}/.res[[2]], {PointSize[0.05],Point[pts[[index]]]}}], region]

What is a dual and why should I care?

A general minimization problem is:

By using Lagrange multipliers, we can write the optimization problem as:

The Lagrangian dual function (or just the dual problem):

Exploring beyond the obvious: Interpreting 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=ConicOptimization[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,cones}=ConicOptimization[t,{p1∈(Subscript[Ω, 1]),p2∈(Subscript[Ω, 2]),Norm[p1-p2]≤t},{p1,p2,t},{"DualMaximizer","ConicConstraintConeSpecifications"}]

The hyperplane is constructed as:

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

Exploring beyond the obvious: Decision making

A company makes three products. Each product is made using a combination of four machines.

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

Each machine can run a max of 2400 min/week. 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;
profit = {45.,60.,50.}.x;

How does varying the maximum capacity change the profit?

pfun=ParametricConvexOptimization[-profit,{timeConstraint,x\[VectorLessEqual]{c1,c2,c3}},x,{c1,c2,c3},"PrimalMinimumValue"]

Exploring beyond the obvious: Decision making

Let the current capacity of each product be 100, 40 and 60 units respectively:

currentCapacityRules={c1->100,c2->40,c3->60};

Current profit is:

currentProfit=-pfun[c1,c2,c3]/.currentCapacityRules

Varying the capacity gives insight into state of profit:

currentCapacity={c1,c2,c3}/.currentCapacityRules;
(-pfun@@(currentCapacity+#))&/@{{10,0,0},{0,10,0},{0,0,10}}

Exploring beyond the obvious: Decision making

Taking derivatives of the parameters gives insight into the change of profit:

dpfun=D[pfun[c1,c2,c3],{{c1,c2,c3}}];
-(dpfun/.currentCapacityRules).{0,0,10}

Manually confirm it:

pfun@@{100,40,60} - pfun@@{100,40,60+10}

This property is also HIDDEN inside the dual:

TakeDrop[Flatten[ConvexOptimization[-profit,{timeConstraint,
x\[VectorLessEqual]{100,40,60}},x,"DualMaximizer"]],4]

(dpfun/.currentCapacityRules)

Optimization with Uncertainty.

Consider a linear optimization problem:

Clear[x,y];
ConvexOptimization[-(x+y),{0≤x≤1,0≤y≤1},{x,y}]

What if we are not sure about the upper bounds on x and y? How do we find the best possible solution taking that uncertainty into account?

RegionPlot[{0≤x≤1&&0≤y≤1,0.9≤x≤1.1,0.9≤y≤1.1},{x,-0.1,1.1},{y,-0.1,1.1}, Epilog->{{PointSize[0.04],Point[{1,1}]}}]

Optimization with Uncertainty.

Taking uncertainty into account we have:

constraints = {0≤x≤(1+(Subscript[δ, x])),0≤y≤(1+(Subscript[δ, y]))};
uncertainty={-0.1≤(Subscript[δ, x])≤0.1,-0.1≤(Subscript[δ, y])≤0.1};

We want to find an optimal solution ((SuperStar[x]),(SuperStar[y])) that will be valid for all values of (Subscript[δ, x]),(Subscript[δ, y]).

One approach can be to sample (Subscript[δ, x]),(Subscript[δ, y]) from the uncertainty set and generate an expanded set of constraints.

newConstraints=Table[constraints /. {(Subscript[δ, x])->dx,(Subscript[δ, y])->dy},{dx,-0.1,0.1,0.01},{dy,-0.1,0.1,0.01}];
ConvexOptimization[-(x+y),newConstraints,{x,y}]

This is HIGHLY inefficient and sub-optimal!!!

Optimization with Uncertainty.

Fit a cubic curve to discrete data whose output contains interval uncertainty ±0.1`:

{input,output} = ;
dataplot=ListPlot[Transpose[{input, Interval[{-0.1,0.1}]+ output}]]

Construct the input matrix using DesignMatrix:

basis={1,x,(x^2),(x^3)};
mat=DesignMatrix[Transpose[{input,output}],basis,x];

Optimization with Uncertainty.

The objective is to minimize the error between the input and output subject to the interval uncertainty α:

objective=(InactiveHead[Norm])[mat.coeffs+-(output+α)];

Find the coefficients associated with each of the basis:

res=RobustConvexOptimization[objective,ForAll[α,{-0.1\[VectorLessEqual]α\[VectorLessEqual]0.1,α∈Vectors[Length[output],Reals]}],coeffs]

Find fit models for 200 possible realizations of the uncertainty:

n=Length[output];
pfun = ParametricConvexOptimization [objective, {},coeffs, α∈Vectors[n,Reals]];
models = Table[pfun[RandomReal[{-.1,.1},n]],{200}];

Show the envelope of the fit models with the robust fit (black) and the data:

Show[Plot[{(coeffs/.res).basis,Min[(coeffs/.models).basis],Max[(coeffs/.models).basis]},],dataplot]

Multiobjective Optimization - sneak peek

Minimize {(x^2)+(y^2), (((x-1))^2)+(y^2)} subject to the constraints -2<=x<=2,-2<=y<=2:

Multiobjective Optimization - sneak peek

Find the distribution of capital to invest in six stocks to maximize return while minimizing risk:

{μ,Σ}=;

The objective is to get a distribution of stock investment that gives a good balance between maximizing return while simultaneously minimizing risk:

objectives=Association[{"Minimize"->w.Σ.w, "Maximize"->μ.w}];

The percentage of investment w must be greater than 0 and must add to 1:

weightConstraints={w\[VectorGreaterEqual]0,Total[w]==1};

Compute the Pareto information:

res=MultiobjectiveOptimization[objectives,weightConstraints,w∈Vectors[Length[μ],Reals],Method->{Automatic,"MinSamplePoints"->20}]

Multiobjective Optimization - sneak peek

A desirable solution that balances risk vs return can be a solution that produces a point on the Pareto curve closest to the utopia point. Plot the Pareto front curve along with some critical points:

Show[res["ParetoFront","MeshRegion"],ListPlot[{{res["UtopiaObjectiveVector"]},
{res["ParetoFront","IdealUnscaledObjectiveVector"]}}, ],]

The weight distribution associated with the nearest point to the the utopia point is:

res["ParetoFront","IdealDecisionVariableRules"]

Find the point on the Pareto front where the return is 20.6:

res["ParetoFront",{"Intersection",{_,20.6}}]

Find the point on the Pareto front where the risk is 1.5:

res["ParetoFront",{"Intersection",{1.5,_}}]

Thank you!