User:Ksenija

Ksenija Bestuzheva
Development branch: ksenija:

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

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

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 should'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 cardinatily 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 unambigious 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.

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

• Analisys

There are now two problem analisys 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)