only the lambdas for inequality constraints should be >= 0
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 = Cases[stdcons, 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.21184, {0, {{Subscript[x, 1] -> 1, Subscript[x, 2] -> 1,
Subscript[x, 3] -> 1, \[Lambda][1] -> 0}, {Subscript[x, 1] ->
AlgebraicNumber[
Root[3 + 2 # + 2 #^2 + #^3& , 1, 0], {0, 1, 0}],
Subscript[x, 2] -> AlgebraicNumber[
Root[3 + 2 # + 2 #^2 + #^3& , 1, 0], {0, 1, 0}],
Subscript[x, 3] -> AlgebraicNumber[
Root[3 + 2 # + 2 #^2 + #^3& , 1, 0], {0, 1, 0}], \[Lambda][1] -> 0}}}}