User:LeonKwek

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.

Goals

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

• 25/05/2011: 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: 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
• 20/05/2011: Move the flagfn stuff over to the IDA integrator data struct to clean up the prepare_integrator functions : 10 min

Medium priority

• 23/05/2011: "Learn how to use matplot and run without the GUI" : 3 hours
• 20/05/2011: Fix the issue of Non incident var with non incident derivative : unknown

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

• Non incident var with non incident derivative thing
• Check if OTHERWISE to sort default boolean_var values is ok. Should the complier pick up if CASE list is non-exhaustive?

Completed Tasks

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

• 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

• 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

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

File:Example.jpg

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 libconsub3.so 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

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

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)
CASE(TRUE, TRUE)
END WHEN

or alternatively

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

WHEN(b)
CASE(TRUE)
END WHEN

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 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.