# MOLPg-SUM1e.mod
# Original AMPL coding by Sven Leyffer, Argonne National Laboratory
#
# MPEC formulation of multi-objective MOLPg problem.
# Aim is to find a good description of the Pareto set.
#
# Formulation SUM (find convex sum representation)
# one weght == 1
#
# A generic MOLP (matrices are upper case; note transpose!)
#
#	minimize   C^T x
#	subj. to   A^T x <= b
#
# for data files, see ...
# 
# filename		reference
# ------------------------------------------------------------
# MOLPg-001.dat		[Steu:86] p. 411, Exercise 13-4C 
# MOLPg-002.dat		[Steu:86] p. 415, Exercise 13-14C
# MOLPg-003.dat		[Steu:86] p. 416, Exercise 13-15C
# 			... there are many more MOLPs in there!
# ------------------------------------------------------------
# References:
# ~~~~~~~~~~
# [Steu:86]	R.E. Steuer, Multiple Criteria Optimization: 
#		Theory, Computation and Applications, John 
#		Wiley, 1986, New York.
# ------------------------------------------------------------

param ModelNo >= 1, <= 3, integer, default 1;

# ... problem dimensions & index sets
param n >= 2, integer;	# ... number of variables
param m >= 2, integer;	# ... number of constraints
param p >= 2, integer;	# ... number of objectives
set N  := 1..n;	
set M  := 1..m;
set P  := 1..p;
set P0 := 2..p;

# ... problem data
param C{P,N};				# ... objective gradients
param A{M,N};				# ... constraint Jacobian
param b{M};				# ... upper bounds on A^T x

# ... multi-obnjective stuff
param nK integer, >= 0, default 10;     # ... number of efficient points
set K := 1..nK;	                	# ... index to efficient points
set KxK within K cross K;      		# ... pairs for cross distances
param lz{P0};
param uz{P0};

# ... variables
var x{N,K};				# ... variables
var y{M,K};                     	# ... multiliers
var z{j in P0,K} >= lz[j], <= uz[j];	# ... multi-objective GOALs
var u{P0,K};                            # ... multipliers of GOAL constraints
var eta;

# ... objective functions
var f{j in P,k in K} = - sum{i in N} C[j,i]*x[i,k];

maximize min_dist: eta;

subject to

   # ... definition of eta (all cross distances in objective space)
   ubd_eta { (k1,k2) in KxK}: eta <= sum{j in P} (f[j,k1] - f[j,k2])^2;

   # ... linear constraints & complementary slackness
   lin{j in M, k in K}: sum{i in N} A[j,i]*x[i,k] <= b[j]
                        complements  y[j,k] >= 0;

   # ... KKT conditions
   KKT{i in N, k in K}: 0 <= - C[1,i] - sum{j in P0} C[j,i]*u[j,k]
                                      + sum{j in M}  A[j,i]*y[j,k] 
                        complements x[i,k] >= 0;

   # ... GOAL constraints
   GOAL {j in P0, k in K}: 0 <= u[j,k]  complements  f[j,k] <= z[j,k];

   # ... GOAL constraints
   GOALz{j in P0, k in K}: f[j,k]  = z[j,k];
   GOALu{j in P0, k in K}: u[j,k] >= 0;

############################################################################
data; let ModelNo := 3;
if ModelNo == 1 then data MOLPg-001.dat;
if ModelNo == 2 then data MOLPg-002.dat;
if ModelNo == 3 then data MOLPg-003.dat;
data;
# ... uncomment the appropriate pay-off matrix
# pay-off matrix for ModelNo == 1
#  f1        f2         f3
# -26.0167   3.87073    5.82437
# -0.5       -12        0.583333
# -11.3305   2.78814    -9.11864
##  param   : lz,       uz :=
##   	   2 -12       3.87073
##   	   3 -9.11864  5.82437 ;
# pay-off matrix for ModelNo == 2
#  f1        f2         f3
# -13.375    -0.208333  -1.39583
# -4.375     -14.25     -6.5
# -4.83333   -5.27778   -12.7778
##   param   : lz,       uz :=
##	   2 -14.25    -0.208333
##	   3 -12.7778  -1.39583  ;
# pay-off matrix for ModelNo == 3
#  f1        f2         f3
# -13.4      -3.2       2.9
# -4.22883   -10.6108   1.72063
# 2.375      0          -13.375
   param   : lz,       uz :=
   	   2 -10.6108       0    
   	   3 -13.375        2.9 ;
############################################################################

let nK := 10;
# ... definition of cross distance indices
let KxK := { };
for {k1 in {1..nK}}{
   let {k2 in {1..nK}: k1 <> k2} KxK := KxK union { (k1,k2) };
};
display KxK;

# ... distribute the payoff levels
fix {k in K} z[2,k] := lz[2] + k/nK*(uz[2]-lz[2]) + (nK-k)*Uniform01();
fix {k in K} z[3,k] := lz[3] + k/nK*(uz[3]-lz[3]) + (nK-k)*Uniform01();
display z;

# ... NCP model
problem NCP: x, y, u,                                   # variables
	     lin, KKT, GOAL;                            # constraints

unfix z;

# ... MPCC model
problem MPCC: min_dist,                                 # objective
              x, y, u, f, eta, z,                       # variables
	      lin, KKT, GOALz, GOALu, ubd_eta;                  # constraints