Lotka Experimental Design (AMPL)

This page contains a discretized version of the MIOCP Lotka Experimental Design in AMPL format. You should be aware of the comments regarding discretization made on the AMPL overview page.

AMPL

The model in AMPL code for a fixed control discretization grid with a collocation method, using an implicit Euler. We need a model file vplotka.mod,

# ----------------------------------------------------------------
# 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 * G21*G21
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[2,1,i] * G[2,1,i]);
 
# F2: w1 * G11*G12 + w2 * G21*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[1,2,i] + w[2,uidx[i-1]] * G[2,1,i] * G[2,2,i]);
 
# F3: w1 * G12*G12 + w2 * G22*G22
subject to dynamicF3{i in 1..nt}:
	F[3,i] = F[3,i-1] + h * ( w[1,uidx[i-1]] * G[1,2,i] * G[1,2,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];

a data file vplotka.dat,

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

and a running script vplotka.run,

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