# c-reload-q.mod   OOR2-AN-342-284
# AMPL coding: S. Leyffer, University of Dundee, April 2000.
#
# Randomly generated core reload optimization problem for 
# quadrant of 24 x 24 reactor core (total 104 fuel rods)
#
# Data files: c-reload-104-1.dat

set I;				# ... set of nodes
set L ordered;			# ... set of batches
set M;				# ... set of bundles in batches
set T ordered;			# ... set of time steps
set ND := 1..2;			# ... number of diagonal pairs
set EoC := { last(T) };		# ... end of cycle definition
set Gnz within I cross I;	# ... non-zeros of G

# ... scalar parameters
param conspw;			# ... conspw = Pc in the article
param kfresh;
param alfa;
param flim;
param cytime;			# ... Delta t = cytime/(card(T)-1) = 70 here

# ... array parameters
param G{Gnz};			# ... matrix G[i,j] sparse

# ... continuous variables
var r {I,T}    >= 0, := kfresh/14;	# ... phi in paper	
var kinf {I,T} >= 0, := kfresh;
var keff {T}   >= 0;
var phi  {I,T} >= 0;

# ... x[i,l,m] = 1 if node i contains bundle age l from trajectory m
#var x {I,L,M} binary, := 0.5;		
var x {I,L,M} binary;		


### OBJECTIVE FUNCTION ###
maximize goal: sum{e in EoC} keff[e];

subject to

  # ... every bundle is assigned in 1 node
  sumi{l in L, m in M}: sum{i in I}  x[i,l,m] = 1;

  # ... every node contains 1 bundle
  sumlm{i in I}: sum{l in L, m in M} x[i,l,m] = 1;

  # ... standard tri-linear reloading equation
  plac{i in I, e in EoC}: 
    kinf[i,1] = sum{l in L: l > 1}
                ( sum{m in M} x[i,l,m]*( sum{j in I} x[j,l-1,m]*kinf[j,e] ) )
                + ( sum{m in M} x[i,1,m] ) * kfresh;

  # ... standard kernel equation definition
  kern{i in I, t in T}: 
       sum{(i,j) in Gnz} (G[i,j]*kinf[j,t]*r[j,t]) = r[i,t]*keff[t];

  # ... burnup described as decay of kinf
  kinff{i in I, t in T diff EoC}: 
       kinf[i,t+1] = kinf[i,t]
                     - alfa*conspw*cytime/(card(T)-1) * kinf[i,t]*r[i,t];

  # ... normalize the power density to 1
  cpow{t in T}: sum{i in I} kinf[i,t]*r[i,t] = 1;

  # ... power peaking constraint
  peak{i in I, t in T}: kinf[i,t]*r[i,t] <= flim / card(I);