#######################################################################
# ex-4-GS.mod: Gauss-Seidel formulation of small EPEC with two leaders
# 
# Example 4 from Fukushima and Pang, "Quasi-Variational Inequalities,
# Generalized Nash Equlibria, and Multi-Leader-Follower Games", to
# appear in Computational Management Science.
# 
#######################################################################

set I := 1..2;
param a{i in I} := if (i==1) then 1 else -1;

# ... leader variables
var x{I} >= 0, <= 1;      # ... leader 1 & 2

# ... follower variables
var s >= 0;
var y >= 0;

# ... leaders' objective functions
minimize f{i in I}: a[i]*( x[i]/2 + y );

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

#################################################################
###   GAUSS-SEIDEL 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-Seidel
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 leader{i in I}: 
	f[i],  			                # ... objective
	x[i], y, s,			 	# ... variables
	slack, compl;	                        # ... constraints

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

   # ... solve leader 1
   problem leader[1];
   fix x[2] := x0[2];
   solve leader[1]; 
   display _varname, _var.lb, _var     , _var.ub; 
   display _conname, _con.lb, _con.body, _con.ub;
   let x1[1] := x[1];
   let s1[1] := s;
   let y1[1] := y;
  
   # ... solve leader 2
   problem leader[2];
   fix x[1] := x1[1];  ### Gauss-Seidel
   solve leader[2]; 
   display _varname, _var.lb, _var     , _var.ub; 
   display _conname, _con.lb, _con.body, _con.ub;
   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-Seidel)
   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