Improved DAE syntax

From ASCEND
Jump to: navigation, search
This article is about planned development or proposed functionality. Comments welcome.

New syntax is proposed for ASCEND for dynamic modelling involving DAEs, as currently solved by LSODE, IDA, and DOPRI5. See also discussion on event handling and integration of conditional models.

Proposal 1: "der(x)" and "INDEPENDENT t"

Proposal is currently

  • add "der(x)" quasi-operator to the language, to allow use to refer to the derivative of a variable with respect to the independent variable
  • "INDEPENDENT t" declarative statement, to allow a variable in the model to be identified as the variable with respect to which integration is to take place.
  • modify current expression parser to record above relationships using LINK semantics implemented by Dante Stroe.
  • modify current DAE solvers to use a shared API to deal with the above new way of representing DAE systems (replacing the current "ode_id" property method).

It is proposed that the "der(x)" operator would be dealt with in the parser by silently replacing "der(x)" with another variable, perhaps named "x_der" which would be automatically declared, LINKed to "x" and then inserted in the appropriate place in the "struct Expr" tree.

This would be a convenience layer on top of the currently-in-preparation LINK functionality. One can always revert back to that implementation for some or all of the variables in a model, if that greater flexibility is required. But the goal is the make the "der(x)" and "INDEPENDENT t" syntax work for 95% of user models.

It is proposed that the TypeDefinition system would be extended to include a mechanism by which types to be used for automatically-generated derivative instances could be declared. For example, type "velocity" would be defined as "derivative_of :== 'velocity'. This is yet to be discussed.

Discussion is occurring on this topic currently (Jul 2009) in the ascend-developers mailing list.

Under this proposal,

MODEL my;
   t IS_A time;

   INDEPENDENT t;
   y IS_A distance;
   der(y) = -y;
END my;

Pseudocode for what might happen in the parser:

if parseexpr == INDEPENDENT && name:
    var = find_var_by_name(name)
    if simulation->indepinstance:
        fail("multiple indep vars specified in simulation")
    simulation->indep_instance = var;

if parseexpr == DER_TOK && LPAREN && name && RPAREN:
    var = find_var_by_name(name)
    varname = get_der_name(name)
    derivvar = find_var_by_name(varname)
    if derivvar is not None:
         if not var_is_linked_as_deriv(derivvar, var):
             fail("incorrectly declared derivative",varname)
         return derivvar
    else:
         if context_is_procedural():
             fail("derivative referred to in procedure was not present in declarative section of model")
         type = get_type_of_var(var)
         dertype = get_deriv_type(type)

         if not dertype:
             fail("no type defined for derivatives of type",type)
         newvar = create_var_of_type(varname, model, dertype)

         link_var_as_derivative(newvar, var)
         return newvar

Higher derivatives

  • how do we deal with second and higher derivatives?

It is proposed that second and higher derivatives would not be permitted in an expression (i.e. no nesting of "der" operators). This would force the user to declare "xdot IS_A solver_var; xdot = der(x);", which actually matches well with the interface provided by DAE solvers like IDA, which don't permit second derivatives, anyway.

Name collisions

  • how to deal with name collisions, eg if the user creates their own "x_der" and then later uses "der(x)" in an expression?

We can check that if "x_der" is NOT linked to "x" as derivative. If it's not, then we can assume that the user has defined their own "x_der" variable, and we can tell them to rename (parser fail -- syntax error).

Partial derivative models

  • according to Ben, there are a range of problems for which it would be desirable to be able to switch the independent variable from, for example, time to a spatial direction, and vice versa.

For problems like this to work, we would need to have syntax that allows the with-regard-to variable being specified, for example "der(x,t)" or "der(s,z)". When integrating with regard to "t", "der(x,t)" would then have to be treated as a derivative, but "der(s,z)" would have to be treated as just a normal variable.

But can this approach coexist with the proposed der(x) plus INDEPENDENT t idea?

Dimension checking for der(x)?

  • how can we ensure that dimension checking works correctly for derivative variables?

If we propose that derivative variables are automatically created, how do we automatically declare their dimensionality? We can work out the dimensionality using dimension arithmetic, but how do we store the result? Normally, dimensions are associated with declared TypeDefinitions, so we might find that we need to create type definitions on the fly, just so that we can enforce dimension checking. Is that plausible? Or should we just throw an error if no suitablely-matching TypeDefinition is found?

INDEPENDENT declaration

  • what if the user doesn't define INDEPENDENT? (we can't do it automatically and still have units of measurement behaving properly)
  • what if the user declares more than one INDEPENDENT? (would we allow this, and just insist that all are consistent?)

A possible solution here is to add a new instance child as part of the simulation instance. The INDEPENDENT statement will result in storing an instance pointer in that location. If multiple INDEPENDENT statements are found, we would require that the instances named are equivalent (in a clique, etc).

Ben pointed out that we could make use of the UNIVERSAL construct to make the time variable visible from everywhere, without having to pass it as a parameter.

As far as integrators are concerned, all we need to implement is "integrator_find_indep_var", which has to store the indentity of the independent variable into our 'IntegratorSystem' object. Any approach in the instance tree that results in a well-defined independent variable as a member of the solver's var list, is fine for the Integrator.

For a first cut, it is prosed that only one INDEPENDENT statement would be allowed in a simulation. Any subsequent INDEPENDENT statement would throw an error and be ignored. For a second cut, multiple INDEPENDENT statements would be allowed. These would be stored as a list in the SimulationInstance. At the time of instantiating the integrator, these would be checked for equivalence -- the variables should all end up being be members of a single clique.

We would require that INDEPENDENT statements came before the first der(x) statement, such that when parsing der(x), the correct link could be created. An error would be like "Unable to parse "der(x)"... not INDEPENDENT variable has yet been assigned.

Does INDEPENDENT have to be declarative? Because otherwise the dimensionality of der(x) instances can't be defined?

Ben's benchmark

Ben says... "i will insist as one of the test cases (if any of this gets implemented) that the index-3 pendulum ode be modelable and solvable by all published methods without gymnastics in the resulting language constructs."

  • does this require that we implement DAE index reduction?
  • does any of our current integrators solve this problem?

Other issues

  • how to handle correctly deal with "der(x)" in methods. What if "x_der" isn't already in the model in that particular case?
  • can we deal with this on-line expr substitution successfully with yacc/bison?
  • do we care about being able to switch the independent variable in a model? is this a potential use case that we need to deal with? (it's not something that other modelling languages seem to accommodate)
  • what are the implications of not being able to explicitly refer to the derivative of a variable except via the "der(x)" syntax?
  • what if a user short-cuts the naming system, and refers directly to "x_der" instead of "der(x)"?

Proposal 2: "der(x,t)"

Proposal is:

  • Provide syntax for specifying derivatives as der(x,y), where both the base variable as well as the independent variable as specified.
  • in general, allow the user to specify derivatives w.r.t. any variable they choose.
  • would still need a way of telling the solver which independent variable to use, eg INDEPENDENT t.
  • When integrating, all derivatives that are not w.r.t the assigned indep var would be treated as normal algebraic variables.
  • It would be possible to switch independent variables at will, eg to solve a BVP for the initial conditions, then march time for IVP solution. (not sure how this work really work...)

This is the approach Ben would like to see, because it facilitates possible change of independent variable.

This approach still requires the independent variable to be specified. It would require an INDEPENDENT statement to be present somewhere.

In this approach, we would be able to generate the necessary LINK statement very easily, but we would still have challenges with dimensional correctness and unit-checking.

Given that this approach still requires the INDEPENDENT statement, I (JP) don't see a problem with the "der(x)" shorthand coexisting with this. For the case where the INDEPENDENT statment has not yet been seen, you could throw an error. Otherwise, you would use the current accepted independent variable. In the case where this was then illegally redefined with another, inconsistent INDEPENDENT statement, you would just fail the parsing process.

Along with event handling

An argument against der(x,t) is that the event handling code has implied 'time' variable; its operation depends on there being a well defined independent variable. Modelica-like operators like "reinit" and "pre" don't have to state the 't' variable all the time, so perhaps it doesn't make sense to always have to declare it in 'der', either?

Support for PDE modelling

If we would like ASCEND to be extensible to modelling partial differential equations, what would be required?

  • some way of defining the boundaries of the domains on which PDEs are defined
  • some way of defining node points within the domains
  • way of defining what variables are defined for each node point
  • automation of calculation of finite difference equations
  • some way of presenting the results
  • significant changes to the ASCEND parser, because array limits would not be defined until the mesh had been read in?

It seems that this would be a major undertaking, and would take ASCEND into a fairly large new area.

Proposal 3: Ben's special properties

Attributes with the prefix '$' would be reserved for specially supported attributes of atoms. Currently all "supported" attibutes are only supported by *solvers* and not the compiler itself. The compiler allows only scalar attributes. This solution gets rid of ode_type in 99% of cases.

The extension looks like parameterized atomic types:

MODEL foo;
  X IS_A something(t);
  t IS_A domain of time over (min, max);
  U.$d(t) = -2*U;
  x[1..3] IS_A domain OF distance WHERE
    x[1] SPANS (min,max),
    x[2] SPANS (min,max),
    x[3] SPANS (min,max);
END foo;

You can ARE_THE_SAME things as much as you like, except for X, t, because they have a detectably special relationship:

ATOM something (t willbe real);
  (* the following will exist lazily (no instantiation until referred. *)
  $d[i] IS_A DERIVATIVE OVER (i) IN t;
  $fixed IS_A boolean; (* time to promote fixed to the compiler build ... or not*)
END something;

And extension:

ATOM pde3_something(
    indexset isa set of integer_constant;
    direction[indexset] willbe real, time willbe real
) REFINES something(time);
    $pd[i,j,k] IS_A DERIVATIVE OVER (i,j,k) IN direction;
END pde_something;

and the extension to relations (for those who like modeling *algorithms* in ascend directly) is

RELATION
  $pd IS_A PARTIAL_DERIVATIVE OVER VARS;
  $pd2 IS_A PARTIAL_DERIVATIVE OVER (VARS,VARS);
END RELATION;

Please note, the types on the RHS of IS_A inside atom definitions are never models and are not definable within the model syntax. There is a yet-to-be-defined table (extendable-- but only at build time of the compiler) of extraordinary beasts that are allowed in atoms. $ is what kicks us into the table evaluation mode.

Please note that domains are essentially (in this scheme) other parameterized atoms -- they end up (in the case of x given above) with

$range_low[indexset];
$range_high[indexset];

and potentially all sorts of additional special parts to maintain relationships with external discretizations (or internal to ASCEND ones).

So all that external pde-stuff john cataloged is easily *named* in ascend notation. As for creating mesh generators, paraview, openfoam, and all the rest real PDE users need: no, that's for the specialists to create, and for us to reuse when building system models with lots of embedded PDEs.