(* Name of this file : MutiplierMethod.ma *)
(* Examples of MultiplierMethod.m  use *)
(* The problem solved here is Min f  s.t. g = 0 andh <= 0*)
(* Note that f can be 0. and g or h can be empty lists {} *)

 (* Any comments welcome!   J.-C. Culioli, culioli@cas.ensmp.fr  *)

Needs["Optimization`MultiplierMethod`"];

f = x1^2 + (x2-4)^2 ; g ={x1-x2} ; h = {5 x2};
x = {x1,x2}; xk = {1,2}; MultiplierMethod[f,g,h,{x1,x2},{1.,2.}]

MultiplierMethod[f,g,h,{x1,x2},{1.,2.},
ConvergenceGoal -> -5, DualParameter -> True]

u = {u1,u2,u3}; g = {{1, 1,1},{0,1,-2}}. u + {5,-1}; h = {};
j = (u1-u2)^2 + u3;
MultiplierMethod[j,g,h,u,{1.,2.,0.},DualParameter -> True]

f = -x2; g =  { x1^2 + x2^2 - 1 }; h = {-x1, -1 - x1 + 2 x2};
MultiplierMethod[f,g,h,{x1,x2},{10.,10.},Augmentation -> 5]

(* Unconstrained optimization...*)
MultiplierMethod[50 (x1^2-x2)^2 + (x1-1)^2,{},{},{x1,x2},{-2.,3.}]

FindMinimum[50 (x1^2-x2)^2 + (x1-1)^2,{x1,-2},{x2,3}]

Other Tests

Options[MultiplierMethod]

Trivial 1

f = x1^2   ; g ={} ;
h = {-x1+1}; MultiplierMethod[f,g,h,{x1},{5.}] // Timing

Trivial 2

f = (x1-2)^2   ; g ={x1 -1} ;
h = { }; MultiplierMethod[f,g,h,{x1},{5.}] // Timing

Trivial 3

f = x1^2 + (x2-2)^2  ; g ={x1 + x2 -5/2} ;
h = {-x1 + .25 , -x2 + 2};
MultiplierMethod[f,g,h,{x1,x2},{5.,10.}] // Timing

Quadratic programming

f = x1^2 + (x2-4)^2 ; g ={x1-x2} ;
h = {5 x2}; MultiplierMethod[f,g,h,{x1,x2},{1.,2.}] // Timing

Quadratic programming /  inequality only

f = 2 x1^2 + x1 x2 + x2^2 - 12 x1 - 10 x2;
g = {};  h = {x1 + x2 - 4, -x1, -x2};
MultiplierMethod[f,g,h,{x1,x2},{0,0},
DualParameter -> True]  // Timing

Quadratic programming with 2 equalities and one inequality...

u = {u1,u2,u3}; j = (u1-u2)^2 + 5 u3;
g = {{1, 1,1},{0,1,-3}}. u + {5,-1};  h = {-u1};
MultiplierMethod[j,g,h,u,{1.,2.,0.},
DualParameter -> True]  // Timing

Tests from the books by  W. Hock & K. Schittkowski (1980) and  K. Schittkowski (1987)...

TP1 p. 24

f = 100 (x2 - x1^2)^2 + (1-x1)^2; g = {}; h = {-x2 -1.5};
MultiplierMethod[f,g,h,{x1,x2},{-2.,1.}] // Timing

TP2 p. 25

f = 100 (x2 - x1^2)^2 + (1-x1)^2; g = {}; h = {-x2 +1.5};
MultiplierMethod[f,g,h,{x1,x2},{-2.,1.}]  // Timing

TP216 p. 40

f = 100 (x1^2 - x2)^2 + (x1 - 1)^2;
g =  { x1 (x1-4) - 2 x2 + 12 };
h = {};
MultiplierMethod[f,g,h,{x1,x2},{-1.2,1.},
DualParameter -> True]  // Timing

TP217 p. 41

f = -x2;
g =  { x1^2 + x2^2 - 1 };
h = {-x1, -1 - x1 + 2 x2};   (* 28 s *)
MultiplierMethod[f,g,h,{x1,x2},{10.,10.}] // Timing

f = -x2;
g =  { x1^2 + x2^2 - 1 };
h = {-x1, -1 - x1 + 2 x2};   (* 32 s *)
MultiplierMethod[f,g,h,{x1,x2},{10.,10.},
ConvergenceGoal -> 3, Augmentation -> 1] // Timing

TP230 p. 54  Chamberlain Problem

f = x2;  g =  {};
h = - {-2 x1^2 + x1^3 + x2 , -2(1-x1)^2 + (1-x1)^3 + x2};

MultiplierMethod[f,g,h,{x1,x2},{0.,0.}]   // Timing

TP234 p. 58

f = (x2-x1)^4 -(1-x1); g =  {};
h = {x1 - 2, x2 - 2, -x1 + .2, -x2 + .2, x1^2 + x2^2 -1};

MultiplierMethod[f,g,h,{x1,x2},{0.,0.}] // Timing

TP321 p. 142      ...difficult...

f = (x1-20)^2 + (x2 + 20)^2;
g =  {x1^2 /100. + x2^2 -1}; h = {};
MultiplierMethod[f,g,h,{x1,x2},{0.,0.}] // Timing

f = (x1-20)^2 + (x2 + 20)^2;
g =  {x1^2 /100. + x2^2 -1}; h = {};
MultiplierMethod[f,g,h,{x1,x2},{0.,0.},OuterIterations->100] // Timing

TP322 p. 143    ...difficult...

f = (x1-20)^2 + (x2 + 20)^2;
g =  {x1^2 /100. + 100. x2^2 -1}; h = {};
MultiplierMethod[f,g,h,{x1,x2},{0.,0.},OuterIterations->150] // Timing

TP339  p. 160 Container Problem  (difficult...)
solution : {3.36168, {x1 -> 2.380, x2 -> 0.3162, x3 -> 1.943}}

f = 1/(5. x1 x2 x3) + 4./x1 + 3./x3;   (* 28 s *)
h =  {2. x1 x3 + x1 x2 -10., -x1,-x2,-x3,x1-5,x2-5,x3-5};
g = {}; MultiplierMethod[f,g,h,{x1,x2,x3},{1.,1.,1.}] // Timing

TP340   p. 161  Pascal Problem

f = -x1 x2 x3;
h =  {-1.8 + x1 + 2 x2 + 2 x3, x1-1}; g = {};
MultiplierMethod[f,g,h,{x1,x2,x3},{1.,1.,1.},
DualParameter -> True]

Another quadratic problem

criterion =  1/2 u0^2 + 1/6 (u1^2+u2^2+u3^2) ;

Lag = criterion  +         (* The Lagrangian (not used)  *)
p(u0-v0^2 - 2) +           (* u0    2  *)
p1(u0 + u1 -v1^2 - 5) +    (* u0 + u1    5  *)
p2(u0 + u2 -v2^2 - 2) +    (* u0 + u2    2  *)
p3(u0 + u3 -v3^2 - 1);     (* u0 + u3    1  *)

f = 1/2 u0^2 + 1/6 (u1^2+u2^2+u3^2) ;
u = {u0,u1,u2,u3};   g = {};  (* no equality constraint *)
h = {-u0 + 2, - u0 - u1 + 5, - u0 - u2 + 2 , - u0 - u3 + 1};
MultiplierMethod[f,g,h,u,Table[2,{4}]] // Timing