Jump to: navigation, search

Leon Kwek is an undergraduate student at ANU working with John Pye on adding support for integration of conditional models to ASCEND via modification of the IDA solver to add boundary-crossing/DAE reinitialisation support.

Browse my code: leon:


Implement code that completes the support for boundary-crossing in DAE integration:

  • Check support for boundary detection is complete/robust. (80%: Satisfactory for now)
  • Using LRSlv, solve the logical problem and determine if any logical states have changed (50%: Work around solution, want to fix later)
  • Determine if as a result any WHEN statements have changed state (100%: was already implemented)
  • If so, reconfigure the model: turn of disused equations, turn on new equations (100%: already implemented).
  • Determine strategy to solve the initialisation problem, possibly bearing in mind the variables/relations that have recently been added to/removed from the model (75%: Current strategy is to trash the whole IDA object and rebuild from scratch. According to Ben, not worth fancy footwork unless re-init is identified as a bottleneck).
  • Solve the initialisation problem either using IDA's built-in routine or else through direct use of QRSlv (50%: Currently using IDACalcIC. Works for simple models but not extensively tested. Not tried with index-switching problems yet)
  • Restart IDA and continue integration (75% No hitch with the complete trash and rebuild strategy).

Test the resulting code on a set of sample models given by Jan Huckleheim (in the first instance)

If possible: write some new solar energy simulation models that make use of the new feature, eg thermostat controller for a hot water system.

Current task list

(Date Added): (Task) : (Estimated time / difficulty)

High / immediate priority

Medium priority

  • 20/05/2011: Fix the issue of non-incident var with incident derivative : unknown
  • 21/06/2011: Remove the need to analyse twice at start up to set the equations correctly

Low priority

  • 14/04/2011: Figure out a better / the correct solution to the LRSlv messing with the compiler issue. Possibility: work with Ben to properly file a "design document" + follow through with the solution: unknown

Stuff to discuss with John

  • scons help with SUNDIALS_VERSION_MINOR
  • By reporting a point whenever we cross a boundary, the size of the data set is unknown until after integration. Might cause issues with plotting scripts?

Completed Tasks

(Date Added -- Date Completed): (Task) : (Estimated time -- Actual time) Comments:

  • 20/05/2011 -- 22/06/2011 : Get a conditional "bouncing ball" model running. Want to rigorously test stopping, reporting and restarting after boundary crossings Expect issues with IDACalcIC when using very fine time steps / boundaries that occur on or very close to time steps :2 hours
  • 25/05 -- 19/06: Sort the issue of IDA stopping short of boundaries, switching the system, then actually crossing the boundary and switching the system back to the wrong equation set.
  • 20/05/2011 -- 25/05/2011: Move the flagfn stuff over to the IDA integrator data struct to clean up the prepare_integrator functions : 10 min, 10 min
  • 23/05/2011: "Learn how to use python to drive ASCEND without GUI" : 3 hours, 3 hours

Learned a whole lot about the GUI in the process. Could probably debug the GUI so it work mores robustly with boundary switching when the time comes

  • 24/05/2011 -- 24/05/2011: Add boundary crossing hysteresis: no crossing a boundary twice in the same direction: unknown

This is implemented, but ended up not being the solution to the stopping short of boundaries problem I thought it might. EDIT 26/05 doesn't actually work as thought. Consider this problem not solved

  • 20/05/2011 -- 23/05/2011 Learn how to use matplot and run without the GUI

Turned out to be a non-issue. Can just use pylab, since I'm already familiar with matlab

  • 20/05/2011 -- 23/05/2011: Get a conditional charged particle model up and running to iron out the new boundary switching code

Model works (as far as I can tell) with use of "dummy" equations to work around the non-incident incident derivative issue. Issues with the pygtk GUI: sometimes it works perfectly and produces nice boundary switching plots, other times the integration fails with exactly the same starting conditions. Boundaries are sometimes detected twice, which causes the wrong equations to activate. This does not occur when driving ASCEND with CUnit or basic python scripts. Don't understand the operation of the python GUI enough to debug. Added the "Learn python, or at least enough to understand the GUI" task to figure this one out.

Progress Log

Aug 10 Merged with the latest trunk version to pickup various fixes with the GUI, all tests appear to be working fine.

Managed to track down mysterious inconsistent crashes following restart: the function integrator_get_ydot leaves passes null entries in the integratorsystem->ydot list. Hence IDA would end up with random values in it's ydot NV_vector depending on however N_V_new_serial assigns memory values. These were sometimes too large and were causing IDACalcIC to explode. Fixed for now by assigning 0 to the null values in the ydot NVector. Might be a better way though.

Merged with dante stroe's branch to pickup the link stuff. Difficult to test how successful the move was -- none of the with_LINK test models compile properly, I think it's because they aren't using solver_vars...I've managed to get my models working using DER(ydot, y) syntax in methods, but

DER(v, y);
DER(a, v);

is not working, have to setup a dummy variable ydot although Dante's comments everywhere say this shouldn't be the case. Hard to tell if feature not fully implemented or changes in the last two years have broken. Will have to ask john

2 Aug Meeting with John

  • Strange errors seen when running, no bounce is being observed! Why is that? Is the bouncing being turned off somehow?
  • Possibility of a virtual machine as another tool for tracking down these problems
  • Valgrind on the Cunit test code to track down and fix any memory leaks coming from the IDA stuff specifically (trying to ignore any other ones that don't seem important).
  • Because of work done to track down memory leaks, it might be worth merging the trunk code changes into your branch ("svn merge -r:OLD:NEW svn://trunk .")
  • Regardless of the bug-fixing effort, we need to finish this project... C code by itself is working (we hope) so let's use that as far as possible.
  • Follow up on the Jan Huckleheim sample models and think about any additional (thermal?) ones that we can add.

'Runs on the board':

  • boundary switching works (when you started, it was only boundary detection)
  • directional boundary detection implemented
  • forcing SATISFIED statements to use IDA's internal boundary state implemented
  • some demo models working as expected
  • developed some understanding about needed language features and types of problems that can be modelled

Things that we'd like to incorporate:

  • possibly incorporate Dante's work into Leon's branch
  • theme idea: expanding/improving support for dynamic models in ASCEND
  • a real live thermal model of something... to be discussed.
    • idea of a solar boiler, see Jose
    • idea of hot water system with dynamics introduced from data reader.
  • discussion of language features for switching and what others have done

22 Jun Had the bouncing ball model written up a while ago, but finally got it to simulate correctly today: Bball.png

the hitch was a sneaky combination of two errors, one of which has a certain random element that was confusing the beejeesus out of me. First though, I was stringing up a second derivative as

(* Some relations *)
v = ydot;      (* "dummy" equation for 2nd derivative *)
F_s = -1*k2*v; (* damping in the bounce *)

METHOD ode_init;
y.ode_type 	:= 1; v.ode_type	:= 2; 	
y.ode_id 	:= 1; v.ode_id	:= 1;
ydot.ode_type 	:= 1; a.ode_type 	:= 2;	
ydot.ode_id 	:= 2; a.ode_id 		:= 2;

using the above Ascend was complaining about being 1 variable sort of relations. After many hours of torture, switching the position of v and ydot in ode_init solved the issues (though I can't remember what inspired me to try that).

At the same time I was getting errors with "IDA:IDASolve: At t = 2.47421 and h = 3.8147e-06, the error test failed repeatedly or with |h| = hmin. (error -6)" and variables getting set above or below upper this point I think it has something to do with sharp changes at boundary crossings. Sometimes the error comes during IDACalcIC and sometimes during IDASolve, but the kicker is that the errors don't occur consistently. That is to say, can often just rerun the simulation without touching the settings and it works.

Also: getting different results (with the same parameters) depending on how I launch the tests! So far no errors when running the CUint tests through my ida_test.c and the afore mentioned inconsistencies using the GUI (GUI also crashes when I try to open the model using the drop down menu, have to call from the command line).

Finally, trying to drive using a python script, which is my vehicle of choice for its plotting capability, I get consistent errors with variables set above their bounds when integrating with exactly 62 or greater time steps. 61 steps is ok.

My sense is that the explanation of what is going on here lies somewhere in the guts of either the function evaluators or in IDA's routines, neither of which I'm especally inclined to dig into, will wait for comments from John. In the mean time however, I'm out of short term goals, which this sort of testing and bug hunting. I'm going to look into trying to model a slightly more complex conditional system to keep things interesting.

21 Jun Spent hours last week trying to get to the bottom of the SUNDIALS_VERSION_MINOR business but again to no effect. Decided in the end to do a complete fresh install of ubuntu 11 to go along with the installation of a new CPU. On the new system with the sundials 2.4.0 package the scons TryRun of IDA now fails altogether. The test works fine when I run it separately from the ascend built environment but I haven't been able to pin down the issue. I tried changing the sundials_version_text test code to just main(){return 0} but no dice.

Everything compiles and runs fine with the check disabled, so I'll leave it at that for now. IDASetRootDirection does the job fine and issues with double crossings seem to be fixed. Noticed some strange results while testing however: found the issue stems again from problem analyzing(.c) without conditional systems in mind. Boolean vars get assigned default to true so that the initial analysis of the WHENs will generally choose the wrong set of initial equations. As a temporary workaround can just do a LRSolve if boundaries are present in the model and re-analyse if some_disvars_change. A better solution will come about along side dealing with the non-incident var issue, which will have to deal with the deeper design issue of somehow flagging the model as conditional during analysis.

Also begs some design questions about "assuming IDA is always correct" when it comes to boundary crossings. If the wrong set of equations mange to slip in at some point, then IDA happily continues through the whole integration while sending back nonsense data. The user has to have responsibility for running sanity checks on any results they get, but as my experience with the wrong initial equations has shown, these errors can be very hard to pin down because the output looks otherwise reasonable.

8 Jun Discussion with John, Leon. Agreed that the IDA boundary direction function from SUNDIALS 2.4.0 is needed to resolve issue of possible repeated boundary detection. Need to resolve issue of stray values of SUNDIALS_VERSION_MINOR from SCons as part of that. Decided that the problem of solving systems with non-incident x but incident x_dot will need revisiting of the system/analyze.c code; there may not be an easier way. So we will leave that issue for later. Aiming to get the bouncing ball model solved for next week!

26th May Spent the whole day unsuccessfully trying to track down why my SUNDIALS_VERSION_MINOR remains equal to 3 after updating. Scons appears to detect the version change, namely if I print minor right before the line

context.env['SUNDIALS_VERSION_MINOR'] = minor

in the SConstruct I get minor = 4. I did a full purge of both ASCEND and sundials on my system to see if I had any old copies of libascend or libida_ascend lying around, but to no effect. Checked all env variables...checked that the correct headers are referenced. I'm out of ideas for the moment, will move onto "Non incident var with non incident derivative" instead in the mean time.

25th May 5pm: IDASetRootDirection was introduced in version 2.6.0 of IDA, so we'll have to upgrade, and place version checks / warnings to ensure some backward compatibility. Completed the "Move the flagfn stuff" task without issue

3pm: Aha! Turns out IDA already has a function IDASetRootDirection. Obvious that it should exist now that I know about it. Assuming this works as advertised, should only take 20 minutes or so to have the switching algorithm rock solid.

Decided to play around a bit while thinking of the solution to the double switching:

Charged particle.png

A charged particle crosses from an electric to a uniform magnetic field. The magnetic field is then "switched off" at t = 35. Used a python script to run the model, and a little reporter was setup to record the data. Plot using matplotlib. Proof of completing the "Learn how to use python to drive ASCEND without GUI" task.

24th May While playing with the python wrapper, discovered the cause of boundaries being "double crossed." As anticipated before, under the right circumstances IDA will tret at 4.999999995 for a boundary at x < 5, then cross the boundary again. Obvious solution is the hysteresis thing: dis-allow boundary crossings in the same direction twice in a row. Implementation design: Add array length integ->nbnds. Flag +1 / -1 when a boundary is crossed increasing or decreasing (or 0 if never crossed). Issue: how to detect which direction a crossing happens? Added new high priority task to add this feature

All right this solution is implemented useless durr: see explanation by way of example: The boundary is y < 5. IDA wpproaces the boundary while integration, stops at y = 4.999999999979 and boundary rebuild is triggered. Logically we are on the other side of the boundary and just made a negative crossing. At the immediate next step, move from y = 4.999999999979 to 5.000000000012 detect root again. Switch states, making a positive crossing. Now we are in the region y > 5 with the wrong set of equations.

If the boundary variable was time we could just move IDA onto the boundary, but that might have adverse effects on derivatives...hmmm. Still thinking of a solution. Going to leave the positive / negative crossing checks in for the moment since they appear to work and may be useful later.

23 May: Completed task "charged particle model", but as intended, uncovered issues with the switching process.

20 May, later that day....

Attempting once again to setup a conditional model of a charged particle accelerating through a uniform electric into a uniform magnetic field. Last attempt thwarted after trying to use .ode_type = 3 for high order derivatives.

Issues encountered:

  • BUG: Seg fault when trying to clean up integrated models without boundaries. Caused by trying to ASC_FREE the array of stored boundary states (bnd_cond_states in ida.c). Straight forward fix
  • Both the electric and magnetic field equations (and I suppose a whole lot of F = ma models) make no reference to the position variables, x and y in this 2d case. In this case x and y are dropped from the var list. Ida analyse picks this up, error message is "Non-incident var with an incident derivative. ASCEND can't handle this case at the moment, but we hope to fix it." The time to fix this is probably now (depends how difficult...) Short term solution: add a dummy variable and equation x = xdummy, y = ydummy. ASCEND seems to accept the dummies and appears to integrate ok. Obviously not desirable though.
  • boolen_vars were defaulting to True, which causes problems if that happens to be a CASE that doesn't occur in the model so you don't write it in. Learnt about the OTHERWISE keyword, did the trick.
  • BUG: Was writing output with tret and yret as null vectors when skipping over a time step if a boundary occurred too close.

Progress / Features added

  • IDA now outputs data exactly at the point of a boundary crossing, since this is probably of interest to the user.

Secret undocumented coding April 14 -- May 20

  • Implemented the complete teardown/rebuid of IDA at boundary crossings. Re-arranged ida.c to split single use setup stuff from that which needs to be repeated every time IDA is rebuilt. Broke the process into logical parts (malloc, set parameter inputs, calc IC etc.) then stuck in sequence in a new function ida_prepate_integrator.
  • Removed the has_crossed flag since the IDA root detect flag works OK
  • Added a check to see if a boundary crossing is detected close to (or on) the current requested time step. Necessary because we have to feed the value to IDACalcIC.

20 May discussion with John, Leon

  • Implemented first option below, "double-purpose the current 'perturb' stuff using a special value of 'perturb_value' to indicate the IDA usage, and pass a list of boundaries plus IDA's reported truth values through in place of the perturb list."
  • Testing that integration is working
  • Some issue about IdaCalcIC going wrong when initialising too close to a recently-passed boundary. A time-step is needed by that routine, and we have to guess(?) a reasonable value of that.
  • Some glitchy behaviour possible occurring just before (?) or at a boundary when plotting the response of the 'boundaries' model.
  • Perhaps try to do a bouncing-ball model as a clearer way to test.
  • In ida.c, code that passes flagfn pointers can be simplified by adding fields in the IntegratorIdaDataStruct that is accessible via the 'integ' object.

14 Apr discussion with John, Leon

  • want to avoid letting LRSlv evaluate any boundary residuals, since IDA always knows their value already and will tell us if they have ever been crossed.
  • IDA actually doesn't know the truth-value of a boundary expression, so we need somehow to establish the truth-value before integration commences.
  • we anticipate needed to support crossing of multiple boundaries at once (simultaneously, as far as IDA is concerned). IDA can provide a list of boundaries that has been crossed, so we need to be able to deal with that returned information correctly.
  • CMSlv drives LRSlv to 'test' different truth values on boundaries. It does this by flagging boundaries that it wants to 'imagine' crossing, and LRSlv, via 'perturb' code in the SATISFIED evaluation routines, returns the altered truth value of boundaries where requested.
  • IDA in the other hand wants to drive LRSlv without letting LRSlv evaluate any of the boundary residuals. So we need a way of making SATISFIED evaluations refer back to some external data instead of visiting the relevant relation and evaluating it directly. Possible options:
    • double-purpose the current 'perturb' stuff using a special value of 'perturb_value' to indicate the IDA usage, and pass a list of boundaries plus IDA's reported truth values through in place of the perturb list.
    • remove existing perturb stuff and pass a function pointer plus data pointer through to the evaluation code, and allow the solver to implement its own boundary evaluation method. CMSlv+LRSlv would use this to artificially tweak the design boundary(ies) and IDA+LRSLv would use this to look up the desired value from some kind of list passed to LRSlv by IDA (generated by looking at the flags attached to bnd objects).
    • write a parallel set of evaluation functions that work differently, with no change to current CMSlv-related functionality.
    • do some dirty tricks using direct calls to ascend/compiler code from the solver, and flag boundary relations with their 'forced' truth-values in the case where they are being evaluated via IDA+LRSlv.
    • implement another dirty trick in which a global variable forces switched behaviour whenever an attempt to evaluate a SATISFIED statements comes along. Or is it possible to store some flag in some other high-level object in the compiler? We would still need a way of passing the required boundary states through to the compiler, as well.
    • change LRSlv so that it never evaluates any boundaries, and instead ALWAYS picks up the value of the residual and uses that to give the boundary's truth value without any new real-valued calculation occurring. We would then need to change CMSlv so that did the required evaluation before calling to LRSlv, and IDA would need to poke the needed values into to the boundary relations before calling LRSlv as well.
  • If we modify LRSlv then we need to assure that CMSlv is not broken. Leon has a copy of CONOPT now for testing CMSlv still works. Models that you can check with are at Conditional modelling#Examples.
  • Compiling CONOPT needs gfortran; you should have a at the end of the build process. Run 'scons install' as root to install to /usr/local/lib which ASCEND should then pick up.

29th March

logical states A way of flagging boundaries so that LRSlv looks up the desired state of boundary condition (as opposed to evaluating) in the event of a crossing detected by IDA. The issues described below was dealt with by adding flags to the compiler side relation structure -- most definitely not an ideal solution. Looking to seek advice from Ben Allan on the dangers of this move / some other way solver-side around.

WHENs A function written for CMSlv seems to handle everything we need for our purposes here, however I haven't yet found an efficient way to get variables to update once the appropriate equations have been switched on and off. Calling QRSlv does the trick, but that seems like overkill.

The bigger challenge at the moment however -- reinitialising IDA given that the number of variables in the problem may have changed as a result of a crossing. Ben has suggested that the smart approach is to simply bang up the whole system from the beginning and forget about nifty optimisations until they become glaringly necessary. I'm partial to this angle, largely because it's the most straight forward and obvious but also because it's intuitive and therefore robust in that it's harder to make programming mistakes.

9th March

logical states Decided to abandon the perturb stuff as on reflection it refers to a different sort of problem. When numerically integrating we always certain which boundary we are crossing and where we will end up after the event.

Current understanding of the logical problem: we need to recognise an arbitrary setup of logical states. For example, say we want to model a valve that opens when both a given mass flow rate exceeds a value m1 and a temperature reaches T1

c1: m > m1;
c2: T > T1;

Then we want to be able to express opening of the valve equally as

b1 == SATISFIED(c1, 0.001{kg/s})
b2 == SATISFIED(c2, 0.1{K})

WHEN(b1, b2)

or alternatively

b == ( SATISFIED(c1, 0.001{kg/s}) && SATISFIED(c2, 0.1{K}) )


Now for the latest issue: LRSlv sorts the boolean variable simply by calling the slv_common function slv_direct_log_solve. This in turn uses ASCEND's internal relation and logical relation evaluator, which can handle the required branches of arbitrary length. The scope of these functions however, is such that we don't have access to the boundary or logrel structures, either of which would allow the flagging of boundary crossings detected by IDA.

7th March

boundary detection Issues with multiple detection during integration found to actually be the result of the incorrect assumption that IDA's 'internal' independent variable always matches the requested return value (tret). Re-initialising IDA with the desired start time solves this issue.

logical states Seems that simply calling LRSlv's solve function bypasses the boundaries and evaluates logical relations directly. This is to be expected -- it is a logical not boundary solver -- however this means that we cannot get the desired switch in conditional relations without modifying LRSlv. All this worry may be for nought however. The following from logrel_util.h :

*  The easier way to test if the crossing of a boundary affects the current
*  configuration is to solve the logical relations assuming such a
*  crossing and finding if the configuration of the problem changes.
*  This task is performed here by using the "perturb" flag and a list
*  of pointers to the relation instance which constitutes the boundary(ies).
*  The activation of the perturbation flag causes the calculation
*  routines to INVERT the truth value of the SATISFIED terms containing
*  the relation instances pointed by the list of pointers.
*  The flag and the list of pointers are needed because the SATISFIED term
*  is not an instance by itself and we have access to it only trough the
*  logrel instance.

As is to be expected, CMSlv already has a means for handling this issue, though I can't make sense of how it works at the moment.

23rd Feb (at revision 3253)

boundary detection IDA's built-in root finding engine seems to do the job in the case of boundary crossings of independent, dependant, derivatives and combinations of the three. Existing code suggests that IDA will return a positive or negative flag for a respective increasing or decreasing root however I cant get this to work.

ideas for further testing: minimum separation of boundaries IDA will detect. Speed of the rootfinding.

Logical States Setting LRSlv to work is straightforward however there is the issue of setting the model parameters for the logical solution. For example, IDA may indicate a root at t = 3.99983 for the boundary t >= 4, in which case calling LRSlv immediately does not trigger a change in the appropriate boolean variable. I think IDA also returns a minimum step size used during the integration which could be used to 'bump' a variable over to the correct side of a boundary if the issue of increasing/decreasing roots described above can be resolved. What about bumping multi variable boundaries?

Another thing to consider: boundaries that occur on the time samples of the integration loop are picked up twice by the rootfinding algorithm. If a root occurs on or very close to the next time step tnext, upon restarting IDA will need to run to tnext +1 before back to the top of the loop. Notes on handling conditional branches from the documentation of a similar modelling program[1] mention that boundaries that depend only on the independent variable have the bonus of being known ahead of time. In that case it is possible to handle the discontinuity by arranging for an integration time step to fall on the boundary and switching the model appropriately. However the need to control the time step probably adds more complexity than it saves (especially for nested components).

WHEN statements have changed Will look next into how the WHEN conditions are handled.


  1. Elmegaard and Houbak, 2004, DNA – A General Energy System Simulation Tool, Technical University of Denmark