User:Ksenija

From ASCEND
Jump to: navigation, search

Ksenija Bestuzheva
Development branch: ksenija:, ksenija2:

Student in the State University of Management (Moscow), studying applied mathematics and information technology.

I have some knowledge of mathematical analysis, linear algebra, complex variable theory, differential equations, mathematical statistics, control theory. I have experience with C, C++, Java programming through my studies.

GSOC2013

Possible ways of improvement of our conditional syntax

The first possible modification is defining the condition right in the event statement:

EVENT (x > 0)
...
END EVENT;

Some details like tolerances could be set, for example, inside brackets. But I also like the current definition of events and whens using variables. So we may leave the whole structure (cases, etc.) of the statement as it is now, but allow not only variables, but also relations as conditions. If we decide to do this for events, then, I think, we should do the same for whens.

The second is defining relations inside events:

a IS_A submodel1;
...
EVENT(...)
  CASE TRUE:
    USE a;
    x0 = x;
END EVENT;

This would make the models shorter but still allow submodels referenced by event. In this case it seems that not so much changes in the implementation would be needed. And, as for the previous suggestion, if we do this for events, the same should be done for whens, I think.

Example models

ksenija2:models/ksenija/test_event2.a4c - the simplest model with events. When x crosses some given values, it increases discontinuously, and then integration continues in a usual way.

ksenija2:models/ksenija/thermostat.a4c - a model of a thermostat. The heater is swithed on when the temperature rises above 23 degrees and is switched on when the temperature falls below 19 degrees.

ksenija2:models/ksenija/bball_event.a4c - a bouncing ball. When it hits the ground, it instantly changes its speed.

ksenija2:models/ksenija/bball_event3.a4c - another version of the bouncing ball model.

ksenija2:models/ksenija/hysteron_event.a4c - a model of a simple hysteron. This is an example from Krasnosel'skii, M.; Pokrovskii, A.. Systems with Hysteresis. A piston is moving inside a cylinder, the position of the piston is the input u, the position of the cylinder’s right end is the output x.

A sketch of the cilinder


ksenija2:models/ksenija/hysteron_event2.a4c - one more model of a hysteron. This is an example from Krasnosel'skii, M.; Pokrovskii, A.. Systems with Hysteresis.

ksenija2:models/ksenija/simultaneous.a4c - a model with two boundaries crossed simultaneously.

ksenija2:models/ksenija/vessel_with_siphon.a4c - a model of auto-oscillations of fluid level in the vessel with simultaneous inflow and outflow from Ju.I.Neimark, "Mathematical models in natural science and engineering".

A sketch of the vessel


ksenija2:models/ksenija/clock.a4c - a model of the Galileo-Huygens clock from Ju.I.Neimark, "Mathematical models in natural science and engineering".

Project plan

  • Community bonding period
    • May 23 – June 17. Discuss the project in a more detailed way. Think over the plan, discuss some complicated cases like nested events, events inside whens, defining the direction of the boundaries, etc.
  • Parser changes
    • June 17 – June 25. Create the parser rule and the structure for storing the new statement. Write code for verifying the new statement, create the event type description, generate the event names.
    • June 25 – July 8. Write test models and test cases, fix bugs (if any).
  • Changes in the instantiator
    • July 8 – July 13. Process the event statement in the instantiator. Make all checks which can't be done before. Make sure that the needed variables and relations exist.
    • July 13 – July 24. Create the event instances. Find all needed conditions, create pointers between relations, assignments and events. This may require a new compiler pass or may be it is better to do it in one of the already existing passes (together with WHENs, may be?).
    • July 24 – July 27. Finish the work with the compiler. Make sure that all cases are processed correctly.
    • July 27 – August 2. Test, fix bugs (if any).

The needed changes in the compiler should be done by mid-term evaluation.

  • Changes in the problem analysis function
    • August 2 – August 16. Create the lists of events in the problem data structure. Add the needed flags to relations and assignments. Write code for switching the flags on and off at the boundaries.
    • August 16 – August 20. Test, fix bugs (if any).
  • Changes in the solver
    • August 20 – September 6. Change the boundary crossing function. After solving the logical relations go through the events list and change the flags using the written before function. Solve the system and reset the flags. Some more checks for system consistency are likely to be required.
    • September 6 – September 16. Test, fix bugs (if any). Add more models which use the new functionality.
  • Completing the project
    • September 16 – September 23. Final testing, writing documentation.

Weekly reports

September, 09 - September, 15.

This week I have been working on GUI features.

The points where an event was triggered are marked red on the plot. Plots for each event can also be created showing how the variable changes at each step of processing the system at the boundary.

Several equal values of the independent variable mean that there was an event
Red points indicate events
The user may choose an event and plot the values found at the chosen boundary


September, 02 - September, 08.

  • Changed cond_config.c so that when the OTHERWISE case of the event is switched on at the boundary, the method associated with the event is not called.
  • Switched off event analyzing in configure_conditional_problem (which is called during problem building). The events are analyzed later, so there is no need for that call.
  • Set the event_active flags FALSE before analyzing the system at the boundary.
  • Report errors if the analysis of the system after solving the boundary equations fails.
  • Until this week the integrator was reinitialized after each root found. This didn't influence the speed much, but for some models this reinitialization could lead to incorrect solutions. Without it when solving the hysteron model IDA reported an error because of two very close roots.

To solve this problem I changed the boundary crossing code. Now after solving the boundary equations for the first time it is checked if values of some discrete variables have changed. If yes, the system is analyzed one more time. These actions are repeated in a loop until no events are called or no discrete variables change their values. So the case when an event causes another event is handled properly and no unneeded reinitialization of the integrator is done.

  • Thought about boundary emission and read about the implementation of this feature in Modelica and Mathematica.
  • I decided to spend the last week of GSOC on improving the plotting tools for models with events. There are several points at each boundary which can be useful for debugging of the models, and I am going to plot them in such a way that the user would be able to see those points and distinguish them from other points.

August, 26 - September, 1.

  • Removed some unneeded code.
  • Changed the ida_cross_boundary function so that previous of discrete variables values are now corrected only if a double crossing has been detected.
  • Added a warning message which is shown if relation instances are referenced by both events and whens.
  • Set the flags indicating if there is a method associated with the event. If there isn't, the search for the method is skipped.
  • Changed the integrator_ida_solve function crossing one and the same boundary twice in one direction is avoided not only if IDASolve is recalled up to the same tindex (which happens if the root is found far from the desired output t), but also if the timestep is advanced.

August, 19 - August, 25.

  • Implemented one more (experimental) approach to event handling. Mostly because of the problem with unwanted multiple boundary crossings in option 2 I think that the better approach is option 1. Here is the description of the two algorithms and their comparison:

Option 1.

This code uses flags set by IDA for solving the systems of logical relations at boundaries. Except when finding the initial values of boundary flags, for this algorithm there is no difference between strict and non-strict inequalities.

The advantages of option 1:

The most important: currently this algorithm is more reliable. The main problem with the second option is that in some models the system gets stuck at some boundary, crossing it again and again because of small unwanted variation in the values of variables on which the boundary depends. This may result in wrong solution of the model or infinite loops. This can be avoided by some changes in models (for example, excluding the equation which causes the unneeded variation from the boundary system), but option 1 works without such changes too.

It is a bit more fast than option 2, though I can't say that it is easy to see the difference in speed.

Option 2.

With this option LRSlv doesn't use the flags set by IDA and finds the values of boolean variables only by current values of variables in the model.

The advantages of option 2:

In option 1 there is an assumption that if the boundary is crossed, then after solving the system this boundary won't be crossed in the opposite direction. But this assumption makes the code more reliable and helps avoid the problem which I described in option 1. And as far as I understand, all models can be written in such a way that they will satisfy this assumption (for some of them the conditions will need to be a bit modified, but it is possible and not so difficult for the user).

In option 2 strict and non-strict inequalities are different. This works in the following way: when IDA report a root found, only non-strict inequalities are satisfied. ida_cross_boundary is called, all necessary actions are performed. But strict inequalities at this point are not satisfied. So IDA makes one more very small timestep, and ida_cross_boundary is called one more time. These additional calls of ida_cross_boundary make option 2 slower. But the question is: is it correct to trigger the events and/or switch to a new region after this small timestep, but not at the point where IDA reported a root? On the other hand, the advantage is that option 2 enables the user to set the order of events/regions. For example, if some region should be switched on before the event is triggered, and not after. But I haven't yet seen such models in which the order matters.

  • Added a button for viewing PreInfo to the var properties window.
  • Added one more version (the one from the Modelica specification) of the bouncing ball model: ksenija2:models/ksenija/bball_event3.a4c. It works fine, so now the problem with the bouncing ball may be considered solved.
  • Fixed some bugs in event handling:
    • In the case of double crossings two calls to log_solve result in previous values of boolean variables being equal to their current values, while actually they may differ. So before solving the logical system the values of discrete variables are saved in an array and then the previous values are restored.
    • Reset the boundary flag not only for the last processed boundary, but for all boundaries, because multiple boundaries may be crossed simultaneously.
    • Flag the boundaries which are crossed as a result of solving the system at the boundary.
  • Fixed a problem with wrong initial evaluation of boundaries.

For example, in the bouncing ball model the initial speed is 0. The model contains a conditional relation v >= 0. Before starting integration this condition is evaluated as true. The boundary flag is 1. The ball starts falling, the speed is negative, but the boundary flag is still 1.

To fix this problem I created an array of flags indicating if the boundary still needs to be evaluated. At the next timestep an attempt to evaluate the flagged boundaries is made, and this is repeated until all boundaries have the correct flags. This code is more reliable when the boundaries are not crossed during the first timestep.

August, 12 - August, 18.

  • Thought about possible ways of improving the conditional syntax. After trying several options of the parser rule found the one which didn't cause conflicts. Thought about how the data structures will need to be changed so as to handle new syntax. (All this was before I received a letter from Ben saying that I shouldn't do it now).
  • Started working on reinit() syntax. The idea was that the reinit(x,expr) operator would evaluate expr to a value, make x unknown and introduce the equation x = value. I added a parser rule and a new flag to the relation statement structure. Then I spent some time (not so little, I should say) reading the code, especially that deals with relations, and thought about the implementation of reinit. I came to a conclusion that reinit seems to be rather ambiguous. Currently I think that in models it can be easily replaced by assignments or a combination of boundary equations and FIX/FREEs.
  • So I moved on to calling methods from events. Until this week they were found by standard names event and end_event, but they couldn't be associated with particular events. Now they are found using the name of the event (for example, if the name of the event is my_event, then the method called before solving the system at the boundary would be called my_event, and the method called after solving the system would be called my_event_end). But I am also considering leaving also the old approach because in some models the same methods need to be called regardless of which event is triggered. May be some optimisation of calling the methods is possible (need to think about it). One more question is: how to associate events declared in loops with methods?
  • Made some minor changes in cond_config.c (where the events which need to be triggered are found and the system is reconfigured).
  • Rewrote one more model using the pre() syntax.

And finally I went to the country for the weekend and got lost in a forest. And that forest was a bit more frightening then, for example, a forest consisting of undirected cycle-free graphs :)

July, 29 - August, 11.

During these two weeks I have been working on pre() variables. Here is the description of their implementation:

  • When creating the type definition of the model the symchars for pre() variables are automatically generated using the name of the argument. The type of the argument is found and the pre() variable is given the same type. pre() variables are added to the model childlist.
  • At the same time some checks are done: it is checked that pre() variables are contained only in the autogenerated IS_A statements, relations and ALIASES statements; that there are no nested pres (because there is no sense in variables like pre(pre(x))); that there are no pres inside derivative variables (but derivatives inside pres are ok).
  • The structures for real atom instances contain a structure PreInfo which stores the information about pre() arguments (for pre() variables) and about pre() variables (for variables for which pres are declared).
  • Implemented merges of variables which have pres. If such variables are merged, their pre() variables are also merged. Also pre() variables are refined if their arguments are refined.
  • In the analyzer a list of pre() variables is created in the problem_t structure. Each pre() variable has the pointer to its argument. It is checked that pre() variables are incident only in those relations which are referenced by an event. All pre() variables are automatically fixed. At the end of problem analysis the list of pre() variables is moved to the diffvars structure which is later used by the integrator.
  • On the solver side everything is rather simple. At each boundary the value of the pre() variable is updated and becomes equal to the value of its argument. All the other time pres remain unchanged. So the following relation: x = pre(x)+1 solved at the boundary means that the new value of x (the value of x after the boundary) will be equal to the value of x before the boundary plus one.
  • I have updated the following models: the model of a bouncing ball and one of the test models.

July, 22 - July, 28.

  • Changed the observer so as to observe boolean variables as well as real variables.
  • Added two more models: one more model of a hysteron and a model of auto-oscillations of fluid level in a vessel with a siphon.
  • Disallowed integers and symbols in the list of event conditions.
  • Fixed the problem with non-active variables having active derivatives.

Started working on pre() variables:

  • Remembering some problems with implicit declarations of der() variables I am going to make declarations of pre() explicit, like declarations of derivatives.
  • Created a new structure inside a Name structure, a structure for the new statements and wrote all the needed functions and defines for them.
  • Added a parser rule for a pre() variable and for a statement declaring pres. This statement automatically generates IS_A statements.

July, 15 - July, 21.

  • Added a model of a bouncing ball which instantly changes its speed when hitting the floor.
  • Added a model of a thermostat.
  • Wrote code for solving the boundary equations before integration if needed. If the boolean variables match the values in an event case before integration is started, the event will be triggered.
  • Added a model with two boundaries crossed simultaneously.
  • Fixed some bugs.
  • Read more about a model with a sliding mode.
  • Searched for new models with interesting kinds of boundaries. (Yu.I. Neimark, Mathematical Models in Natural Science and Engineering; M. Krasnosel'skii, A. Pokrovskii, Systems with Hysteresis.)
  • Read about the Pantelides algorithm. (http://jpye.dyndns.org/pantelides/; E. Hairer, G. Wanner, Solving ordinary differential equations II. Stiff and differential-algebraic problems)

July, 8 - July, 14.

This week I have been working mostly on solving the model of a hysteron.

  • The IVP is solved only when the function which analyses events returns 1 which indicates that some events have been activated.
  • Changed the code which sets the boundary flags (indicating if the boundary is satisfied) in order to handle crossing one and the same boundary twice in one direction.
  • The boundary flags are reset only after processing all events and whens at the boundary, so if the boundary is satisfied at the moment when the root is detected, it will be considered satisfied while solving all the boundary equations. Will need to think if it is ok to make such an assumption.
  • Fixed some bugs.
  • Added more test models and wrote a test for integrating models containing events.

Now the hysteron model is solved correctly.

July, 1 - July, 7.

  • Added some more checks and warnings.
  • Added a test model of a simple hysteron using the EVENT syntax.
  • Wrote code for finding events which should be activated. First it is checked if some discrete variables have changed their values. If yes, we go through the list of events, check if the discrete variables which are used by these events have changed their values. An event is activated only if its condition changes its value. Then the values of the discrete variables are compared to the values listed in the case of an event. If they match, the event becomes active.
  • Started working on the solution of the test model. After some events become active the system is solved with QRSlv and then again with LRSlv. The first two or three boundaries are crossed correctly (and one of them contains an event) but then a wrong region is chosen and after the next event the solution becomes incorrect. In this model boundary behavior is mixed with region, so I decided to put the solution of this model off till some simpler models work.
  • So I added one more test model. It may resemble a thermostat controller but the equations are arbitrary. The only requirement for the functions was that one should increase and the other should decrease. When I tried solving the model three events were triggered correctly and then the values of boolean variables became incorrect and since then no events were activated. I found out that the reason was the following: the decreasing function stopped right at the boundary. An event was triggered, the decreasing function was switched off and the increasing one was switched on. The value of the output variable was already above the boundary, but the boundary flag showed that the boundary was only crossed from above to below. The next time it was crossed the boolean variable took such a value as if this time the boundary was crossed from below to above, which was wrong. So I added a flag to the boundaries indicating in which direction it was last crossed. Using the value of this flag and the value returned by IDA's function IDAGetRootInfo I check for crossings of such kinds and now all the boundaries are processed correctly. I will commit as soon as I put this code in order.


Error creating thumbnail: File missing
A screenshot of the plot for the test model


June, 24 - June, 30.

  • Finished the work on the compiler part: arrays of EVENTs are now created correctly; I changed CopyInstance so as it handles instances containing events, wrote code for merging instances referenced by events, changed some switches to handle event instances.
  • Changed the C++ layer so as to handle event instances. Events are shown correctly in the PyGTK GUI now.
  • Added one more model and some asserts and output to the test.
  • Added the OTHERWISE case to events which will define the model behaviour at continious-time frame when no events are triggered.
  • Started working on problem analysis. Created the needed data structures and added events to lists in the structures which represent the problem being analyzed. Added flags indicating if the relation/logical relation/etc. is referenced by some events.
  • Added a simple test for the analysis part.
  • Wrote code for analyzing events in the configure_conditional_problem function which is called at the end of problem analysis. Wrote code for passing the events lists into the slv_system_t structure.
  • Events inside WHEN statements are switched on and off in configure_conditional_problem depending on the values of logical variables used by WHENs.

June, 17 - June, 23.

  • I created a new structure StateEVENT for storing the information about the statement and wrote the functions for working with this structure;
  • wrote a new parser rule for the new statement;
  • created a new type EVENT (events will be implemented as instances so they should have a type description);
  • created a new type of instance for events EventInstance;
  • added a new field to structures of those instances which can be referenced by events. This field will store pointers to all events which reference the instance;
  • wrote code for instantiating events. It includes the following: checking the events and making sure that all needed instances exist, creating an event instance and creating lists of pointers to relations, submodels, etc.
  • wrote code for destroying event instances.

Final report

Event syntax.

New event syntax has been implemented:

event_name: EVENT (bvar)
 CASE TRUE: USE rel1;
 OTHERWISE: USE rel2;
END EVENT;

bvar is a boolean variable, rel1 and rel2 are relations or logrelations or submodels. When bvar changes its value from TRUE to FALSE, the corresponding case is activated. Then the model is solved with LRSlv and QRSlv, and after that rel1 is deactivated.

rel2 is switched on during the continious-time frame or at those boundaries at which this event is not triggered.

If the model contains methods with names event_name and/or event_name_end, these methods will be called before and after solving the system at the boundary correspondingly.

New boundary crossing algorithm.

Here is the algorithm of the boundary crossing function:

  • For each boundary:
    • If this boundary has been crossed:
      • Has it been crossed twice in one direction?
      • No:
        • change the boundary flag.
      • Yes:
        • if it has not been done yet, save current values of discrete variables.
        • Change the boundary flag.
        • Call the logical solver (it could be called not it this loop, but it would make the code more complicated and it would not make the algorithm much more efficient because double crossing of two boundaries at the same time is rather improbable. In all other cases there is no difference.)
        • Change the boundary flag one more time (so the value of the flag will be the same as at the beginning of the function.
  • Call the logical solver.
  • If there has been a double crossing, set the correct previous values for those variables for which it is necessary.
  • Have some discrete variables changed?
    • No:
      • reset the boundary flags, continue integration.
    • Yes:
      • reanalyze the system: switch on the correct regions (WHEN cases), call the methods associated with events, switch on the boundary equations.
      • Solve the system with QRSlv.
      • Check if some boundaries have been crossed as a result of solving the system with QRSlv. Set the correct boundary flags.
      • Solve the system with LRSlv.
      • If some discrete variables have changed, repeat from reanalyzing the system.
      • Reanalyze the system: switch on the correct regions, call the methods associated with events, switch off the boundary equations.
      • Reset the boundary flags, continue integration.

Here boundaries are data structures corresponding to the relations declared inside CONDITIONAL. Boundary flags can have two values: 0 and 1. If a flag is 0 (or 1), this means that when solving the system of logical equations the statement: SATISFIED(cond_rel), where cond_rel in the boundary condition, will return 0(or 1). The flags are set at the beginning of integration and changed at boundary crossings.

Pre() syntax.

The aim of pre() variables is to access the value of the variable before the boundary. Pre() variables are FIXed and change their values only at boundaries. An example:

y IS_A solver_var;
PREVIOUS y;
bnd_eq: y = pre(y)*0.9;

Variable structure modeling in ASCEND. Plan for summer of code 2013 and my final-year project

The goal

General description of the problem

Currently ASCEND solves problems with region switching behavior pretty well. But it still doesn't have good means for describing and solving models with changes happening exactly at the boundaries (which may also be one-sided - e. g. hysteresis). These may be solving some boundary equations or assigning values to some variables.

Example problems

Some example problems are given on this page: Event handling. Some more interesting examples may concern encoding analog signals into digital signals (see Delta-sigma modulation on Wikipedia). I am still searching for more examples.

Syntax

Requirements

The syntax should enable the user to 1) state which actions should be performed at the boundary 2) associate them with an exact boundary 3) state in which direction the boundary should work (for one-sided boundaries). And, as usual, the syntax should be unambiguous, flexible and clear.

Proposed syntax.

So, I think that the new statement will look something like:

EVENT (condition)
  (* Some equations solved or assignments done when the condition becomes true *)
END EVENT

We will need a way to define in which direction the boundary is crossed. For the user it can be something like this: condition in the form x == 0 means two-sided boundaries; conditions in the form x >= 0 or x <= mean one-sided boundaries: the actions associated with these boundaries will be performed not every time when x crosses 0, but only when the corresponding condition from false becomes true.

One of the use cases involves accessing the value of a variable before the boundary and after the boundary in one equation (for example, v(+) = v(-)*0.9). My suggestion is to use auxiliary variables for such models. Then the model would be written in the following way:

v, v1 IS_A speed;
t IS_A time; 
g IS_A acceleration;
DERIVATIVE OF v WITH t;
x IS_A distance;
v = der(x,t);
der(v,t) = g;
v1 = v;
EVENT (x == 0)
   v := v1*0.9;
END EVENT;

Implementation

ASCEND is efficient with large models, and this work shouldn't break this characteristic. So we need the code to be fast and efficient.

ASCEND already has some functionality for variable structure systems, so some code should be reused. As in WHENs, the relations will be compiled regardless of which values the logical variables have at any given moment. One more boolean child will be added to relations: it will indicate if the relation is a boundary equation or an ordinary one. Boundary equations should also store pointers to the boundaries which they are associated list. Still need to think what to do with assignments. Or may be we can use constant assignments so as not to break the rule that assignments are used only in methods.

In IDA we can use the code for boundary crossing as the basis. But instead of just reconfiguring the system and reinitializing IDA we should detect if there are any events associated with the boundary and, if there are, perform the needed actions.

GSOC2012

Project description

The goal of the project is to add new syntax for derivatives which will improve the capabilities of ASCEND in dynamic modelling and increase usability. Currently defining derivatives in ASCEND is not enough intuitive and convenient: for the user the process consists of, firstly, defining usual variables by using the IS_A statement, and then stating that one variable is the derivative of the other. With the new syntax the derivative of x in respect to t would look like der(x,t) and by using this expression new variables would be created automatically. Corresponding changes to the solver API would need to be made.

Project plan

  • Parser changes
    • Generate names of the type and of the new variable
    • Create a unique type description for a derivative
    • Add all needed checks
    • Append a new IS_A statement to the model’s statements.
  • Changes in the instantiator
    • Add DerInfo to RealAtomInstance.
    • Process the ‘der’ expression.
    • Do all checks which can’t be done when creating the type description.
  • Changes in the problem analysis function
    • Change the analysis function so that it would use the new relationship between variables and their derivatives.
  • Changes in the solvers.
    • Change IDA so that it would use the results of the work of analysis function fully.
    • Write the analysis function for LSODE using the results of the work of analysis function.
    • Write the analysis function for DOPRI5 using the results of the work of analysis function.

Todos for the near future

May, 9.

  1. Create a type description for a derivative atom instance. I have already started writing some code in my working copy a few days ago. I have written a function CreateDerivAtomTypeDef which is rather similar to CreateAtomTypeDef, but a derivative atom type description is always real (the types of the state and independent variables are checked before calling CreateDerivAtomTypeDef), refines solver_var and keeps pointers to the type descriptions of state and independent variables. The dimension is defined by using DiffDimensions(dimension of the state variable, dimension of the independent variable. CreateDerivAtomTypeDef is called only if the corresponding type with the name that is generated for it is not found in the type library, so the types will be unique. I think that a function which finds a type by the value of the string which a pointer points to, not by a pointer itself, will be needed for finding this type and the type solver_var. The childlist of the derivative type can’t be formed in the usual way, because this type is created when creating the model type is still in process and the derivative atom’s children would mess up with the model’s children. I suppose that the children can be just copied from solver_var – now I can’t think of any reasons why this may be not the right way to form the childlist.
  2. Make der(der(x,t),t) expressions possible. May be it should have been done earlier, but doing it now is also ok, I hope. The parser rule will be changed to der(expr,varlist). Expr can be either of type e_var or e_der. The function DoDer in typedef.c will check the type of expr, and if it is e_der, it will call DoDer for this expr and then replace the name of the state variable with the generated name.
  3. Change the independent statement. Add a special flag ‘independent’ to the solver_var.
  4. Add DerInfo to RealAtomInstance. I will need to think about the implementation.

June, 10

  1. Finish changing the functions from typedef.c so that they would use the new names.
  2. Create arrays of derivatives. My idea is that such variables as der(x[i],t[j]) can be instantiated as der(x,t)[i][j]. So that to make everything clearer for the user the information about which of the arguments is an array could be added to the string form of the name, for example: der(x[set],t[set]). So it would be easily distinguished from der(x[i][j],t), and so on. There shouldn’t be any problems with setting the DerInfo because the arguments of the derivatives are also stored in the varlist in the NameUnion. I will need to draw attention to the case when some of the derivative arguments are created inside loops.
  3. Test everything that has been done on this stage.
  4. Change the instantiator in order to process new names.

July, 13

  1. Change integrator_find_indep_var: with the use of the diffvars code it could find the independent variable in a much simpler way.
  2. Change integrator_analyse_ode so that it would use the diffvars structure. After it is done the new derivatives should work with all the currently available in ASCEND ODE/DAE solvers.
  3. Think about some more possible refinements of problem analysis?

Weekly reports

May, 21 - May, 27.

  • I had some problems with my laptop, so I installed Linux Ubuntu and everything that is needed for ASCEND on my desktop.
  • Created the derivative type description.
  • Generated an IS_A statement.
  • Nested der expressions are now parsed successfully.

May, 28 - June, 3.

  • Added some checks for ders in the methods section. Only those ders are accepted which are also used in the declarative section.
  • Added a case for expressions of type e_der to EvaluateExpr and to some other functions, mostly those which deal with relations. Now models with ders are instantiated successfully.
  • Modifyed DoDer (the function where the name and the IS_A statement are generated) so as to support derivatives of array elements. Now if a model contains such a derivative, an array of derivatives is created having the same cardinality as the array containing the state variable does.

June, 4 - June, 10.

  • Thought about derivatives in those statements where names (not expressions) are required. Discussed with mentors how to modify the grammar rule so as to make it unambiguous and to cover all possible cases. Also discussed the naming of the derivatives.
  • Added the new element to NameUnion which stores information about the arguments of the derivative. Updated name.c, name.h, nameio.c, nameio.h.
  • Started making changes to the functions from typedef.c which process derivatives using the new names.

June, 11 - June, 17.

  • Finished changing functions from typedef.c using the new names.
  • Implemented arrays of derivatives.
  • Added a test for the parser changes.
  • Fixed some bugs.

June, 18 - June, 24.

  • Improved derivatives of arrays so that they would contain only those variables which are mentioned in the model. So, if derivatives of some array elements are not used in the model, they wouldn't be created.
  • Epic fail! After a discussion with my mentors we decided that the derivatives should be declared explicitly.
  • Created a new statement: DERIVATIVE fvarlist WITH fname (WITH fname is optional). For each variable from fvarlist a derivative name is created, added to the new varlist, and an IS_A statement for this new varlist is created. After the derivatives are declared they can be accessed by names like der(x) and der(x,t).
  • Added a function to typedef.c which generates the type for the IS_A statement. Modified the function which previously generated the der variables: now the only thing it does is extending the list of derivative arguments if it comes across a derivative with a single argument and generating the symchar.
  • Added a function which finds the type for nested derivatives (or generates one).

June, 25 - Jule, 1.

  • Created a compound statement ISDER which contains a list of IS_A statements. It is needed so as to allow variables of different types in the varlist.
  • Added some new test models, updated the tests.
  • Created struct DerInfo which stores pointers to state and independent variables (for derivatives), and to derivatives (for state and independent variables). All RealAtomInstances contain this struct. If a variable is not a derivative, state or independent variable, its DerInfo is NULL.
  • Wrote functions which create DerInfo and get some information from it.
  • Added a search to MakeInstance. If the variable that is going to be instantiated is a derivative, try to find a derivative of the same state variable with respect to the same independent variable. If found, the instance is not created and the instance which was found is linked to parent. If not found, instantiate the variable and set DerInfo.
  • Implemented merges of independent variables. Still need to work on one case (which requires merging of derivatives).

July, 2 - July, 8.

  • Implemented merges of state variables and derivatives.
  • Changed FindInstances. If a model contains merges of state or independent variables, the derivative symchar may be different from the one which is added to the childlist. So, a derivative is now found by the derinfo of its state and independent variables, not by its symchar.
  • In typedef.c replaced errors when a derivative can't be found in the childlist with warnings.
  • Changed the search for a derivative which is done before creating the derivative instance. The new search makes it possible to avoid some problems with merges of arrays.
  • Fixed some problems.

July, 9 - July, 15.

  • Removed some unneeded code. Added more comments.
  • Fixed some problems in the compiler. Added more tests for the new syntax.
  • Added user's documentation for the der syntax.
  • Created new files k12_analyze.c and k12_diffvars.c with functions which build the diffvars lists and the derivative sequences using the new pointers between instances.

July, 16 - July, 22.

  • Changed functions from integrator.c: now they use the diffvars structure.
  • Added new models. Rewrote some existing models with the new syntax.
  • Added more tests.
  • Added switching between der and ode_id syntax to the GUI.
  • Converted the documentation file into LyX.

July, 23 - July, 29.

  • Added a test for DOPRI5.
  • Added more models using the new syntax.
  • Fixed a few bugs in the integrators.
  • Read the documentation of SWIG and some books about C++ and Python.

July, 30 - August, 05.

  • Added a 'der' function which can be called from the Python layer to access derivatives from Python scripts given the object wrappers on the derivative arguments.
  • Fixed some bugs in the compiler.
  • Wrote one more test model for the compiler.
  • Added some error messages.

August, 06 - August, 12.

  • Added a button for viewing DerInfo to the var properties window.
  • Worked on updating the derivative types if the derivative arguments are refined or merged.
  • Updated the tests: added some new ones and changed the relations so that they are now dimensionally correct.
  • Fixed a few minor problems in the compiler.
  • Added more comments.

Final report

Summing up the result of the work, the following goals have been achieved during this summer:

New syntax

New derivative syntax has been implemented:

DERIVATIVE OF x,y WITH t;

or

DERIVATIVE OF x,y;
INDEPENDENT t;

As a result of the above statements new variables are created: der(x,t) and der(y,t) which automatically become derivatives of x and y correspondingly. The use of new syntax is described in the documentation in a more detailed way: ksenija:doc/howto-dersyntax.lyx

There are just a few changes which will need to be made so as to support partial and higher-order derivatives. I will describe them because it may be useful for those who may want to implement partial and higher-order derivatives in the future.

Higher-order derivatives are now fully implemented on the compiler level as nested derivatives. For example, der(der(x,t),t) stores pointers to der(x,t) and t and is a second order derivative of x.

This summer I didn't care about syntax like der(x,t,t). My thought is that in this case a list of independent variable instances in DerInfo could be replaced with a list of lists of instances, where the length of each list is equal to the order of the derivative.

The implementation of the new syntax contains some first steps towards supporting partial derivatives on the compiler level. Arrays of derivatives will need to be created for derivatives with respect to arrays. The same is already done for arrays of state variables, so actually there will be nothing new. Functions DoDerIsa and DoDer from typedef.c and ExecuteISA and MakeInstance from instantiate.c will need to be changed.

Pointers between instances

All RealAtomInstances have now a pointer to a struct DerInfo. It stores gl_lists of state, independent variables, derivatives with respect to the instance and derivatives of the instance.

When some instances are merged or refined, the derivatives connected with them are also merged or refined. There can never be two derivative instances for the same state and independent variable, and the type of the derivative always matches the types of the arguments. So, the user doesn't have to care about this.

Analysis

There are now two problem analysis functions: the old and the new one. The new function builds the diffvars structure using the pointers between instances. Everything that is needed for the new analyze function is contained in files k12_analyze and k12_diffvars (ksenija:ascend/system/k12_analyze.h and ksenija:ascend/system/k12_diffvars.h). The global variable g_use_dersyntax defined in ksenija:ascend/system/system.h can be used to choose between the two functions.

The integrators now use the diffvars structure. The actions needed for building the derivative sequences are no longer repeated.

UI changes

Switching between the old and the new derivative syntax can be done using Tools - Use the der() syntax button in the PyGTK GUI. The needed option should be chosen before pressing Solve or Integrate. Otherwise the diffvars structure won't be built and user will have to reload the model.

The 'Use the der() syntax' switch in the PyGTK GUI

The DerInfo can be viewed by pressing the DerInfo button in the variable properties window.

The DerInfo window

Derivatives can be accessed from Python scripts by using the 'der' function, e.g.:

x = M.x
y = M.y
print "Der(x,y) = ",float(der(x,y))

where M is a simulation Python object (see ksenija:models/ksenija/example.py)

Ideas for the new syntax for models with hysteresis

The difficulty which I came across when trying to write some models with hysteresis was that I couldn't set the state of the system, on which its behavior depends. It should change with time, but the WHEN statement is not suited for this, I think.

I have three ideas about the new syntax. I will illustrate the use of the proposed statements with a simple house heating model suggested by John Pye. Indoor temperature 'set' to ~20 °C, and then thermal losses causing heat to leak out according to a thermal resistance equation. Then, a heater turning on when temperature is below 18 °C, and off when above 20 °C.

The first idea is to make a statement that makes it possible to change the value of a variable or to switch the status of a relation (active or inactive) when some conditions are satisfied. For example, in the house heating model it may be setting the initial temperature T0 to current temperature T when T reaches 20 °C, and then T0 will stay the same until T drops to 18. Whether the heater is turned on or off will depend on T0. Or it may be possible to change the state of the heater and the expression for T directly using the same statement without changing T0. I think that such syntax can help to describe any system with hysteresis and it is intuitive. However, this idea may have some disadvantages.

The second idea is to save the value of the last extremum (or root, there is not much difference) of some function. So T0 will be the last extremum of T, and the equation for T will change when T is below 18 or above 20, and when T is between 18 and 20 the equation will depend on T0. Though I suppose that the first variant is better because it is more general.

The third possibility is an operator that returns the value of a variable after the last boundary crossing. If it differs from the current value, it is also considered as a boundary crossed. This operator can be used to form a condition on which the current value of T0 will depend. For example, when T0 after the last boundary crossing is below 18 and T is above 20, T0 = T is used. When T0 after the last boundary crossing is above 20 and T is between 18 and 20, T0 is equal to T0 after the last boundary crossing, and so on.

Response from John

I think that ASCEND does have the notion of state -- CONDITIONAL statements, combined with the logical variables and relations, do provide a way to infer the state of a model from the current values of the conditional variables.

What seems to be lacking currently is the ability to hold any additional state information anywhere else, such as with 'sticky' (or 'memory') logical variables that retain their value from previous steps, and are only triggered to change when events happen.

Your first idea, I think, corresponds to what is already possible with WHEN. That function uses the values of logical variables to turn relations on and off. The syntax is fairly horrible, but the semantics are there (and I would like to try to fix the syntax one day).

Your second and third ideas essentially relate to adding some form of 'memory' function in ASCEND. I think that the ideas you put forward would work with the thermostat use-case, but are possibly not general enough to warrant being implemented into the language.

One possibly-general approach that to these problems that we previously put forward was the idea an 'EVENT' statement that was triggered by some kind of boundary-crossing, that causes a METHOD to be run. In that method, we could potentially do all sorts of different things, such as reversing the velocity of a ball when it bounces, etc. We could also switch the value of boolean state variables, such as "heating_on := TRUE".

One thing that is lacking in that approach, however, is the ability to write boundary equations, eg "fuel_rate(+) = fuel_rate(-) + 5 {g/s}". In other words, it would be great to be able to access the values of the variables *before* the boundary, and use those to write equations that allow us to set the state *after* the boundary.

The idea of adding boundary equations makes the whole thing much harder. In effect, one would need to launch a mini QRSlv or similar to evaluate and solve all of the boundary equations. The "IDACalcIC" solver does this task currently, but is quite limited on the kinds of initial conditions that it supports. A more flexible set would require us to write our own initial conditions code.

Here are some use-cases for this stuff think I think are good to ponder -- basically the whole topic of 'event handling' and 'hybrid simulation'.

  • thermostat controller, as previously discussed
  • a bouncing ball that *instantly* reverses its velocity (or v(+) = -v(-) * 0.9, perhaps) when hitting the ground. Note that Leon's approach uses a springy floor, instead of the instant velocity reversal approach. Sometimes we don't want to have to add this extra physics to our simulation.
  • simulation of a logic circuit, with flip-flops, and gates, and delays -- discrete state simulation
  • a flow rate controller that increments flow in fixed steps when certain thresholds are passed.
  • a tank becoming full and starting to overflow
  • a vessel with an inlet in the side, and an outlet protruding into the pipe from above; if the level is above the outlet, liquid comes out; if the level is below, gas (or nothing) comes out. 'sliding mode' is when a system like this gets stuck on the boundary or oscilates rapidly across it. how do we deal with that?

Note also that in our current ASCEND, we have all boundaries being explicitly stated, through CONDITIONAL statements. In future, we would like at least some boundaries to be automatically created, eg when you write "y = abs(x-5) + 1", you would like your solver to add a boundary at "x - 5 = 0", so that you ensure an accurate solution as the solution passes by the cusp.

-- Jpye 05:05, 10 February 2012 (EST)

I may have explained my idea not clear enough. I didn't mean that the statement which I described first should be similar to WHEN, the difference is that after the status is switched it remains the same even if the condition becomes false. So it is also some sort of 'remembering' previous states. Though the EVENT statement seems to be more general and to give the user more opportunities.

As far as I understand, writing boundary equations will require one more new statement, won't it?

So is the variant you described already accepted and can I start working on the implementation?

Today I have looked through ida.c and some other files one more time in order to learn more about how initial conditions are calculated.

--Ksenija 00:21, 11 February 2012 (EST)