Message Boards Message Boards

5
|
10523 Views
|
5 Replies
|
5 Total Likes
View groups...
Share
Share this post:

Solve the Karush-Kuhn-Tucker equations with Reduce

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}}}}
POSTED BY: Frank Kampas
5 Replies

Frank,

You should consider putting this code up in the new Wolfram function repository! It's a great post!

Regards,

Neil

POSTED BY: Neil Singer

I believe that the complexity of Reduce is double exponential. In other words, the time required to do a calculation increases very rapdily with the number of variables and constraints.

POSTED BY: Frank Kampas

Thanks again. I'll try to run it in a more powerful machine. Happy New Year!

Patricio

POSTED BY: Patricio García

Hi

Thanks 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 again

Patricio

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]
POSTED BY: Patricio García
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract