Some years ago I published a short article in the Mathematica Journal describing solving the Karush-Kuhn-Tucker equations with Reduce, to do symbolic optimization. I was pleased to see that the approach subsequently used by several people. However, the code in that article has the problem that it gives all local minima. I've recently updated the code to only give global minima. The new code has the advantage over Minimize that it gives multiple global minima and also provides the values of the Lagrange multipliers, which give the sensitivity of the objective function to changes in the constraints. The code is shown below with copious comments. I've also given two examples in which the code returns a result but Minimize does not, even though this is an unusual circumstance.
Code
In[1]:= (* Generate the Karush-Kuhn-Tucker Equations *)
KTEqs[obj_ (* objective function *), cons_List (* constraints *), vars_List (*
variables *)] :=
Module[{consconvrule = {GreaterEqual[x_, y_] -> LessEqual[y - x, 0],
Equal[x_, y_] -> Equal[x - y, 0],
LessEqual[x_, y_] -> LessEqual[x - y, 0],
LessEqual[lb_, x_, ub_] -> LessEqual[(x - lb) (x - ub), 0],
GreaterEqual[ub_, x_, lb_] -> LessEqual[(x - lb) (x - ub), 0]} ,
x, y, lb, ub , stdcons, eqcons, ineqcons, lambdas, mus, lagrangian, eqs1,
eqs2, eqs3, alleqns, allvars },
(* Change constraints to Equal and LessEqual form with zero on the right-
hand side *)
stdcons = cons /. consconvrule;
(* Separate the equality constraints and the inequality constraints *)
eqcons = Cases[stdcons, Equal[_, 0]][[All, 1]];
ineqcons = Ctdcons, LessEqual[_, 0]][[All, 1]];
(* Define the Lagrange multipliers for the equality and inequality \
constraints *)
lambdas = Array[\[Lambda], Length[eqcons]];
mus = Array[\[Mu], Length[ineqcons]];
(* Define the Lagrangian *)
lagrangian = obj + lambdas.eqcons + mus.ineqcons;
(* The derivatives of the Lagrangian are equal to zero *)
eqs1 = Thread[ D[lagrangian, {vars}] == 0];
(* Lagrange multipliers for inequality constraints are \[GreaterEqual]0 to \
get minima *)
eqs2 = Thread[mus >= 0];
(* Lagrange multipliers for inequality constraints are 0 unless the \
constraint value is 0 *)
eqs3 = Thread[mus*ineqcons == 0];
(* Collect the equations *)
alleqns = Join[eqs1, eqs2, eqs3, cons];
(* Collect the variables *)
allvars = Join[vars, lambdas, mus];
(* Return the equations and the variables *)
{alleqns, allvars}
]
In[2]:= (* Convert logical expressions to rules *)
torules[res_] := If[Head[res] === And, ToRules[res], List @@ (ToRules /@ res)]
In[3]:= (* Find the global minima *)
KKTReduce[obj_(* objective function *), cons_List (* constraints *),
vars_List (* variables *)] :=
Block[{kkteqs, kktvars, red, rls, objs, allres, minobj, sel, ret, minred,
minredrls},
(* Construct the equations and the variables *)
{kkteqs, kktvars} = KTEqs[obj, cons, vars];
(* Reduce the equations *)
red = LogicalExpand @
Reduce[kkteqs, kktvars, Reals, Backsubstitution -> True];
(* Convert the Reduce results to rules (if possible ) *)
rls = torules[red];
(* If the conversion to rules was complete *)
If[Length[Position[rls, _ToRules]] == 0,
(* Calculate the values of the objective function *)
objs = obj /. rls;
(* Combine the objective function values with the rules *)
allres = Thread[{objs, rls}];
(* Find the minimum objective value *)
minobj = Min[objs];
(* Select the results with the minimum objective value *)
sel = Select[allres, #[[1]] == minobj &];
(* Return the minimum objective value with the corresponding rules *)
ret = {minobj, sel[[All, 2]]},
(* Else if the results were not completely converted to rules *)
(* Use MinValue to find the smallest objective function value *)
minobj = MinValue[{obj, red}, kktvars];
(* Use Reduce to find the corresponding results *)
minred =
Reduce[obj == minobj && red, kktvars, Reals, Backsubstitution -> True];
(* Convert results to rules, if possible *)
minredrls = torules[minred];
ret = If[
Length[Position[minredrls, _ToRules]] == 0, {minobj, minredrls}, {minobj,
minred}];
];
(* Remove excess nesting from result *)
If[Length[ret[[2]]] == 1 && Depth[ret[[2]]] > 1, {ret[[1]], ret[[2, 1]]},
ret]
]
In[4]:=
Examples
In[5]:= Minimize[{x^2 - y^2, Cos[x - y] >= 1/2, -5 <= x <= 5, -5 <= y <= 5}, {x, y}]
Out[5]= Minimize[{x^2 - y^2, Cos[x - y] >= 1/2, -5 <= x <= 5, -5 <= y <= 5}, {x, y}]
In[6]:= KKTReduce[x^2 - y^2, {Cos[x - y] >= 1/2, -5 <= x <= 5, -5 <= y <= 5}, {x, y}]
Out[6]= {-25 + 25/9 (-3 + \[Pi])^2, {{x -> -(5/3) (-3 + \[Pi]),
y -> 5, \[Mu][1] -> (20 (-3 + \[Pi]))/(3 Sqrt[3]), \[Mu][2] ->
0, \[Mu][3] ->
1/9 (9 + 6 Sqrt[3] Sin[5 + 5/3 (-3 + \[Pi])] -
2 Sqrt[3] \[Pi] Sin[5 + 5/3 (-3 + \[Pi])])}, {x -> 5/3 (-3 + \[Pi]),
y -> -5, \[Mu][1] -> (20 (-3 + \[Pi]))/(3 Sqrt[3]), \[Mu][2] ->
0, \[Mu][3] ->
1/9 (9 + 6 Sqrt[3] Sin[5 + 5/3 (-3 + \[Pi])] -
2 Sqrt[3] \[Pi] Sin[5 + 5/3 (-3 + \[Pi])])}}}
In[7]:= TimeConstrained[
Minimize[{(Subscript[x, 1] - Subscript[x, 2])^2 + (Subscript[x, 2] -
Subscript[x, 3])^4, (1 + Subscript[x, 2]^2) Subscript[x, 1] + Subscript[
x, 3]^4 - 3 == 0}, {Subscript[x, 1], Subscript[x, 2], Subscript[x,
3]}], 60]
Out[7]= $Aborted
In[8]:= AbsoluteTiming @
KKTReduce[(Subscript[x, 1] - Subscript[x, 2])^2 + (Subscript[x, 2] -
Subscript[x, 3])^4, {(1 + Subscript[x, 2]^2) Subscript[x, 1] + Subscript[
x, 3]^4 - 3 == 0}, {Subscript[x, 1], Subscript[x, 2], Subscript[x, 3]}]
Out[8]= {1.67203, {0, {{Subscript[x, 1] -> 1, Subscript[x, 2] -> 1,
Subscript[x, 3] -> 1, \[Lambda][1] -> 0}, {Subscript[x, 1] ->
AlgebraicNumber[Root[3 + 2 #1 + 2 #1^2 + #1^3 &, 1], {0, 1, 0}],
Subscript[x, 2] ->
AlgebraicNumber[Root[3 + 2 #1 + 2 #1^2 + #1^3 &, 1], {0, 1, 0}],
Subscript[x, 3] ->
AlgebraicNumber[
Root[3 + 2 #1 + 2 #1^2 + #1^3 &, 1], {0, 1, 0}], \[Lambda][1] -> 0}}}}