# MOQP-GOALe.mod
# Original AMPL coding by Sven Leyffer, Argonne National Laboratory
#
# MPEC formulation of multi-objective QP problem.
# Aim is to find a good description of the Pareto set.
#
# Formulation GOAL (from GOAL programming)
#
# A generic MOQP (matrices are upper case)
#
#	minimize   { x^T G_k x + g_k^T x , for k=1,...,p }
#	subj. to   A^T x >= b
#                      x >= 0
#
# for data files, see below.
# ------------------------------------------------------------

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

# ... problem data
param G{P,N,N}  default 0;      # ... objective Hessian
param g{P,N}    default 0;	# ... objective gradients
param A{M,N}    default 0;	# ... constraint Jacobian
param b{M}      default 0;	# ... lower 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};

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{pp in P,k in K} =   sum{i in N, j in N} G[pp,i,j] * x[i,k] * x[j,k] / 2
                        + sum{i in N} g[pp,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
   c{j in M, k in K}: sum{i in N} A[j,i]*x[i,k] >= b[j]
                        complements  y[j,k] >= 0;

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

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

############################################################################
# ... choose one of following data sets
data data/moqp-0001.dat;
data;
  param   : lz,       uz :=
   	   2 -44.4324 127.791
   	   3 -34.226  347.572 ;
##data data/moqp-0002.dat;
##data;
## param   : lz,       uz :=
##   	   2 -80.6163  27.7028 
##   	   3 -21.7589  245.881 ;
##data data/moqp-0003.dat;
##data;
## param   : lz,       uz :=
##   	   2 -65.2195  -13.5559
##   	   3 -4.2092   147.144 ;
##data data/moqp-0004.dat;
##data;
## param   : lz,       uz :=
##   	   2 -31.2347  67.8435
##   	   3 -9.11777  41.1228 ;
##data data/moqp-0005.dat;
##data;
## param   : lz,       uz :=
##   	   2 -67.8913  184.357
##   	   3 -47.1084  383.309 ;
##data data/moqp-0011.dat;
##data;
## param   : lz,       uz :=
##   	   2 -80.8264  44.9205
##   	   3 -70.9535  52.8758
##data data/moqp-0012.dat;
##data;
## param   : lz,       uz :=
##   	   2 -34.7929  80.8069
##   	   3 -35.8924  105.814
##data data/moqp-0013.dat;
##data;
## param   : lz,       uz :=
##   	   2 -41.1619  82.9196
##   	   3 -51.9199  80.572
##data data/moqp-0014.dat;
##data;
## param   : lz,       uz :=
##   	   2 -59.4905  58.579
##   	   3 -40.5348  87.7915
##data data/moqp-0015.dat;
##data;
## param   : lz,       uz :=
##   	   2 -55.307  10401.3
##   	   3 -232.556 88.5839
############################################################################

data;
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;

# ... generate distribution of weights
# ... 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, f, u,                                  # variables
	     c, KKT, GOAL;                                # constraints

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