Message Boards Message Boards

GROUPS:

Assuming lambdas >= 0 for KKT script

Posted 2 months ago
554 Views
|
2 Replies
|
0 Total Likes
|

Hey everyone! I'm currently working on some scripts for different elements of optimization including KKT. I've just made a switch to be using Mathematica and

I found the following script from this forum https://community.wolfram.com/groups/-/m/t/193644

The problem though is, that it does not properly assume lambdas being >= which I need. I've tried inserting Assumptions -> when I'm using the application, but I don't think it correctly registers the lambda in that, as the lambdas are only being declared in module code.

Therefore, I would like help on how to properly force the application to assume all lambdas being >= 0

Thanks! - Nik

POSTED BY: Nik Aa
Answer
2 Replies

Hard to say much without having the actual code (in copy/paste format).

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}}}}
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