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

<< MultiplierMethod.m

Needs["Optimization`MultiplierMethod`"];

Version /. Options[MultiplierMethod]

$Version

Date[Round[TimeZone[]]]

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

(* Check for inappropriate convergence goal,
   and DualParameter processing. *)
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]

(* The last attempt found an inadmissible solution
   in pre Mma 5.0 versions;
   try it again with a better starting point. *)

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

(* Or ramp the penalty to reduce infeasibility *)

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

(* 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

(* increase penalty to try and reduce 
   optimality violations present 
   in pre Mma 5.0 versions; 
   notice that we get the same answer *)

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},
	AugmentationGrowth->2,
	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

(* use a better starting point 
   to find a better solution *)

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};
MultiplierMethod[f,g,h,{x1,x2},{10.,10.}] // Timing

(* Try again with a better starting point; 
   better result and faster in pre Mma 5.0 versions. *)

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

(* ramping augmentation takes about same time
in pre Mma 5.0 versions *)

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

(* Try again with reduced convergence goals *)

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

(* Try again with a better starting point; 
   gives better results in pre Mma 5.0 versions. *)

f = -x2;
g =  { x1^2 + x2^2 - 1 };
h = {-x1, -1 - x1 + 2 x2};
MultiplierMethod[f,g,h,{x1,x2},{1.,1.},
	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

(* increase iterations to reduce infeasibility *)

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

(* ramp penalty
   to reduce infeasibility in pre Mma 5.0 versions,
   and run faster *)

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

(* ramp penalty 
   to reduce infeasibility in pre Mma 5.0 versions,
   iterate former answer to speed convergence,
   limit augmentation growth to help satisfy KKT *)

f = (x1-20)^2 + (x2 + 20)^2;
g =  {x1^2 /100. + x2^2 -1}; 
h = {};
MultiplierMethod[f,g,h,{x1,x2},{x1,x2}/.ans[[2,2]],
	Augmentation->10, AugmentationGrowth->1.2
] // 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.}] // Timing

(* increase iterations to reduce infeasibility *)

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

(* ramp penalty to reduce infeasibility 
   in pre 5.0 versions, and run faster *)

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

(* notice bigger penalty 
   is not always better solution 
   (in pre Mma 5.0 versions) *)

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

(* notice bigger penalty yet 
   may yield better solution 
   (in pre Mma 5.0 versions) *)

f = (x1-20)^2 + (x2 + 20)^2;
g =  {x1^2 /100. + 100. x2^2 -1}; 
h = {};
ans = MultiplierMethod[f,g,h,{x1,x2},{0.,0.},
	AugmentationGrowth->4] // 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;
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]

(* as stated, problem is not really well formed;
   with x1 << 0, x2 << 0, and x3 == -x2, the
   minimum is unbounded, with no violations.
   Here, only KKT is violated, indicating we're not
   yet at the optimum.  *)

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]

(* make problem well formed
   by requiring non negative x *)

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

(* still can get into trouble 
   with poor starting point in pre Mma 5.0 versions;
   notice how KKT almost satisfied; complementary
   slack violation indicates non-optimal solution *)

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

Another quadratic problem

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

Problem WRIGHT9

Taken from Yinyu Ye's SOLNP Users' Guide, page 12.  SOLNP is an easy-to-use, fast, powerful, self-contained nonlinear (bound, equality, and/or inequality) constrained optimizer which runs in MATLAB and Octave; see http://dollar.biz.uiowa.edu/col/ye/matlab.html and ftp://dollar.biz.uiowa.edu/pub/yyye/matlab/manual.ps.

x = {x1,x2,x3,x4,x5};
f =	10*x1*x4 - 6*x3*x2^2 + x2*x1^3 +
	9*Sin[x5-x3] + x5^4*x4^2*x2^3  ;
h = {
       x.x - 20,
       -( x1^2*x3 + x4*x5 +2 ),
       -( x2^2*x4 + 10*x1*x5 -5 )
} ;

Use the default settings and the first starting point, and get a reasonable solution.

MultiplierMethod[f, {}, h, x, {1,1,1,1,1},
	DualParameter->True]

(* tweaked to clean up most of the violations
   in pre Mma 4.0 versions *)
MultiplierMethod[f, {}, h, x, {1,1,1,1,1},
	Augmentation -> 10,
	AugmentationGrowth -> 1.1,
	DualParameter->True]

Use other starting point, NB: must boost penalty to avoid infeasible solution.

x0 ={
       1.091,
        -3.174,
        1.214,
        -1.614,
        2.134
        } ;

(* tweaked to find a good value *)
y = MultiplierMethod[f, {}, h, x, x0,
	Augmentation->10*10^3, 
	DualParameter->True
	]

(* alternate tweaked recipe to get good results *)

z = MultiplierMethod[f, {}, h, x, x0, 
	Augmentation->10^3,
	AugmentationGrowth->1.25,
	DualParameter->True,
	CheckConvergence->True
	]

Notice that we can't make the penalty too big: if the penalty is too big, then there is miniscule advancement in pre Mma 5.0, and suboptimal solutions in all cases.  Note how KKT condition violation indicates that there is a problem.

f  /.  (Thread[x  ->  x0])

y = MultiplierMethod[f, {}, h, x, x0, 
	Augmentation->10^8,
	AugmentationGrowth->10,
	DualParameter->True
]

In pre Mma 5.0, we were able to locate another solution we happened upon when starting from an optimal solution.

z2 = { 1.4796 , -2.6366, 1.0547, -1.6115, 2.6739 } ;

{f, h}  /.  Thread[x  -> z2]

(* tweaked to locate other root in pre Mma 5.0*)
zz = MultiplierMethod[f, {}, h, x, z2, 
	Augmentation->100,
	AugmentationGrowth->1.1,
	OuterIterations->60,
	DualParameter->True
	]

Demonstrate that we can get silly results; true solution is unbounded.

Clear[x];
MultiplierMethod[x, {}, {}, {x}, {0}]

The End