reset;

option randseed'';

param I := round(Uniform(35,45));	# number of students
param P := round(Uniform(10,12));	# number of projects

param w {i in 1..I,k in 1..3} := #round(Uniform(0.5,P+0.5));
					 round((Uniform(0.5,(P+0.5)^(1/0.7)))^0.7); # k-th choice of student i
									# not yet distinct
					 
param K := 10;						# minimal number of projects

var z {i in 1..I,p in 1..P} binary; # whether student i gets project p
var x {i in 1..I,k in 1..4} = 		# whether student i gets choice k (k=4: none top 3)
				if k<4 then z[i,w[i,k]] else 1-sum{j in 1..3} z[i,w[i,j]];	
									
var y {p in 1..P} = sum{i in 1..I} z[i,p];	# number of students for project p
var m {p in 1..P} binary;			# for linearization 0,3,4,5,6


maximize LP: sum{i in 1..I} (100*x[i,1]+10*x[i,2]+1*x[i,3]-1000*x[i,4]); # objective

subject to NB1{i in 1..I}: sum{k in 1..4} x[i,k] = 1;		# unique k-th choice

#subject to NB2{i in 1..I,k in 1..3}: x[i,k] = z[i,w[i,k]];	# coupling x and z

#subject to NB3{p in 1..P}: sum{i in 1..I} z[i,p] = y[p];	# coupling z and y

subject to NB4{p in 1..P}: y[p] >= 3*m[p];					# 0 or 3-6 students/project
subject to NB5{p in 1..P}: y[p] <= 6*m[p];					# linearized

subject to NB6{i in 1..I}: sum{p in 1..P} z[i,p] = 1;		# each student gets exactly one project

subject to NB7:            sum{p in 1..P} m[p] >= K;		# at least K projects


option solver cplex;

solve;														# solution

display w;

display x;
display z;
display y;

display I,P,K;