Difference between revisions of "Lotka Experimental Design (AMPL)"

# ----------------------------------------------------------------
# Lotka Experimental Design problem with collocation (implicit Euler)
# (c) Sebastian Sager
# ----------------------------------------------------------------
param nt > 0 integer;
param nu > 0 integer;
param ntperu > 0 integer;
param nx > 0 integer;
param np > 0 integer;
param npar > 0 integer;
param nF > 0 integer;
 
param tf > 0;
param h;
param initialization;
 
set I:= 0..nt;
set U:= 0..nu-1;
param uidx {I};
 
param p{1..np};
param M;
 
param x0{1..nx};
param G0{1..nx,1..npar};
param F0{1..nF};
 
var x{1..nx,0..nt};
var G{1..nx,1..npar,0..nt};
var F{1..nF,0..nt};
var w{1..nx,0..nu-1} >= 0, <= 1;
var u{0..nu-1} >= 0, <= 1;
 
### Calculate trace of 2x2 inverse matrix by hand and add penalization
minimize covariance: (F[1,nt] + F[3,nt]) / ( F[1,nt]*F[3,nt] - F[2,nt]*F[2,nt])
 + p[9] * 0.5 * h * ntperu * 
   (   w[1,0] + w[2,0] + w[1,nu-1] + w[2,nu-1]
     + 2 * sum {i in U diff {0,nu-1} } ( w[1,i] + w[2,i] )
   );
 
### Predator prey dynamics
subject to dynamicX1{i in 1..nt}:
	x[1,i] = x[1,i-1] + h * ( p[1]*x[1,i] - p[2]*x[1,i]*x[2,i] - p[5]*x[1,i]*u[uidx[i-1]] );
 
subject to dynamicX2{i in 1..nt}:
	x[2,i] = x[2,i-1] + h * ( -p[3]*x[2,i] + p[4]*x[1,i]*x[2,i] - p[6]*x[2,i]*u[uidx[i-1]] );
 
### Variational differential equations
### Calculated by hand for par = (p2, p4)
# fx11 = p[1] - p[2]*x[2,i-1] - p[5]*u[uidx[i-1]]
# fx12 = - p[2]*x[1,i-1]
# fx21 = p[4]*x[2,i-1]
# fx22 = -p[3] + p[4]*x[1,i-1] - p[6]*u[uidx[i-1]]
# fp11 = -x[1,i-1]*x[2,i-1]
# fp22 = x[1,i-1]*x[2,i-1]
 
# G11: fx11*G11 + fx12*G21 + fp11
subject to dynamicG11{i in 1..nt}:
	G[1,1,i] = G[1,1,i-1] + h * ( 
      ( p[1] - p[2]*x[2,i] - p[5]*u[uidx[i-1]] ) * G[1,1,i]
    + ( - p[2]*x[1,i] ) * G[2,1,i]
    + ( -x[1,i]*x[2,i] )
    );
 
# G12: fx11*G12 + fx12*G22 
subject to dynamicG12{i in 1..nt}:
	G[1,2,i] = G[1,2,i-1] + h * ( 
      ( p[1] - p[2]*x[2,i] - p[5]*u[uidx[i-1]] ) * G[1,2,i]
    + ( - p[2]*x[1,i] ) * G[2,2,i]
    );
 
# G21: fx21*G11 + fx22*G21
subject to dynamicG21{i in 1..nt}:
	G[2,1,i] = G[2,1,i-1] + h * ( 
      ( p[4]*x[2,i] ) * G[1,1,i]
    + ( -p[3] + p[4]*x[1,i] - p[6]*u[uidx[i-1]] ) * G[2,1,i]
    );
 
# G22: fx21*G12 + fx22*G22 + fp22
subject to dynamicG22{i in 1..nt}:
	G[2,2,i] = G[2,2,i-1] + h * ( 
      ( p[4]*x[2,i] ) * G[1,2,i]
    + ( -p[3] + p[4]*x[1,i] - p[6]*u[uidx[i-1]] ) * G[2,2,i]
    + ( x[1,i]*x[2,i] )
    );
 
### Fisher information matrix 
# F1: w1 * G11*G11 + w2 * G12*G12
subject to dynamicF1{i in 1..nt}:
	F[1,i] = F[1,i-1] + h * ( w[1,uidx[i-1]] * G[1,1,i] * G[1,1,i] + w[2,uidx[i-1]] * G[1,2,i] * G[2,1,i]);
 
# F2: w1 * G11*G21 + w2 * G12*G22
subject to dynamicF2{i in 1..nt}:
	F[2,i] = F[2,i-1] + h * ( w[1,uidx[i-1]] * G[1,1,i] * G[2,1,i] + w[2,uidx[i]] * G[1,2,i] * G[2,2,i]);
 
# F3: w1 * G21*G21 + w2 * G22*G22
subject to dynamicF3{i in 1..nt}:
	F[3,i] = F[3,i-1] + h * ( w[1,uidx[i-1]] * G[2,1,i] * G[2,1,i] + w[2,uidx[i-1]] * G[2,2,i] * G[2,2,i]);
 
### Summing up measurements
subject to wSum{k in 1..nx}:
	sum{i in 0..nt-1} h*w[k,uidx[i]] <= p[6+k];
 
### Initial values
subject to iniX{k in 1..nx}:	x[k,0] = x0[k];
subject to iniG{k in 1..nx, j in 1..npar}:	G[k,j,0] = G0[k,j];	
subject to iniF{k in 1..nF}:	F[k,0] = F0[k];

# ------------------------------------
# Data: Lotka Experimental Design problem
# ------------------------------------
 
# Change values here
param ntperu := 10;
param nu := 50;
param nt := 500;
param tf := 12;
param nx := 2;
param np := 9;
param npar := 2;
param nF := 3;
param initialization := 1e-1;
 
let x0[1] := 0.5;
let x0[2] := 0.7;
let {k in 1..nx, j in 1..npar} G0[k,j] := 0;
let {k in 1..nF} F0[k] := 0;
let F0[2] := 0;
 
let p[1] := 1.0;
let p[2] := 1.0;
let p[3] := 1.0;
let p[4] := 1.0;
let p[5] := 0.4;
let p[6] := 0.2;
let p[7] := 4.0;
let p[8] := 4.0;
let p[9] := 0;
 
# Set indices of controls corresponding to time points
for {i in 0..nu-1} {
	for {j in 0..ntperu-1} {
		let uidx[i*ntperu+j] := i; 
	}
}
let uidx[nt] := nu-1;
 
 
# Initialization
let h := tf / nt;
 
let {i in 0..nu-1} u[i] := 0.3;
let {k in 1..nx, i in 0..nu-1} w[k,i] := 1.0;
let {k in 1..nx, i in 0..nt} x[k,i] := x0[k];
let {k in 1..nx, j in 1..npar, i in 0..nt} G[k,j,i] := initialization; # G0[k,j];
let {k in 1..nF, i in 0..nt} F[k,i] := initialization; # F0[k];
let {k in 2..2, i in 0..nt} F[k,i] := 0;

# ------------------------------------
# Solve Lotka Experimental Design problem
# ------------------------------------
 
model vplotka.mod;
data vplotka.dat;
 
option solver ipopt;
solve;