#######################################################################
# ex-001-GS.mod: Gauss-Jacobi formulation of small EPEC with two leaders
# 
# leader i's problem is:
#
#     minimize    (x_i+1)^2
#     subject to  s = x_1 + x_2 + y
#                 0 <= s  _|_  y >= 0
#######################################################################

set I := 1..2;

# ... leader variables
var x{I};

# ... follower variables
var s;
var y;

# ... leader objective
minimize leader{i in I}: (x[i]+1)^2;

subject to 

	slack: s = x[1] + x[2] + y;
	compl: 0 <= s   complements   y >= 0;

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

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

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

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

# ... define single leader MPEC
problem Stackelberg{i in I}: 
	leader[i],			        # ... objective
	x[i], y, s,			 	# ... variables
	slack, compl;	                        # ... constraints

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

   # ... solve leader 1
   problem Stackelberg[1];
   fix x[2] := x0[2];
   solve Stackelberg[1];
   let x1[1] := x[1];
   let s1[1] := s;
   let y1[1] := y;
  
   # ... solve leader 2
   problem Stackelberg[2];
   fix x[1] := x0[1];
   solve Stackelberg[2];
   let x1[2] := x[2];
   let s1[2] := s;
   let y1[2] := y;

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

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

}; # end for