SebastianSager: corrected Fisher information matrix

2012-10-18T16:18:22Z

corrected Fisher information matrix

← Older revision		Revision as of 16:18, 18 October 2012
Line 92:		Line 92:

	### Fisher information matrix		### Fisher information matrix
	# F1: w1 * G11G11 + w2 ~~G12~~*~~G12~~		# F1: w1 * G11G11 + w2 G21*G21
	subject to dynamicF1{i in 1..nt}:		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[1,2,i]);		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~~G21~~ + w2 ~~G12~~*G22		# F2: w1 * G11G12 + w2 G21*G22
	subject to dynamicF2{i in 1..nt}:		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-1]] * G[1,2,i] * G[2,2,i]);		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 * ~~G21~~~~G21~~ + w2 G22*G22		# F3: w1 * G12G12 + w2 G22*G22
	subject to dynamicF3{i in 1..nt}:		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]);		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		### Summing up measurements

SebastianSager: two typos, thanks to Henry Kar Ming Chan for pointing this out

2012-10-03T07:25:35Z

two typos, thanks to Henry Kar Ming Chan for pointing this out

← Older revision		Revision as of 07:25, 3 October 2012
Line 94:		Line 94:
	# F1: w1 * G11G11 + w2 G12*G12		# F1: w1 * G11G11 + w2 G12*G12
	subject to dynamicF1{i in 1..nt}:		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]);		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[1,2,i]);

	# F2: w1 * G11G21 + w2 G12*G22		# F2: w1 * G11G21 + w2 G12*G22
	subject to dynamicF2{i in 1..nt}:		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]);		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-1]] * G[1,2,i] * G[2,2,i]);

	# F3: w1 * G21G21 + w2 G22*G22		# F3: w1 * G21G21 + w2 G22*G22

SebastianSager: Initial setup of AMPL files

2012-08-31T09:14:35Z

Initial setup of AMPL files

New page

This page contains a discretized version of the MIOCP [[Lotka Experimental Design]] in [http://www.ampl.org AMPL] format. You should be aware of the comments regarding discretization made on the [[:Category:AMPL|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,
<source lang="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];
</source>
a data file vplotka.dat,
<source lang="AMPL">
# ------------------------------------
# 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;
</source>
and a running script vplotka.run,
<source lang="AMPL">
# ------------------------------------
# Solve Lotka Experimental Design problem
# ------------------------------------

model vplotka.mod;
data vplotka.dat;

option solver ipopt;
solve;
</source>

[[Category:AMPL]]

Lotka Experimental Design (AMPL) - Revision history

SebastianSager: corrected Fisher information matrix

SebastianSager: two typos, thanks to Henry Kar Ming Chan for pointing this out

SebastianSager: Initial setup of AMPL files