TACO initial commit

# Solving ODE/DAE optimal control problems with AMPL and MUSCOD-II

# Christian Kirches, Sven Leyffer

# July 2011

option presolve 0;

suffix scale IN >= 0;

suffix interp_to IN;

suffix type symbolic IN;

suffix slope_min IN;

suffix slope_max IN;

option type_table '\

1	u0			piecewise constant control\

2	u1			piecewise linear control\

3	u1c			piecewise linear continuous control\

4	u3			piecewise cubic control\

5	u3c			piecewise cubic continuous control\

6   dae			DAE algebraic state variable\

7	Lagrange	Force a least-squares objective to be treated as a Lagrange type one\

8   dpc         Treat ODE initial value constraint as decoupled point constraint\

';

function diff;

function eval;

function integral;

#

# Batch distillation model

# See Diehl'2006

#

include OptimalControl.mod;

# time and free end-time

var t;

var tf := 2.5, >= 0.5, <= 10.0;

# constant parameters

param Pur := 0.99;		# percent

param V := 100.0;		# mol/h

param m := 0.1;			# mol

param mC := 0.1;		# mol

# control

var R := 8.0, >= 0.0, <= 15.0;

let R.type := "u1";

let R.scale := 0.1;

# differential states

param NDIS := 5;		# PDE discretization points

var M0;

var x{0..NDIS+1};

var MD;

var xD;

var alpha;

# algebraic expressions eliminated by AMPL's presolve

var L = R/(1+R)*V;

var y{i in 0..NDIS} = (1+alpha)*x[i]/(alpha+x[i]);

var dot0 = ( L*x[1] - V*y[0] + (V-L)*x[0] ) / M0;

var dot{i in 1..NDIS} = ( L*x[i+1] - V*y[i] + V*y[i-1] - L*x[i] )/m;

var dotNDISp1 = V/mC * (-x[NDIS+1] + y[NDIS]);

# objective function

minimize Compromise:

	eval (t - MD, tf);

# terminal constraint

subject to Purity_Constraint:

	eval(xD, tf) >= Pur;

# ODE system

subject to ODE_M0:

	diff(M0, t) = -V+L;

subject to ODE_x_0:

	diff(x[0], t) = dot0;

subject to ODE_x{i in 1..NDIS}:

	diff(x[i], t) = dot[i];

subject to ODE_x_NDISp1:

	diff(x[NDIS+1], t) = dotNDISp1;

subject to ODE_MD:

	diff(MD, t) = V-L;

subject to ODE_xD:

	diff(xD, t) = (V-L) * (x[NDIS] - xD)/MD;

subject to ODE_alpha:

	diff(alpha, t) = 0.0;

# Initial value constraints

subject to IVC_M0:

	eval(M0, 0) = 100.0;

subject to IVC_x_0:

	eval(x[0], 0) = 0.5;

subject to IVC_x{i in 1..NDIS+1}:

	eval(x[i], 0) = 1;

subject to IVC_MD:

	eval(MD, 0) = 0.1;

subject to IVC_xD:

	eval(xD, 0) = 1;

subject to IVC_alpha:

	eval(alpha, 0) = 0.2;

option solver muscod;

option muscod_options "solve=solve_tbox deriv=sparse";

option muscod_auxfiles "cr";

problem batchdist: Compromise, 

                   Purity_Constraint, ODE_M0, ODE_x_0, ODE_x, ODE_x_NDISp1, ODE_MD, ODE_xD, ODE_alpha,

				   IVC_M0, IVC_x_0, IVC_x, IVC_MD, IVC_xD, IVC_alpha,

                   R, M0, x, MD, xD, alpha, t, tf;

solve batchdist;

#

# The classical brachistochrone problem

#

# see e.g. Betts'93

#

# Christian Kirches

# July 2011

#include OptimalControl.mod;

option presolve 0;

suffix scale IN >= 0;

suffix interp_to IN;

suffix type symbolic IN;

suffix slope_min IN;

suffix slope_max IN;

option type_table '\

1	u0			piecewise constant control\

2	u1			piecewise linear control\

3	u1c			piecewise linear continuous control\

4	u3			piecewise cubic control\

5	u3c			piecewise cubic continuous control\

6   dae			DAE algebraic state variable\

7	Lagrange	Force a least-squares objective to be treated as a Lagrange type one\

8   dpc         Treat ODE initial value constraint as decoupled point constraint\

';

function diff;

function eval;

function grid;

function integral;

var t;

var tf := 1.0, >= 0.1, <= 1.0;

let tf.scale := 0.5;		# improves convergence

var x := 0, >= 0, <= 1;

var y := 0, >= 0, <= 1;

var v := 0, >= 0, <= 8;

var a := 0.5, >= 0, <= 1.57079327;

let a.type := "u1";

let a.slope_min := -10.0;

let a.slope_max := +10.0;

param gravity := 32.174;

minimize 

EndTime: eval (t,tf);

let EndTime.scale := 0.1;

subject to 

ODE_x: diff(x,t) = v*cos(a);

ODE_y: diff(y,t) = v*sin(a);

ODE_v: diff(v,t) = gravity*sin(a);

IVC_x: eval(x,0) = 0;

IVC_y: eval(y,0) = 0;

IVC_v: eval(v,0) = 0;	

TC_x:  eval(x,tf) = 1.0;

# treating IVCs as boundary constraints improves convergence

let IVC_x.type := "dpc";

let IVC_y.type := "dpc";

let IVC_v.type := "dpc";

option solver muscod;

option muscod_options "nshoot=20 hess=hess_update stiff=0 atol=1e-6 itmax=20";

problem brac: EndTime, ODE_x, ODE_y, ODE_v, IVC_x, IVC_y, IVC_v, TC_x, t, tf, x, y, v, a;

option muscod_auxfiles "cr";

solve;