# Initial-value modelling

This task is to extend ASCEND's ability to solve initial value differential-algebraic problems. Such problems comprise a mixture of ordinary differential and algebraic equations. The "twist" to this extension is to convert the integration equations for stepping the ODEs into ASCEND models and, initially, to use methods and a Tcl script to do the stepping. In this manner, we do not need to add any new solver code to ASCEND.

Attached (Media:DifferentialAlgebraicEqn.pdf) HOWTO chapter for solving DAE initial value models) is a .pdf file documenting the effort to date. It is in the form of a HOWTO in the ASCEND documentation. The code (in its current state) will be added here soon.

## Current System

The following describes the current system - which is NOT that described in the above .pdf file. It corresponds to the solving of the following DAE models and scripts: dyn_flash.a4l, dyn_column.a4l, dyn_separation_demos.a4s, dyn_tank.a4c and a4s.

To understand the current system, the following 10,000 ft view of the ODE solving process and its history may be useful. The IVP integrators for strict ODEs exist at a level completely separate from the details of solving nonlinear square algebraic systems. All an integrator needs in most cases is ydot=f(y) (FEX) and dydot/dy (JAC). JAC can be produced by a simple trick of linear algebra (this is seen in tcltk/interface/Sensitivity.c). The solution of the time-marching equations (find new y given ydots) is handled in LSODE; the computation of ydot given a fixed state-variable set y is done using the conventional nonlinear solver. If the nonlinear solver has a structural singularity with the states y fixed, either the model is ill-posed or the DAE index of the related DAE problem is greater than 1.

The required math is easy. The required analysis of the model instance DAG to produce flat representations takes place in two stages.

• stage 1: fix the state variables y to form a well-posed problem, free the derivative ydot, and solve the nonlinear algebraic problem for a consistent initial condition. This is the role of the slv_system_t in the diagram.
• stage 2: collect the y and ydot pairs from the instance tree in a format useful for ODE solvers and begin the time-marching iteration. The integrator is allowed to change the values of the y variables in the instance tree as time is moved forward. This is the role of the Integ_system. The NLA solver is used to compute the needed ydot values at each t and test y. The linear sensitivity analysis is used to compute JAC as needed. JAC may be dense or sparse in reality, but for ascend purposes it is always stored as dense though it may be sparse in most PDE/chem-mech-eng problems. This is a potential subject of optimization, but a better use of limited resources would be in making the DAE solving approach that Art is suggesting above more user-friendly.

The user interface (how hierarchical, reusable dynamic models get connected to IVP solvers) is described in <a href="https://pse.cheme.cmu.edu/ascend/ftp/pdfHowto/howto-ivp.pdf" class="external text" title="https://pse.cheme.cmu.edu/ascend/ftp/pdfHowto/howto-ivp.pdf" rel="nofollow">howto-ivp.pdf</a>. The approach taken to pairing y with ydot variables is by assigning a matching integer index value to each using METHODs. This was a "no language changes required" prototype. It can be used to handle DAEs or ODEs. It is slightly cumbersome (and definitely "busy work") to require the user-written state identification and pairing methods; additional syntax has been desired for a long time.

-- Ben Allan - 01 Jun 2006

Syntax extensions require in all but the most trivial cases substantial community negotiation and an eye for where ASCEND is going and has been. The big picture for differential syntax is that we want:

• to be able to handle the declaration of PDEs as well as ODEs (eventually).
• to avoid making a special case of time in the language (though in specific models it's fine).
• to make sure the syntax is compatible with WHEN statements being used to model dynamic discontinuities of all sorts.
• to have an integrator version that handles root-finding and DOF changes so discontinuities are converged properly.

It should be noted that boundary descriptions are already part of the language per the conditional modeling work of Rico-Ramirez.

Some of the basic concepts that have been already devised, but not implemented are:

• Declaring an IVP independent variable of any name and dimensionality, and doing it more than once in the same model as (e.g. time) one variable is not necessarily the only variable that might be useful in ivp modeling as the independent variable.
• Declaring a differential state of any solver_var type.
• Declaring a derivative wrt a state and an independent variable using regular atom types in order to get bounds and scaling correct (yes, derivatives do need bounds and scaling in most physical systems). A tricky case that must be handled in a principled way is the definition of an array of derivatives given an array of states, bearing in mind that ascend sets are unordered sets. In the compiler, sets are really documented as, "arbitrarily and consistently ordered", such that if two set instances contain the same elements, any FOR operation which is executed over one set will process the matching elements in the same order for the other set.
• Defining the domain of an independent variable for IVP or PDE. This likely is different from the bounds on the same variable, because the domain affects the dependent variables in a PDE.
• Leaving a declared state OUT of the set of states to be solved for with the IVP solver.
• Leaving a derivative OUT of the set of variables to be solved for with the NLA solver.
• Defining an "espilon after" (t+) for making statements about discontinuities; the implications (is it an entire shadow model lurking?) are tricky.

-- Ben Allan - 02 Jun 2006

Ben has outlined some very important thoughts we need to discuss in his list above.

I believe there are (at least) two major issues for creating a new DAE solver in ASCEND - if that were our only short term goal. One involves getting the math right. This I hope I have done in writing the .pdf document above and in writing and testing the approach in the models in the folder trunk/models/westerberg/ivpNondimensional in the subversion repository. (Note that these models do not install when creating and installing ASCEND so look for them only in a code checkout.) The second issue is to add to the language the ability to identify the independent variable (e.g., time or distance) and to associate each state with its derivative with respect to the independent variable. This latter could be by adding statements of the type:

INDEPENDENT VARIABLE t; STATE ppp HAS DERIVATIVE dpppdt;

where t, ppp and dpppdt are user names.

Issues that arise when making this type of declaration are:

• several models may declare the independent variable. It has to be the same variable for all models, requiring that it be aliased or a passed parameter from on high or some such thing.
• IF we want the integration equations (e.g., the bdf equations) to be ASCEND models, there has to be a way to associate these equations to the states and their derivatives. My approach in the .pdf document requires me to put these into vectors - not a nice thing to have to do as it requires me to identify in my modeling all the states and their derivatives and map them into these vectors. The strength to my approach is that I can have a variety of integration equations (e.g., bdf or RK4 or whatever) that I write and add to the system strictly as models and not as C code. This will be a sticky issue we need to discuss.

To restate this last point, if we have the DAE solver (which will be an SLV solver in C) collect the states and their derivatives by finding all declaraions and then associate integration equations with each of them, then how are these equations added to the ASCEND model so the compiler creates its solve system for the model plus these equations (and where these equations are my choice as I have written them as an ASCEND model)? Maybe all this just requires some clever ASCEND modeling that I have not yet thought enough about how to do??

-- User:Aw0a - 08 Jun 2006

I'm anti the "x HAS DERIVATIVE y" syntax. I think it's wordy and "y = diff(x)" or "y=der(x)" would be better (it eliminates the need for extra lines in the model declaration). Multiple-word syntax tokens considered harmful: we need to make the ASCEND language more compact (as evidenced by the fact that you felt it necessary to write a macro package for typing models in Emacs).

We have a key problem here: should ASCEND display an instance hierarchy that corresponds to the model that was declared by the user, or should it show its internal instance hierarchy? This latter, in the case of Art's approach will involve lots of addition curve-fit variables and derivatives, etc. I am in favour of showing only what the user has declared: What ASCEND displays should fit with the user's mental model, rather than the other way around.

-- User:Jpye - 13 Jun 2006

I am not for a more compact language. We were deliberately wordy, and I have never regreted it. ASCEND is nice because one can guess meaning from the words used. I like self documenting languages. APL (a long time ago) was so compact an expert in it would always misunderstand any program longer than a couple of lines, even if he had only written it a day before. One line of code could about jump through hoops. It did matrix operations with very compact syntax.

Are you suggesting using der(x) (for example) as the derivative of x with respect to the independent variable? y = der(x) only make y equal to der(x); it is NOT defined as der(x) using an equals statement - unless you overload the equal sign (I am against that).

I would have to think of the implications of using der(x) as the derivative of x. How would I give an initial condition to der(x) (which looks a lot like a function and not a variable in the model)? For steady-state, I would like to fix der(x) to be zero, for example. Ben also noted above that der(x) needs bounds, etc. If der(x) is to be looked at as a variable, then it has a different format from other variables (not good). When I think of DAE models, I really do think of the derivative of x wrt the independent variable as another variable in the model. So I do want it to be a variable that I can use in an equation (for example, to square it), bound it, give it nominal values to, give it initial guesses for, etc. My conclusion from this: der(x) has to be a variable, and der(x) is not a good syntax for it.

-- User:Aw0a - 13 Jun 2006

I have an idea that may do all the things for which we all have been asking to model and solve IVP DAEs. It is motivated by John insisting we use a form like y=diff(x) to relate the derivative of x wrt time to the variable y and Ben and I arguing that we do not use this form. This idea would NOT require us to modify the compiler, it would also allow one to use whatever model -- as an ASCEND model -- one wishes to do the integration, and it would NOT require one to run around the model manually to find all the derivatives and map them into a vector to pass to the integration routine. It should only require us to write a solver in C, as with our other solvers. It looks a lot like the statement that Ben has suggested and that John has suggested - yes, both.

Use statements of the form

``` dxdt IS_A speed;
x IS_A length;
t IS_A time;
```
``` derivx IS_A derivModel(dxdt, x, t);
```

I can use IS_REFINED_TO to upgrade any base derivModel to whatever form I wish it to have, where the refinement would bring in the actual polynomial equation to use - e.g., BDF or RK4. The solver would simply have to find all the instances of derivModel in a model to discover what are state variables when it is time stepping, integration error controlling, next poly order predicting and "stopping condition" checking. I could in principle have a different derivModel refinement for every state variable - a nice twist. So I could use different kinds of polynomials for different variables.

One issue: I wonder how we could enforce that all independent variables are the same one variable in the entire problem - the UNIVERSAL statement (seldom used of late in ASCEND) might help.

We can get away with this approach because integration polynomials only ever involve the state, it derivatives and the independent variable. There are no references to the values for other variables in these polynomials. I suspect this would also let us handle PDEs in a similar way.

The derivModel will set aside the space needed for past values of derivatives and states. The user's model, which includes this statement for each state, will define the physical relationship among the states and derivatives. derivModel will only supply the numerical integration equations.

To handle predicted variables that are not state variables (such as an algebraic variable whose value I wish to predict to aid converging at each time step or that I wish to use as a stopping condition), use a statement of the form

``` algegm IS_A algebModel(m, t);
```

and the system can collect these and use a different polynomial model to use for predicting.

I cannot see any obvious holes in this - yet. If I am not kidding myself on it, then it is actually a demonstration that the ASCEND language is doing a nice job in extending to this new type of problem with the features it already has.

-- User:Aw0a - 16 Jun 2006

I decided to see just how difficult it would be to create the models required for the approach I just suggested. I committed a new model into models/westerberg/ivpNondimensional/ivpNonNew that you may wish to examine.

Basically, for every state variable (i.e., variable y for which one writes an ODE to define dydt for it), include the following statement in one's ASCEND model:

```defineState_Velocity: stateVel IS_A diff(accel, vel, t);
```

There can be a variety of diff models in one's library, and one would use REQUIRE at the top of the model file to select which one wants at this moment.

This is really all there is to it.

The diff model adds polynomials of order 12 to the ASCEND model, for the current point and each of twelve past points, for both relating y to a polynomial in t and relating dydt to a different polynomial in t.

```These are the only equations one needs to add to provide the
```

numerical integration equations.

The solver will control which of these equations to include at the current time step by setting equation include flags and variable fixed flags appropriately. Mimicing how Lsode steps (as we would implement it), the solver has to choose the order polynomial to use for the next step. Let us suppose it chooses order 3. This means one will use third order polynomials. Only coefficients a to a would be nonzero. The solver would set the remaining, a to a, to zero and fix them. Based on whether it is integrating a stiff problem or a nonstiff problem, the solver has to choose whether to include y = poly(t) or dydt = poly(t) equations for the current and up to three past points. Both kinds of polynomials are in the model for all twelve past points. The solver would uninclude all those not wanted.

That is the essence of getting the integration equations into the DAE model directly for solving them with the user's DAE model equations.

The solver (written in C) would have code that corresponds to the methods I wrote to do the error checking, to decide the order to use next, to predict the next point, to check for stopping conditions, etc.

All in all it looks pretty straight forward and a lot simpler than my previous code. Adding the diff model for each state variable made the models a lot less complex. Also always adding 12th order polynomials and killing the higher order terms by setting coefficients to zero and fixing them made what had to be added easier.

## The current view

The current view is that the integrator API as currently shared by IDA and LSODE will more-or-less stand, but there will be more changes to the way in which integration systems are analysed. We currently propose to use syntax something like

```DERIV(dx_dt,x);
```

This will add a 'LINK' of a specific type ('time-derivative', perhaps) between these two variables, which the analysis routines will be able to find and use when preparing the model for solution by the integrator IDA or LSODE.

## Blocking in dynamics models

The following diagram is helpful in understanding the different types of variables and relations that are present within a dynamic model