Dynamic modelling

From ASCEND
Jump to navigation Jump to search
NLA
QRSlv
CMSlv
IPSlv
NLP
CONOPT
IPOPT
TRON
MINOS
Opt
NGSlv
DAE/ODE
IDA
LSODE
DOPRI5
RADAU5
LA
Linsolqr
Linsol
LP
MakeMPS
Logic
LRSlv

Dynamic modelling in ASCEND is where you solve a system of ordinary differential equations (ODEs) or differential/algebraic equations for the time-varying values of variables in the model over a certain interval, given the initial conditions of the model. We also refer to this as integration. Dynamic modelling can be done within the ASCEND language providing the relationship between variables and their derivatives can be established.

About ODE/DAE solvers

ODE solvers of the LSODE and Runge-Kutta class solve problems that can be written:

<math>\dot{\mathbf{y}} = \mathbf{F} \left(\mathbf{y},t \right)</math>

by projection schemes for <math>\mathbf{y}\left(t+dt\right)</math> which require solving some set of equations based on Jacobian matrix of partial derivatives <math>\mathbf{J} = \partial \mathbf{F} / \partial \mathbf{y}</math> (also referred to as <math>\partial \mathbf{\dot{y}}/\partial \mathbf{y}</math>) being nonsingular. Variants of ODE solvers may allow convergence with limited singularity in <math>\mathbf{J}</math> by using singular value decomposition (SVD) and assumptions about the properties of <math>\mathbf{F}</math>, but this is not general or robust enough to be used in arbitrary models created with ASCEND.

The ASCEND representation of differential equations allows the more general

<math>\mathbf{F}\left(\mathbf{y},\mathbf{\dot{y}},t,\mathbf{z}\right)=\mathbf{0}</math>

where <math>z</math> are apparently non-differential variables (variables that appear in the equations without a corresponding derivative). Considering all variables <math>\mathbf{y}</math>, <math>\mathbf{\dot{y}}</math>, <math>t</math>, <math>\mathbf{z}</math> as algebraic, then the matrix <math>\mathbf{J}</math> can be calculated analytically (within some tolerance) if the following is true:

With y fixed, t fixed, and as many z fixed as necessary, a well posed (square, nonsingular) nonlinear problem is obtainable.

This form is precisely the matrix needed at the initial condition, normally. So for each "function evaluation" of LSODE, we solve the nonlinear problem

<math>\mathbf{F} \left( \mathbf{\dot{y}}, \mathbf{z}_\mathrm{free}\right)=\mathbf{0}</math> such that st. <math>\mathbf{y}</math>, <math>t</math>, <math>\mathbf{y}</math>, <math>\mathbf{z}_\mathrm{fixed}</math> are all fixed

and for each <math>\mathbf{J}</math> needed we take the Jacobian of <math>\mathbf{F}</math> and solve to get a column of sensitivity data. (In principle this can be accelerated by the method of Curtis, Powell, and Reid(citation needed), but that so rarely applies we didn't bother).

With the above technique, it may turn out that <math>\mathbf{J}</math> is singular, even under SVD. This is the usual first indicator of a numerical index problem.

It turns out that in many cases the index problem is even visible structurally in QRSlv -- meaning there is no feasible way to form <math>\mathbf{F}\left(\mathbf{\dot{y}}, \mathbf{z}_\mathrm{free}\right)</math> by computing all the variables <math>\mathrm{\dot{y}}</math> expected and fixing all the variables <math>\mathbf{y}</math> expected. In this case either the problem is badly posed or a proper ODE problem can be recovered by moving some pair <math>\dot{y}_i, y_i</math> into the set <math>\mathbf{z}</math> and perhaps adding extra constraints obtainable by differentiation of a subset of <math>\mathbf{F}</math>.

DASSL and IDA work around low index problems by solving all model and stepping simultaneously, but for problems typical in chemical engineering (setting pressure specifications on multistage operations like distillation) DASSL's technique is not enough for reliable solution even if SVD is used under the covers.

So basically, it's great to use DASSL for performance reasons over LSODE and the like. If the reason for using DASSL, IDA and the like is because LSODE fails, then the modeling approach should be carefully reevaluated.

If the reason for using DASSL/IDA is that they are better at finding initial conditions than LSODE (which has no such functionality), then again the main strength of ASCEND's primary solvers (getting solutions to hard nonlinear problems) is being overlooked. Ben Allan has heard many times of people using dassl as simply a nonlinear solver and integrating nothing at all.

Finally, now that memory is essentially free, boundary value problem (BVP) solution techniques should be of much greater interest than historically, and it just happens that ASCEND has a very nice BVP library: models/bvp.a4l

Syntax for dynamic modelling

The current method for specifying derivatives is by including the ivpsystem.a4l library. Full details are given in Design and Use of Dynamic Modeling in ASCEND IV (see our publications page).

REQUIRE "ivpsystem.a4l";

at the start of your model, then creating a special method

METHOD ode_init;
    dx_dt.ode_id := 1;
    x.ode_id := 1; (* ie these variables are related with ID '1' *)

    dx_dt.ode_type := 2; (* ie 'this is a derivative variable' *)
    x.ode_type := 1; (* ie 'this is a differential variable *)

    t.ode_type:= -1; (* 'tag' the independent variable *)
END ode_init;

which must be run before sending your model to either IDA or LSODE (or a future solver...) for integration.

New syntax

Dante Stroe recently added support for LINK semantics that allows this fairly painful syntax to be improved, however the new code is not yet incorporated into the ASCEND 'trunk' (see svn info). Based on the LINK structure the following syntax can now be used to define derivative chains, for example:

INDEPENDENT t;
DER(dx_dt,x)

This syntax can be used in either the declarative part of the model, either in the Methods section.

  • The INDEPENDENT keyword adds an independent variable to the solver. While it is possible to store multiple independent variables in the solver data structures, the user is currently limited to only one independent variable, since there is no implementation of partial derivatives. (Implementation-wise, a LINK entry with the 'indepenendent' key and the independent variable itself is created)
  • The DER keyword is used to declare the derivative chains in terms of the declared independent variable. The parameters of DER are the differential variables in decreasing order (ode_type in the previous syntax). The derivative chains declared in DER can be arbitrarily long, however the solver does not handle more than 2nd order derivatives. (Implementation-wise, a LINK entry with the 'ode' key and the differential variable is created)

The following models using the aforementioned syntax have been added:

See also