# Solve the Karush-Kuhn-Tucker equations with Reduce

Posted 1 year ago
1775 Views
|
5 Replies
|
5 Total Likes
|
 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 = 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.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}}}} 
5 Replies
Sort By:
Posted 1 month ago
 Frank,You should consider putting this code up in the new Wolfram function repository! It's a great post!Regards,Neil
Posted 6 months ago
Posted 6 months ago
 Thanks again. I'll try to run it in a more powerful machine. Happy New Year!Patricio
 HiThanks so much for the code. I've tried on a trivial example and it's very fast. However, when I used it in a more demanding setting, had to abort it since it wouldn't return any result after over an hour running. I'm copying my problem below and I was hoping you cold tell me whether it is the nature of my problem or that I'm mistaken in the use of the code what is causing such a long execution time.Thanks againPatricio T = {{2, 3, 1}, {3, 2, 1}, {1, 4, 2}, {4, 1, 3}}; techniques = Dimensions[T][[1]]; αVec = ConstantArray[4/6, techniques]; qVec = Array[q, techniques]; onesVec = ConstantArray[1, techniques]; needs = Transpose[T].qVec^(1/αVec); FactorNeeds[mix_] := Transpose[T].mix^(1/αVec); qcritical = Table[qVec /. Minimize[{needs[[j]], onesVec.qVec == 1, Thread[1 >= qVec >= 0]}, qVec, Reals][[2]], {j, 1, factors}]; cons = Join[ Thread[needs <= FactorNeeds[Transpose[qcritical].{0, 1, 0}]], Thread[qVec >= 0]] KKTReduce[-qVec.onesVec,cons,qVec]