# outrata4-GJ.mod	QUR-AN-NCP-5-0-4
# Original AMPL coding by Sven Leyffer, Argonne National Laboratory, 2004.
# 
# Gauss-Jacobi method for solving ...
# 
# Multi-Leader-Follower Game (MLF)) derived from outrata31-outrata34, 
# see J. Outrata, SIAM J. Optim. 4(2), pp.340ff, 1994. All Stackelberg 
# players have the same constraints, but different objectives.
# 
#########################################################################

# ... parameters to give number of leaders
param nl integer, default 4;       # number of leaders
set L := 1..nl;                    # set of leaders
set I := 1..4;                     # ... useful index set

# ... Stackelberg leaders' variables
var x{L} >= 0, <= 10;
var xa = sum{l in L} x[l] / card(L);

# ... Followers' variables
var y{I} >= 0;
var s{I} >= 0;

# ... objective functions of leaders
minimize f1: (x[1]-1)^2 + (y[1]-3)^2 + (y[2]-4)^2 + (y[3]-1)^2 + 10*y[4]^2;
minimize f2: 1.5*x[2]^2 + (y[1]-3)^2 + (y[2]-4)^2 + (y[3]-1)^2 + 10*y[4]^2;
minimize f3:     x[3]^2 + (y[1]-3)^2 + (y[2]-4)^2 + (y[3]-1)^2 + 10*y[4]^2;
minimize f4:     x[4]^2 + (y[1]-3)^2 + (y[2]-4)^2 + (y[3]-1)^2 + 10*y[4]^2;

subject to 

   # ... assume followers repond to average of leaders (xa)
   nlcs1: s[1] = (1 + 0.2*xa)*y[1] - (3 + 1.333*xa) - 0.333*y[3] + 2*y[1]*y[4];
   nlcs2: s[2] = (1 + 0.1*xa)*y[2] - xa + y[3] + 2*y[2]*y[4];
   nlcs3: s[3] = 0.333*y[1] - y[2] + 1 - 0.1*xa;
   nlcs4: s[4] = 9 + 0.1*xa - y[1]^2 - y[2]^2;

   compl{j in I}: 0 <= s[j]  complements  y[j] >= 0;

#################################################################
###   GAUSS-JACOBI ITERATION FOR MULTI-LEADER-FOLLOWER GAME   ###
#################################################################
option solver mpec;

# ... save copies of leader's variables
param x0{L} default 0;      # ... iteration k
param s0{I,L} default 0;
param y0{I,L} default 0;
param x1{L};                # ... iteration k+1
param s1{I,L};
param y1{I,L};

# ... stopping tests
param Tol >= 0, default 1E-5;           # ... tolerance for Gauss-Jacobi
param MaxIt >= 0, integer, default 20;  # ... maximum number of iterations
param DiffNorm default 0;

# ... initialize variables
let{i in L} x[i] := x0[i];
let{i in I} s[i] := s0[i,1];
let{i in I} y[i] := y0[i,1];

# ... define single leader MPEC
problem outrata41:
	f1,  			          # ... objective
	x[1], {i in I}y[i], {i in I}s[i], # ... variables
	nlcs1, nlcs2, nlcs3, compl;       # ... constraints
problem outrata42:
	f2,  			          # ... objective
	x[2], {i in I}y[i], {i in I}s[i], # ... variables
	nlcs1, nlcs2, nlcs3, compl;       # ... constraints
problem outrata43:
	f3,  			          # ... objective
	x[3], {i in I}y[i], {i in I}s[i], # ... variables
	nlcs1, nlcs2, nlcs3, compl;       # ... constraints
problem outrata44:
	f4,  			          # ... objective
	x[2], {i in I}y[i], {i in I}s[i], # ... variables
	nlcs1, nlcs2, nlcs3, compl;       # ... constraints

# ... Gauss-Jacobi
for {k in 1..MaxIt}{

   # ... solve leader 1
   problem outrata41;
   fix x[2] := x0[2];
   fix x[3] := x0[3];
   fix x[4] := x0[4];
   solve outrata41;
   let x1[1] := x[1];
   let{i in I} s1[i,1] := s[i];
   let{i in I} y1[i,1] := y[i];

   # ... solve leader 2
   problem outrata42;
   fix x[1] := x0[1];
   fix x[3] := x0[3];
   fix x[4] := x0[4];
   solve outrata42;
   let x1[2] := x[2];
   let{i in I} s1[i,2] := s[i];
   let{i in I} y1[i,2] := y[i];

   # ... solve leader 3
   problem outrata43;
   fix x[1] := x0[1];
   fix x[2] := x0[2];
   fix x[4] := x0[4];
   solve outrata43;
   let x1[3] := x[3];
   let{i in I} s1[i,3] := s[i];
   let{i in I} y1[i,3] := y[i];

   # ... solve leader 4
   problem outrata44;
   fix x[1] := x0[1];
   fix x[2] := x0[2];
   fix x[3] := x0[3];
   solve outrata44;
   let x1[4] := x[4];
   let{i in I} s1[i,4] := s[i];
   let{i in I} y1[i,4] := y[i];

   # ... check for termination
   let DiffNorm := sqrt( sum{l in L}( (x0[l] - x1[l])^2 )
                         + sum{i in I,l in L}( (y0[i,l] - y1[i,l])^2 
                                               + (s0[i,l] - s1[i,l])^2 ));
   display DiffNorm, Tol, k;
   if (DiffNorm <= Tol) then {
      break;
   }; # end if

   display x, x0, x1;

   # ... take the step (Gauss-Jacobi)
   let{l in L} x0[l] := x1[l];
   let{l in L, i in I} s0[i,l] := s1[i,l];
   let{l in L, i in I} y0[i,l] := y1[i,l];

}; # end for