New numerical optimization functionalities in Mathematica

Yifan Hu

Wolfram Research Inc

Abstract

In the talk we will discuss new optimization functionalities in Mathematica, including mixed integer linear programming, quadratic programming and fractional linear programming.

Outline of the talk

Constrained optimization function in Mathematica

New functionalities

Fractional programming

Integer linear programming

Quadratic programming

The constrained optimization problem

Constrained optimization problems in general:

Minimize f(x)
       subject to  h(x) == 0, g(x) >= 0

f --- the objective function

h, g --- the constraint functions, which together defines a feasible region, F.

If all functions involved are linear, the above gives a linear programming problem.

Introduction: local and global minimum

(SuperStar[x]) is a global minimum: f((SuperStar[x])) ≤  f(x) for all x ∈ F.

(SuperStar[x]) is a local minimum: f((SuperStar[x])) ≤  f(x) for all x ∈ N((SuperStar[x]))⋂ F.

Min      f = - (100/((((x-1))^2)+(((y-1))^2)+1))-(200/((((x+1))^2)+(((y+1))^2)+1))
subject to    -3≤x≤3, -3≤y≤3

Algorithms for optimization in Mathematica (Blue = New)

Blue = New

Mathematica has both functions for local and global optimization.

LinearProgramming: for linear programming and integer linear programming.

Algorithm: simplex/revised simplex, interior point. Branch and bound for integer LP.

Minimize: exact/arbitary precision global optimization. Integer linear programming. Linear and fractional programming. Parametric optimization for polynomial systems.

Algorithms: cylindrical decomposition. For equality constrained problems with box constraints, Lagrange multiplier method is used. It also solves integer linear programming problems using lattice reduction and branch and bound.

Algorithms for optimization in Mathematica (Blue = New)

NMinimize: global optimization.

Algorithms: direct search methods (Nelder-Mead; differential evolution; simulated annealing) and random search. Followed by gradient based local optimization. For linear and fractional linear programming, uses the same algorithm as FindMinimum.

FindMinimum/FindFit: local optimization. But can now handle constraints

Algorithms: Newton, quasi-Newton,  LevenbergMarquardt, conjugate gradient, principal-axis (for unconstrained), interior point ( for constrained), quadratic programming, fractional linear programming.

Fractional Linear programming

({{Minimize (((c^T)x + d)/((e^T)x + f))}, {(Subscript[A, ineq]) x ≥ (Subscript[b, ineq]), (Subscript[A, eq]) x== (Subscript[b, eq]), l ≤ x ≤ u}})

This problem can be tranformed into a linear programming problem in more than one way. The following is based on formulation of Adam Strzebonski.

Let z = 1/((e^T)x + f) and y = x/((e^T)x + f), then we have

({{Minimize (c^T)y + d z}, {s.t. (Subscript[A, ineq]) z ≥ (Subscript[b, ineq]) y, (Subscript[A, eq]) y== (Subscript[b, eq])z, l z≤ y ≤ u z, z≥ 0, (e^T)y + f z == 1}, {or (Subscript[A, ineq]) z ≤ (Subscript[b, ineq]) y, (Subscript[A, eq]) y== (Subscript[b, eq])z, l z≥ y ≥ u z, z≤ 0, (e^T)y + f z == 1}})

Fractional Linear programming

This is accessed through Minimize, NMinimize, or FindMinimum

Minimize[{((11+X[1]+2 X[2]+3 X[3]+4 X[4]+5 X[5]+6 X[6]+7 X[7]+8 X[8]+9 X[9]+10 X[10])/(111+101 X[1]+102 X[2]+103 X[3]+104 X[4]+105 X[5]+106 X[6]+107 X[7]+108 X[8]+109 X[9]+110X[10]
)),-X[1]≤0&&-X[2]≤0&&-X[3]≤0&&-X[4]≤0&&-X[5]≤0&&-X[6]≤0&&-X[7]≤0&&-X[8]≤0&&-X[9]≤0&&-X[10]≤0&&X[1]+X[2]+X[3]+X[4]+X[5]+X[6]+X[7]+X[8]+X[9]+X[10]==1},{X[1],X[2],X[3],X[4],X[5],X[6],X[7],X[8],X[9],X[10]}]

FindMinimum[{((11+X[1]+2 X[2]+3 X[3]+4 X[4]+5 X[5]+6 X[6]+7 X[7]+8 X[8]+9 X[9]+10 X[10])/(111+101 X[1]+102 X[2]+103 X[3]+104 X[4]+105 X[5]+106 X[6]+107 X[7]+108 X[8]+109 X[9]+110X[10]
)),-X[1]≤0&&-X[2]≤0&&-X[3]≤0&&-X[4]≤0&&-X[5]≤0&&-X[6]≤0&&-X[7]≤0&&-X[8]≤0&&-X[9]≤0&&-X[10]≤0&&X[1]+X[2]+X[3]+X[4]+X[5]+X[6]+X[7]+X[8]+X[9]+X[10]==1},{X[1],X[2],X[3],X[4],X[5],X[6],X[7],X[8],X[9],X[10]}]

NMinimize[{((11+X[1]+2 X[2]+3 X[3]+4 X[4]+5 X[5]+6 X[6]+7 X[7]+8 X[8]+9 X[9]+10 X[10])/(111+101 X[1]+102 X[2]+103 X[3]+104 X[4]+105 X[5]+106 X[6]+107 X[7]+108 X[8]+109 X[9]+110X[10]
)),-X[1]≤0&&-X[2]≤0&&-X[3]≤0&&-X[4]≤0&&-X[5]≤0&&-X[6]≤0&&-X[7]≤0&&-X[8]≤0&&-X[9]≤0&&-X[10]≤0&&X[1]+X[2]+X[3]+X[4]+X[5]+X[6]+X[7]+X[8]+X[9]+X[10]==1},{X[1],X[2],X[3],X[4],X[5],X[6],X[7],X[8],X[9],X[10]}]

Fractional Linear programming

What if z = 0?

Min ((Subscript[x, 1])/(1+(Subscript[x, 1])+(Subscript[x, 2]))), s.t.. (Subscript[x, 1])≥ 1, (Subscript[x, 1])≤ 0

transformed into z = 1/(1+(Subscript[x, 1])+(Subscript[x, 2])), y = {(Subscript[x, 1]),(Subscript[x, 2])}/(1+(Subscript[x, 1])+(Subscript[x, 2])).

Min (Subscript[y, 1]), s.t. ((Subscript[y, 1]) ≥ z , (Subscript[y, 1]) ≤ 0, z ≥ 0, (Subscript[y, 1])+(Subscript[y, 2])+z = 1) || ((Subscript[y, 1]) ≤ z , (Subscript[y, 1]) ≥ 0, z ≤ 0, (Subscript[y, 1])+(Subscript[y, 2])+z == 1)

The solution is z = 0, (Subscript[y, 1]) = 0, (Subscript[y, 2]) = 1 thus (Subscript[x, 1]) = (Subscript[y, 1])/z = Indeterminate, (Subscript[x, 2]) = Infinity, optimal value 0. But the problem is actually infeasible

So the algorithm always check if the original problem is feasible when the transformed problem gives z = 0.

Fractional linear programming application: data envelope analysis (DEA)

Consider a number of bank branches (http://people.brunel.ac.uk/~mastjjb/jeb/or/dea.html)

bankdata1={{"Branch",     "Personal transactions (`1000s)",
                    "Number of staff"},{"Croydon",125,18},{"Dorking",44,16},{"Redhill",80,17},{"Reigate",23,11}};
Grid[bankdata1,Dividers->All]

We can compare bank efficiency by a ratio

bankdata2={{"Branch",     "Personal transactions (`1000s) per staff","Efficiency"},{"Croydon",6.94,"100%"},{"Dorking",2.75,"40%"},{"Redhill",4.71,"68%"},{"Reigate",2.09,"30%"}};
Grid[bankdata2,Dividers->All]

Fractional linear programming application: data envelope analysis (DEA)

But what if there are more than one output?

bankdata3={{"Branch","Personal transactions  ('000s)","Business transactions ('000s)","staff"},{"Croydon",125, 50, 18},{"Dorking",44,20,16},{"Redhill",80,55,17},{"Reigate ",23,12,11}};
Grid[bankdata3,Dividers->All]

lab={"Croydon","Dorking","Redhill","Reigate"};GraphPlot[{1->1,2->2,3->3,4->4},SelfLoopStyle->False,VertexCoordinateRules->{{125,50}/18,{44,20}/16,{80,55}/17,{23,12}/11},AxesOrigin->{0,0},VertexRenderingFunction->({Point[#1],Text[lab[[#3]],#1+{1,0}]}&),Axes->True,AxesOrigin->{0,0},AspectRatio->1,AxesLabel->{"Personal transation/staff","Business transactions/staff"}]

Fractional linear programming application: data envelope analysis (DEA)

We can use a weight to make the "Personal transactions" and "Bussiness transaction" into one

What weights to choose? We choose weight differently when measuring efficiency of each DMU (decision-making units ), therefore presenting the DMU in the best possible light. This is DEA (Data Envelope Analysis, Charnes, Cooper and Rhodes, 1978)

Fractional linear programming application: data envelope analysis (DEA)

Example: efficiency of Dorking

Max (Subscript[E, dorking])
s.t. 
(Subscript[E, croydon ])= (125 (Subscript[W, per]) + 50(Subscript[W, bus]))/(18 (Subscript[W, staff]))
(Subscript[E, Dorking])=(44(Subscript[W, per]) +20(Subscript[W, bus]))/(16 (Subscript[W, staff]))
(Subscript[E, Redhill]) = (80(Subscript[W, per]) +55(Subscript[W, bus]))/(17 (Subscript[W, staff]))
(Subscript[E, Reigate]) = (23(Subscript[W, per]) +12(Subscript[W, bus]))/(11 (Subscript[W, staff]))
0≤ (Subscript[E, Croydon]) ≤ 1
0≤ (Subscript[E, Dorking]) ≤ 1
0≤ (Subscript[E, Redhill]) ≤ 1
0≤ (Subscript[E, Reigate]) ≤ 1
(Subscript[W, per]) ≥ 0
(Subscript[W, bus]) ≥ 0
(Subscript[W, staff]) ≥ 0

FindMaximum[{((44 (Subscript[W, per]) +20 (Subscript[W, bus]) )/(16 (Subscript[W, staff]))),{125 (Subscript[W, per]) + 50(Subscript[W, bus])<=18 (Subscript[W, staff]),44 (Subscript[W, per]) + 20(Subscript[W, bus])<=16 (Subscript[W, staff]),80 (Subscript[W, per]) + 55(Subscript[W, bus])<=17 (Subscript[W, staff]),23 (Subscript[W, per]) + 12(Subscript[W, bus])<=11 (Subscript[W, staff]),(Subscript[W, per])>=0, (Subscript[W, bus]) >= 0, (Subscript[W, staff])>=0}},{(Subscript[W, per]), (Subscript[W, bus]), (Subscript[W, staff])}]

FindMaximum[{((44 (Subscript[W, per]) +20 (Subscript[W, bus]) )/(16 (Subscript[W, staff]))),{125 (Subscript[W, per]) + 50(Subscript[W, bus])<=18 (Subscript[W, staff]),44 (Subscript[W, per]) + 20(Subscript[W, bus])<=16 (Subscript[W, staff]),80 (Subscript[W, per]) + 55(Subscript[W, bus])<=17 (Subscript[W, staff]),23 (Subscript[W, per]) + 12(Subscript[W, bus])<=11 (Subscript[W, staff]),(Subscript[W, per])+(Subscript[W, bus ])==1,(Subscript[W, per])>=0, (Subscript[W, bus]) >= 0, (Subscript[W, staff])>=0}},{(Subscript[W, per]), (Subscript[W, bus]), (Subscript[W, staff])}]

FindMaximum[{((23 (Subscript[W, per]) +12 (Subscript[W, bus]) )/(11 (Subscript[W, staff]))),{125 (Subscript[W, per]) + 50(Subscript[W, bus])<=18 (Subscript[W, staff]),44 (Subscript[W, per]) + 20(Subscript[W, bus])<=16 (Subscript[W, staff]),80 (Subscript[W, per]) + 55(Subscript[W, bus])<=17 (Subscript[W, staff]),23 (Subscript[W, per]) + 12(Subscript[W, bus])<=11 (Subscript[W, staff]),(Subscript[W, per])+(Subscript[W, bus ])==1,(Subscript[W, per])>=0, (Subscript[W, bus]) >= 0, (Subscript[W, staff])>=0}},{(Subscript[W, per]), (Subscript[W, bus]), (Subscript[W, staff])}]

bankdata3={{"Branch","Efficiency"},{"Croydon",1},{"Dorking",0.43},{"Redhill",1},{"Reigate ",0.36}};
Grid[bankdata3,Dividers->All]

Quadratic programming

In the new version, we will be able to solve quadratic programming more efficiently

Access is through FindMinimum using Method→"QuadraticProgramming"

FindMinimum[{x[1]^2/100+x[2]^2-100,{10*x[1]-x[2]≥10,2≤x[1]≤50,-50≤x[2]≤50}},{x[1],x[2]},Method->"QuadraticProgramming"]

Quadratic programming: an application example

Portfolio optimization: asset allocation

Find the optimal asset allocation so as to minimize the risk, and achieve a preset level of return.

It rely on assets that has little or inverse correlations. E.g., if asset A increases 20% one year and down 10% next, asset B down 10% one year and up 20% next,  and up years for A are down years for B, then holding both in equal amont would result in 10% every year.

Equities are sometimes uncorrelated to bonds, and probably has little correlation to gold/property. Hence a suitable mix of them would reduce the risk while still provide good return.

Data: return of assets between 1973 to 1994

R = {{1.075`,1.084`,1.061`,1.052`,1.055`,1.077`,1.109`,1.127`,1.156`,1.117`,1.092`,1.103`,1.08`,1.063`,1.061`,1.071`,1.087`,1.08`,1.057`,1.036`,1.031`,1.045`},{0.942`,1.02`,1.056`,1.175`,1.002`,0.982`,0.978`,0.947`,1.003`,1.465`,0.985`,1.159`,1.366`,1.309`,0.925`,1.086`,1.212`,1.054`,1.193`,1.079`,1.217`,0.889`},{0.852`,0.735`,1.371`,1.236`,0.926`,1.064`,1.184`,1.323`,0.949`,1.215`,1.224`,1.061`,1.316`,1.186`,1.052`,1.165`,1.316`,0.968`,1.304`,1.076`,1.1`,1.012`},{0.815`,0.716`,1.385`,1.266`,0.974`,1.093`,1.256`,1.337`,0.963`,1.187`,1.235`,1.03`,1.326`,1.161`,1.023`,1.179`,1.292`,0.938`,1.342`,1.09`,1.113`,0.999`},{0.698`,0.662`,1.318`,1.28`,1.093`,1.146`,1.307`,1.367`,0.99`,1.213`,1.217`,0.903`,1.333`,1.086`,0.959`,1.165`,1.204`,0.83`,1.594`,1.174`,1.162`,0.968`},{1.023`,1.002`,0.123`,1.156`,1.03`,1.012`,1.023`,1.031`,1.073`,1.311`,1.08`,1.15`,1.213`,1.156`,1.023`,1.076`,1.142`,1.083`,1.161`,1.076`,1.11`,0.965`},{0.851`,0.768`,1.354`,1.025`,1.181`,1.326`,1.048`,1.226`,0.977`,0.981`,1.237`,1.074`,1.562`,1.694`,1.246`,1.283`,1.105`,0.766`,1.121`,0.878`,1.326`,1.078`},{1.677`,1.722`,0.76`,0.96`,1.2`,1.295`,2.212`,1.296`,0.688`,1.084`,0.872`,0.825`,1.006`,1.216`,1.244`,0.861`,0.977`,0.922`,0.958`,0.926`,1.146`,0.99`}};

({{, "3m Tbill", "long Tbond", "SP500", "Wilt.5000", "Corp. Bond", "NASDQ", "EAFE", "Gold"}, {1973, 1.075`, 0.942`, 0.852`, 0.815`, 0.698`, 1.023`, 0.851`, 1.677`}, {1974, 1.084`, 1.02`, 0.735`, 0.716`, 0.662`, 1.002`, 0.768`, 1.722`}, {1975, 1.061`, 1.056`, 1.371`, 1.385`, 1.318`, 0.123`, 1.354`, 0.76`}, {1976, 1.052`, 1.175`, 1.236`, 1.266`, 1.28`, 1.156`, 1.025`, 0.96`}, {1977, 1.055`, 1.002`, 0.926`, 0.974`, 1.093`, 1.03`, 1.181`, 1.2`}, {1978, 1.077`, 0.982`, 1.064`, 1.093`, 1.146`, 1.012`, 1.326`, 1.295`}, {1979, 1.109`, 0.978`, 1.184`, 1.256`, 1.307`, 1.023`, 1.048`, 2.212`}, {1980, 1.127`, 0.947`, 1.323`, 1.337`, 1.367`, 1.031`, 1.226`, 1.296`}, {1981, 1.156`, 1.003`, 0.949`, 0.963`, 0.99`, 1.073`, 0.977`, 0.688`}, {1982, 1.117`, 1.465`, 1.215`, 1.187`, 1.213`, 1.311`, 0.981`, 1.084`}, {1983, 1.092`, 0.985`, 1.224`, 1.235`, 1.217`, 1.08`, 1.237`, 0.872`}, {1984, 1.103`, 1.159`, 1.061`, 1.03`, 0.903`, 1.15`, 1.074`, 0.825`}, {1985, 1.08`, 1.366`, 1.316`, 1.326`, 1.333`, 1.213`, 1.562`, 1.006`}, {1986, 1.063`, 1.309`, 1.186`, 1.161`, 1.086`, 1.156`, 1.694`, 1.216`}, {1987, 1.061`, 0.925`, 1.052`, 1.023`, 0.959`, 1.023`, 1.246`, 1.244`}, {1988, 1.071`, 1.086`, 1.165`, 1.179`, 1.165`, 1.076`, 1.283`, 0.861`}, {1989, 1.087`, 1.212`, 1.316`, 1.292`, 1.204`, 1.142`, 1.105`, 0.977`}, {1990, 1.08`, 1.054`, 0.968`, 0.938`, 0.83`, 1.083`, 0.766`, 0.922`}, {1991, 1.057`, 1.193`, 1.304`, 1.342`, 1.594`, 1.161`, 1.121`, 0.958`}, {1992, 1.036`, 1.079`, 1.076`, 1.09`, 1.174`, 1.076`, 0.878`, 0.926`}, {1993, 1.031`, 1.217`, 1.1`, 1.113`, 1.162`, 1.11`, 1.326`, 1.146`}, {1994, 1.045`, 0.889`, 1.012`, 0.999`, 0.968`, 0.965`, 1.078`, 0.99`}})

data

TableForm[Transpose[R],TableHeadings->{{1973,1974,1975,1976,1977,1978,1979,1980,1981,1982,1983,1984,1985,1986,1987,1988,1989,1990,1991,1992,1993,1994},{"3m Tbill","long Tbond","SP500","Wilt.5000","Corp. Bond","NASDQ","EAFE","Gold"}}]

Quadratic programming: portfolio optimization

Expected return

{n,nyear}=Dimensions[R];
(* expected return for {"3m Tbill","long Tbond","SP500","Wilt.5000","Corp. Bond","NASDQ","EAFE","Gold"}*)
ER = Map[Total,R]/nyear

The covariant matrix

Covariants=(R-ER).Transpose[R-ER];
Covariants//MatrixForm

Optimal asset allocation: minimize the volitility and achieve 12% return

vars = Map[Subscript[x,#]&,{"3m T-bill","long T-bond","SP500","Wiltshire 5000","Corporate Bond","NASDQ","EAFE","Gold"}];
FindMinimum[{
vars.Covariants.vars,
Total[vars]==1&&
vars.ER≥1.12&&
Apply[And,Thread[Greater[vars,0]]]
},vars]//Timing

FindMinimum[{
vars.Covariants.vars,
Total[vars]==1&&
vars.ER≥1.12&&
Apply[And,Thread[Greater[vars,0]]]
},vars,Method->"QuadraticProgramming"]//Timing

Quadratic programming: portfolio optimization

10% return

FindMinimum[{
vars.Covariants.vars,
Total[vars]==1&&
vars.ER≥1.10&&
Apply[And,Thread[Greater[vars,0]]]
},vars]//Timing

FindMinimum[{
vars.Covariants.vars,
Total[vars]==1&&
vars.ER≥1.10&&
Apply[And,Thread[Greater[vars,0]]]
},vars,Method->"QuadraticProgramming"]//Timing

Larger problems: 100 assets with 50 years of historical data, randomly generated

{n,nyear}={100,50};
SeedRandom[1];
R=Table[Random[Real,{0.7,1.4}],{n},{nyear}];
ER = Map[Total,R]/nyear;
Covariants=(R-ER).Transpose[R-ER];
vars = Array[x,n];
FindMinimum[{
vars.Covariants.vars,
Total[vars]==1&&
vars.ER≥1.08&&
Apply[And,Thread[Greater[vars,0]]]
},vars]//Timing

FindMinimum[{
vars.Covariants.vars,
Total[vars]==1&&
vars.ER≥1.08&&
Apply[And,Thread[Greater[vars,0]]]
},vars,Method->"QuadraticProgramming"]//Timing

Integer linear programming

Pre-6.0, Minimize is able to solve small pure integer linear programming problems (not mixed) using branch and bound, combined with lattice reduction. NMinimize was also able to solve small integer/mixed integer programming problems.

Minimize[{x+y,x>0&&y>0&&x+2y>3},{x,y}]

Minimize[{x+y,x≥0&&y>0&&x+2 y>3&&{x,y}∈Integers},{x,y}]

Integer/mixed integer linear programming

In the coming version, (machine) integer and mixed integer linear programming problems will be solved much more efficiently

Through LinearProgramming, FindMinimum and NMinimize

FindMinimum[{x+y,x>0&&y>0&&x+2y>3&&{x,y}∈Integers},{x,y}]

LinearProgramming[{1,1},{{1,2}},{3.},{Integers,Integers}]

NMinimize[{x+y,x>0&&y>0&&x+2y>3&&{x,y}∈Integers},{x,y}]

Integer linear programming

A nontrivial problem

Import["flugpl.mps","MPS"]

FindMinimum@@%

Integer linear programming: a larger problem

See what elements this data contain

Import["mod010.mps",{"MPS","Elements"}]

Check how many constrains/variables

Import["misc06.mps",{"MPS","ConstraintMatrix"}]

Solve it

FindMinimum@@Import["misc06.mps","MPS"];//Timing

Aplication: sudoku

This is a code to solve Sudoku problem (from Dan Lichtblau)

sudoku[sparseArray[data_,d_],diags_: False]:=Module[{x,vars,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,constraints,known,soln,n=d^3,sd=Sqrt[d],newvars,j1,k1,removables,augment,newknown,others,fvars,PrintSudoku},
PrintSudoku[string_,soln_]:= Print[string,Grid[Normal[soln]/.{0->""},Dividers->{{{True,False,False}},{{True,False,False}}}]];

vars=Array[x,{d,d,d}];
PrintSudoku["Input=",SparseArray[data,{d,d}]];
fvars=Flatten[vars];
c2=Table[Total[vars[[All,j,k]]]==1,{j,d},{k,d}];
c3=Table[Total[vars[[j,All,k]]]==1,{j,d},{k,d}];
c4=Table[Total[vars[[j,k,All]]]==1,{j,d},{k,d}];
c5=Table[Sum[vars[[j,k+k2*sd,m+m2*sd]],{k,sd},{m,sd}]==1,{k2,0,sd-1},{m2,0,sd-1},{j,d}];
c6=If[diags,{Table[Sum[vars[[k,j,j]],{j,d}]==1,{k,d}],Table[Sum[vars[[k,j,d-j+1]],{j,d}]==1,{k,d}]},(*else*){}];
c7=Map[Table[x[j,#[[1,1]],#[[1,2]]]->If[j==#[[2]],1,0],{j,d}]&,data];
c8=Map[Table[x[#[[2]],j,#[[1,2]]]->If[j==#[[1,1]],1,0],{j,d}]&,data];
c9=Map[Table[x[#[[2]],#[[1,1]],j]->If[j==#[[1,2]],1,0],{j,d}]&,data];
c10=Map[(j1=sd*Quotient[#[[1,1]]-1,sd];
k1=sd*Quotient[#[[1,2]]-1,sd];
Table[x[#[[2]],j,k]->If[j==#[[1,1]]&&k==#[[1,2]],1,0],{j,j1+1,j1+sd},{k,k1+1,k1+sd}])&,data];
known=Dispatch[Union[Flatten[{c7,c8,c9,c10}]]];
constraints=Flatten[({c2,c3,c4,c5,c6}/.known)/.True->Sequence[]];
newvars=Complement[fvars,Variables[known[[1]]/.Rule->List]];
(*reduction to core constraints loop*)removables=Cases[constraints,x[_,_,_]==1];
augment={};
While[removables=!={},newknown=removables/.Equal->Rule;
constraints=(constraints/.newknown)/.True->Sequence[];
others=Cases[constraints,1+_==1];
If[others=!={},others=Thread[Variables[others/.Equal->List]->0];
constraints=(constraints/.others)/.True->Sequence[];];
removables=Cases[constraints,x[_,_,_]==1];
augment={others,newknown,augment};];
If[augment=!={},newvars=Complement[newvars,Variables[Flatten[augment]/.Rule->List]];
known=Dispatch[Union[Flatten[{known[[1]],augment}]]];];
(*end core equality constraint loop*)If[Length[newvars]==0,soln={0,Table[vars[[All,j,k]].Range[d],{j,d},{k,d}]/.known};PrintSudoku["Solution=",soln[[2]]];Return[soln]];
constraints=Join[constraints,{Element[newvars,Integers]},Map[0≤#≤1&,newvars]];
soln=FindMinimum[{1,constraints},newvars];
soln=soln[[2]]/.aa_Real:>Round[aa];
soln=Table[vars[[All,j,k]].Range[d],{j,d},{k,d}]/.soln/.known;
If [!MatrixQ[soln],soln=soln[[2]]];
PrintSudoku["Solution=",soln];

soln];

sudoku problem ({{"", 6, "", 1, "", 4, "", 5, ""}, {"", "", 8, 3, "", 5, 6, "", ""}, {2, "", "", "", "", "", "", "", 1}, {8, "", "", 4, "", 7, "", "", 6}, {"", "", 6, "", "", "", 3, "", ""}, {7, "", "", 9, "", 1, "", "", 4}, {5, "", "", "", "", "", "", "", 2}, {"", "", 7, 2, "", 6, 9, "", ""}, {"", 4, "", 5, "", 8, "", 7, ""}})

mat1=sparseArray[{{1,2}->3,{1,7}->9,{2,3}->6,{3,4}->2,{3,5}->4,{3,6}->1,{3,8}->3,{4,4}->9,{4,7}->7,{5,6}->2,{5,9}->4,{6,2}->8,{6,5}->7,{6,8}->2,{7,1}->8,{7,2}->5,{8,2}->9,{8,4}->7,{8,6}->4,{9,6}->6,{9,9}->1},9];
Timing[res=sudoku[mat2];]

mat2=sparseArray[{{1,2}->6,{1,4}->1,{1,6}->4,{1,8}->5,{2,3}->8,{2,4}->3,{2,6}->5,{2,7}->6,{3,1}->2,{3,9}->1,{4,1}->8,{4,4}->4,{4,6}->7,{4,9}->6,{5,3}->6,{5,7}->3,{6,1}->7,{6,4}->9,{6,6}->1,{6,9}->4,{7,1}->5,{7,9}->2,{8,3}->7,{8,4}->2,{8,6}->6,{8,7}->9,{9,2}->4,{9,4}->5,{9,6}->8,{9,8}->7},9];
Timing[res=sudoku[mat2];]

mat3=sparseArray[{{1,2}->10,{1,5}->16,{1,8}->14,{1,13}->15,{1,14}->13,{2,2}->11,{2,3}->6,{2,5}->2,{2,7}->1,{2,10}->4,{2,14}->9,{3,12}->12,{3,16}->7,{4,1}->13,{4,6}->10,{4,11}->15,{4,13}->5,{4,14}->3,{4,15}->8,{5,4}->8,{5,5}->6,{5,8}->11,{5,9}->4,{5,10}->3,{5,12}->14,{6,2}->15,{6,7}->8,{6,8}->2,{6,9}->13,{6,13}->14,{6,14}->1,{6,15}->4,{7,3}->16,{7,4}->7,{7,6}->5,{7,10}->1,{7,14}->2,{7,16}->13,{8,2}->14,{8,5}->12,{8,9}->16,{8,10}->9,{8,11}->10,{9,6}->13,{9,7}->7,{9,8}->3,{9,12}->4,{9,15}->11,{10,1}->5,{10,3}->12,{10,7}->15,{10,11}->8,{10,13}->10,{10,14}->14,{11,2}->2,{11,3}->9,{11,4}->14,{11,8}->10,{11,9}->3,{11,10}->16,{11,15}->12,{12,5}->14,{12,7}->4,{12,8}->5,{12,9}->1,{12,12}->6,{12,13}->16,{13,2}->16,{13,3}->13,{13,4}->9,{13,6}->8,{13,11}->14,{13,16}->6,{14,1}->12,{14,5}->10,{15,3}->4,{15,7}->9,{15,10}->11,{15,12}->7,{15,14}->15,{15,15}->3,{16,3}->1,{16,4}->3,{16,9}->5,{16,12}->2,{16,15}->7},16];
Timing[res=sudoku[mat3];]

Conclusions

New functionalities in the coming version, including 

     fractional linear programming
     
     quadratic programming
     
     integer/mixed integer linear programming