ASCEND - User contributions [en]

User:RaymondJoseph

2015-09-08T01:37:22Z

Ksenija: A bit of formatting

Hi There!

Thanks for stopping by.

I'm Harry Raymond Joseph, a master's student in EECS at IIT Madras, India. This page is intended to serve as a journal to keep track of my development work at Ascend as part of GSoC 2015. Updates and some relevant aspects of the GSoC proposal I submitted to Ascend are made available below.

==Final Report on the Project==

===About the Project===

----

The project was to better organize the entire ida solver module in Ascend. Specific sub-tasks included:

a) Writing a separate function to pre-process the integrator. Before the project this was packaged as part of the ida.c file.
b) Work on ida.c and make it a lot more readable and understandable.
c) Handle the events syntax, with particular requirements and constraints on time-steps after roots.
d) An overall goal was to make the code handling the solvers/ida moule compact.

===The Hows===

----
Created process_integrator.c

1. Initially, as already mentioned preparing the integrator and allocating memory to it was actually part of the ida.c module. This made it slightly more complicated and difficult to understand the main module. The separate prepare_integrator that was created has the following functionality:

2. Allocate ida_mem : This is done by use of IDA_MALLOC.

3. Pass the integrator configuration settings: several integration configuration settings (eg. Pass tolerance settings, type of solver, etc) need to be obtained from the user and passed as input to the IDA, these are accomplished by using a variety of functions, as listed in the IDA manula.

4. Analyse the equations - An initial analysis of the equations is completed – for instance checking solvability, etc.

5. configure IDA (set up the sufficient memory for y, y_d etc vectors)

Created Process_events_and_whens.c

1. Process_events_and_whens.c is actually the main module that handles all the event processing and boundary checking.

2. It returns a flag peaw (process_events_and_whens) – to the ida.c main function. Based on the value of the flag, the integrator makes several choices on the size of time-steps and direction of roots, etc.

3. In this module, the lrslv solver is called to check for roots.

4. Another important role of this function is that it returns root directions, so that if a boundary is crossed from one side, lrslv need not check that side for boundary crossings again. This is done by passing the array pointer by reference to the Process_events_and_whens.c function.

Made changes to ida.c

1. The complete configuration of the ida.c was changed. The following is a brief outline of how the new ida.c is organized:

2. From within ida.c, the solver calls first prepare_integrator.c to make necessary memory allocations and to load integrator preferences.

3. Once this has been done, the main solver routine calls process_events_and_whnens to solve initial conditions as well as to check for crossed boundaries.

4. An important point to note is that this solving of initial condition happens only once within the loop – this is highly efficient, as this enabled combining together two functions – solve_initial_conditions_and_whens and process_events_and_whens.c into one single file.

5. This was accomplished by means of the flag ‘first_run’ which is initialized to TRUE, but is set to FALSE as soon as the process_events_and_whens function is called for the first time.

6. Now, during the subsequent calls, process_events_and_whens returns ‘peaw’, which is a flag, to ascertain the status if roots were found, and smaller integration time-steps need to be taken by the main solver routine calling process_events_and_whens.

7. Peaw is 143 if roots have been found, if roots have been found it also reconfigures ida, before proceeding to next integration SMALLER timestep.

8. The two cases (normal timestep and short after_root timestep) are handled in the if statements on peaw.

9. subpeaw, is for auxiliary shorter timesteps, if there are further roots closer to the roots first found (indicated by peaw==143)

10. Once peaw is set to 143 (root found), shorter timesteps are being taken. So integration now happens from t0 to t0+tmin instead of from t0 to tout.

11. Now further shorter timesteps need to be taken, to check for changes - at this point we also take into account rootdir.

12. Say, after the first root, we increment time by tmin, and we find another root - satisfying rootdir conditions, so within the loop, the integration now changes from t0+tmin*(auxcount-1) to t0+tmin*(auxcount).

13. Once all the timesteps have been completed, the integration is complete, and the solution is returned.

===To dos in the Project===

----

1. More testing is required to see if time-steps and changing happen properly;

2. Process_events_and_whens needs t be properly debugged.

3. At least the main ida.c and prepare_integrator seem to be working fine, but it might be better to add more comments within the code and perhaps search for places to make it more efficient.

==Timeline==

'''Week 13''' (17th July to 21st August)

-Fixed all errors Ksenija pointed out on prepare_integrator earlier.

-Finishing a final report for end-term review;

-Successful in compiling working code with new main function and new prepare_integrator

-Succeesful in building new process_events_and_whens, although with the older ida.c

-Loaded some models but a syntax error was returned issue to be investigated.

'''Week 12''' (10th August to 16th August)

-Completed debugging prepare_integrator.

-Several errors found when I started debugging process_events_and_whens with the new ida.c

-Errors ranged from complex linking and type incompatibility errors to simple and trivial but difficult t fix beacket_unbalanced errors.

-Also was a good time to revisit process_events_and_whens and see if the logic underlying the code was appropriate

'''Week 11''' (3rd August to 9th August)

-Received detailed feedback from Ksenija on particularly the process_events_and_whens function. In particular
Ksenija pointed out the following issues:
"1. At the beginning you initialize the boundary state variables. But this function will do that *every* time
you call it regardless of whether it is the first call or not. I suggest putting it into an 'if (first_run || all_bnds_set == 0)'. And I think that this
boundary evaluation can be moved to a separate function - this will just be more readable.

2. setup_lrslv is called each time when process_events_and_whens
is called, while it should be called only at the first run.

3. The way you check if a boundary has been crossed is incorrect. IDA_ROOT_RETURN is a constant,
the actual status is returned by IDASolve and it should be compared to IDA_ROOT_RETURN.

4. If
no boundary is crossed (so we are making a small step after a root), you call ida_cross_boundary and then check that it does not change the structure o
f a model. But ida_log_solve can change some logical variables, and some events might be triggered, but most likely only false events, which should not really
be called. So you don't need to call the whole ida_cross_boundary function, but only the logical solver.

5. If a root has been found, you don't need to evaluate
the boundary states here."

-Looked into each of these issues, but was unable to attend to them then, since I was still debugging prepare_intgerator

-There was also an
issue with linking new files on scons. The error detailed an unusual filename. This was because of a string concatenation mistake in the SConscripts file.

-Further
re-analysed if steps being done in prepare_integrator are right;

'''Week 10''' (27th July to 2nd August)

-Fixed previous errors but got stuck at an 'unrecognized symbol' error. The error seems to be very common, but many of the fixes suggested didn't seem to work. Have requested Ksenija for help on this error.

-Fixed many corrections suggested by Ksenija regarding process_events_..

-Some errors fixed include: 'load_rellist' not defined - on running the main solver. There was a problem in linking the files.

'''Week 9''' (20th July to 26th July)
-Faced ida registration problem - Ascend GUI returns the following erroor - 'cannot register ida solver'. Solved the problem - it was due to a wrongly commented command in the prepare_integrator.

-Spoke at length with Ksenija, and shared most of the debugging details and progress by sharing screens on Skype.

-Ksenija suggested I use plots to make sure that everything is alright. Plotting in Ascend requires matplotlib, faced problems getting that sorted out.

-Continued making changes to ida.c, Ksenija pointed out a flaw in testing - was not using ida to solve the models.

-Ran tests integrtor_idaevent and integrator_ida both returned errors so far, and I'm rectifying the same.

'''Week 8''' (13th to 19th July)
-Made further changes to process_event_and_whens.

-Tested both prepare_integrator as well as process_events_and_whens by calling these from the previous versions of the ida.c.

-For this, I commented out several lines of code in both ida.c as well as idaboundary.c to see if it worked well - the returned values for models in the Ksenija branch were correct.

-Continued running tests and mostly making changes to the main solver file.

'''Week 7''' (6th to 12th July)
-Continued debugging and running all tests to fully complete the integration of the codes I wrote into the Ascend framework.

-Faced difficulties debugging. The problem was mainly due to the en bloc approach followed. Instead of debugging small pieces of code, like Ksenija advised, I ended up integrating massive code chunks (eg: ida.c) in one shot. Better code architecture planning could have helped minimize debugging trouble.

-Completed the integration of prepare_integrator.

-Completed the integration of process_events_and_whens.

-Both these functions were tested, seem to be working fine.

-Calls to these from the main function (in the ida.c file) seem to be problematic.

'''Week 6''' (29th June to 5th July)
-Changed the structure of process_events_and_whens and eliminated solve_initial conditions_and_whens.

-process_events_and_whens will now be able to process initial conditions as well. The changed outline is as below:

solve the logrels
if no values changed,
return
if event cases were triggered (ie a positive edge):
write_timestep_data
write_event
if the event has 'USE' statement(s),
run the METHOD,
run QRSlv,
run the _end method, return
write time_step data
if when cases were triggered,
update the solver rel lists
if structure has changed (how to tell?),
reanalyse lists and reconfigure IDA
do the above if first_run == false
else, its:
solve the logrels
if no values changed,
return
if when cases were triggered,
update the solver rel lists
if structure has changed
reanalyse lists and reconfigure IDA

'''Week 5''' (22nd to 28th June)
-Completed implementation of ida_main_solve.
The outline of ida_main_solve is as below:

Boundary flags initialisation: all the condition equations are evaluated
and the values of corresponding boundary flags are set.
Boolean variables initialisation: the system is solved using the logical
equations solver in order to get the correct values of Boolean variables.
Using the logical values found at step 2 the configuration of the system is chosen.
Solver initialisation: setting solver parameters, allocating memory, et cetera.
Loop through all time steps:
If a root has not been found at the previous iteration then then reset all
the root direction settings. It means that IDA will be searching for
roots of condition equations in both directions.
If a root has been found and we are far from desired output then make a very
small time step and solve the system with IDA. Else make a usual time
step which is set by the user.
If we have made a small time step then call a logical solver in order to update
the values of logical variables.
Check if a root has been found at this step. If yes:
Call a function for processing this boundary.
If there have been any changes in the system, then reinitialise the solver.
Set the direction in which the roots of condition equations will be detected:
if a boundary has been crossed, then at the next step (which is an
auxiliary small step) it will not be crossed in the same direction.
If we are still far from desired output then repeat steps.

So basically we have two modes here: the first is when we make an ordinary time step, check all roots and if a boundary is found then call the function for processing it. The second mode is entered.

'''Week 4''' (15th to 21st June)
-Completed implementation of process_events_and_whens.
The outline for process_events_and_whens is as below:

process_events_and_whens would then
solve the boundary relations
if no values changed,
return
if event cases were triggered (ie a positive edge):
write_timestep_data
write_event
if the event has 'USE' statement(s),
run the METHOD,
run QRSlv,
run the _end method, return
write time_step data
if when cases were triggered,
update the solver rel lists
if structure has changed (how to tell?),
reanalyse lists and reconfigure IDA

-Started implementation of ida_main_solve, the main ida solver function.

'''Week 3''' (8th to 14th June)
- Completed implementing prepare_integrator.

- Completed the implementation of solve_initial_conditions_and_whens.c, that's when I got stuck a little.

- Considered changes to the implementation plan below, to try and make it more efficient.

To do: Need to organize header files in a better way and make the overall code in the solvers/ida folder more compact.

'''Week 2''' (1st to 7th June)

- Took more time to read up the IDA Sundials manual, especially the architecture of the IDA solver module.

- Received the following implementation outline from Ksenija and John:

integrator_ida_solve
prepare_integrator
solve_initial_conditions_and_whens
write_timestep_data
while(tsteps_remaining)
integrate_to_next_step
if(root triggered)
process_events_and_whens
advance timestep if we're really close already
else:
write_timestep_data
advance to next timestep
process_events_and_whens would then
solve the logrels
if no values changed, return
if event cases were triggered (ie a positive edge):
write_timestep_data
write_event
if the event has 'USE' statement(s),
run the METHOD,
run QRSlv,
run the _end method, return
write time_step data
if when cases were triggered,
update the solver rel lists
if structure has changed (how to tell?)
reanalyse lists and reconfigure IDA

solve_initial_conditions_and_whens would be

once the continuous part of the problem is solved, we need to evaluate the bnds, and then run LRSlv.
If that results in triggering any changes to WHEN cases, then we should
reanalyse the equations
reconfigure IDA
go back and solver the continuous problem again (loop)
If there are too many iterations, then we should return an error to the user.
we don't want to trigger any EVENTs during this process; only WHENs.

prepare_integrator

allocate ida_mem etc
pass the integrator configuration settings
analyse the equations
configure IDA (set up the sufficient memory for y, y_d etc vectors)

- looking to complete before 16th June

'''Week 1''' (25th to 31st May)

- Setup CUnit and ran some tests

- Spent time understanding IDA issues and implementation aspects

- Looked at ida.c and idaboundary.c in the ksenija2 branch of the svn repository.

- Worked on enhancing this code based on an implementation outline from John.

'''Community bonding period (27th April to 25th May)'''

- Installed Linux, and all other dependencies to prepare for coding.

- Built Ascend on my 64-bit 15.04 Ubuntu system. Log [[File:attached]].

- Reading the code in the Ksenija2 branch of the SVN repository.

==GSoC 2015 Proposal==

Organization: ASCEND

===Abstract:=== The IDA solver requires more refinements both in the numerical as well as user-facing aspects is desirable. One important issue is to do with EVENT syntax - which helps to model whole range of systems, in which parametric constants themselves change values with triggered events, or even more sophisticated systems in which constraint equations are added or deleted when a given event is triggered. This project attempts to strengthen the functionalities of the IDA solver.
(* I'm using the proposal format, obtained from Carlos' GSoC 2012 proposal *)

===About:===
Master's student at TUM Germany/IITM India. Specializing in the numerical aspects in EECS, such as machine learning and probabilistic methods. Research enthusiast, few publications. Amateur poet and writer.

Email:
raymond.harry@tum.de/raymond.joseph.7@gmail.com

Project Topic
Numerical Algorithmic Enhancement: Dynamic Modelling in ASCEND

Basic Information

URL link to your CV: https://drive.google.com/file/d/0B04lAfTBvrd0a1dWc3hPYnZQTkU/view?usp=sharing

URL link to your personal blog/homepage: http://raymondjoseph7.wix.com/ctkotww

Research page: http://dl.acm.org/author_page.cfm?id=92458619092

===Goals of the project===

The main goal (Deliverable 1, D1) is to fully integrate IDA into ASCEND by providing and refining EVENT support syntax elements, especially ensuring that the integrator restarts efficiently after events are triggered with modified constraint equations.
Secondary goals are:

- Improve IDA numerical solver efficiency. (D2)

- Refine existing code in the \ksenija2 branch. (D3)

- Possible GUI improvements in the use of the IDA solver. (D4)

- Event syntax for periodic events. (D5)

- Implement model of the simple solar water heater system as a test case. (D6)

- Make the documentation of migration and faq for users with applications. (D7)

===Schedule===

Week 1: Read up and understand IDA integration issues. Great sources and discussions have been documented on the integration issues.

Week 2: Scrub existing \ksenija2 code.

Week 3, 4 & 5: Implement advanced EVENT syntax, to complete integration of IDA into ASCEND. (D1)

Week 6: Explore possible syntax architectural changes for the IDA solver problems. (D3)

Week 7: Read latest cutting-edge research papers and see if new methods are available to increase solver efficiency. (Hint: Yes, more optimal methods are available, that perform better than the IDA solver adopted as part of the SUNDIALS project) (D2)

Week 8: Implement solver improvements to make the IDA solver state-of-the-art. (D2)

Week 9: GUI improvements. (D4)

Week 10: Periodic events support syntax. (D5)

Week 11: Buffer week for finalizing code, and running tests. Fix any bugs discovered.

Week 12: Implement test case: solar water heater system. (D6)

Week 13: Documentation and FAQs. (D7)

===Coding experience===
Programming languages you have learnt, and how many lines of code, approximately, you have written in each.
C++ >1500 LoC (Lots of numerical methods experience. Monte-carlo methods, Newton-Raphson method, ODE simulation, etc.)
MatLab >500 LoC (Lots of experience with advanced numerical methods. Example: water-filling, convex optimization using CVX, modified gradient-descent methods, linear algebra systems in communication, etc.)
Java >400 LoC
Python >1000 LoC (* Though most of it is classified. I worked as a data analytics programmer at a French startup, Tilkee *)

I also have lex/yacc experience. Used these to create my own 'business artifact modelling language' called BASQL. Based on SQL, the language was relatively successful, a paper on this was published at the IEEE ICES conference jointly with Youakim Badr. The language achieved great success and further research on this is being carried out by Youakim Badr at LIRIS-CNRS, INSA Lyon, France.

Some of these coding samples is provided in the URL.

MatLab: https://drive.google.com/folderview?id=0B04lAfTBvrd0fjhuLWFBNHprOXpsdlhDbU5XSmd3WnNfa3JqakdXRFRScmxIYkdfeTJKWGc&usp=sharing

C++: https://drive.google.com/folderview?id=0B04lAfTBvrd0TWl2T1BBSXlFNDQ&usp=sharing

Have you written any ASCEND models? Please provide copies.

I have been able to write two successful (albeit, very simple models):

1. The first one is on the 'dynamics of love affairs' paper by (Felmlee and Greenberg). The model can be found here:
https://drive.google.com/file/d/0B04lAfTBvrd0SjU2aHZTZ2lwYkE/view?usp=sharing
The snapshot of the solved model is below:
https://drive.google.com/file/d/0B04lAfTBvrd0RmRmTzh5MVQ2WmM/view?usp=sharing

2. The second model is the simple non-flow domestic solar water heater system. The model can be found here:
https://drive.google.com/file/d/0B04lAfTBvrd0SjU2aHZTZ2lwYkE/view?usp=sharing
The snapshot of the solved model is below:
https://drive.google.com/file/d/0B04lAfTBvrd0TERrcG5nbS1kckU/view?usp=sharing

[[Category:GSOC2015]]

File:Derinfo.png

2014-09-14T18:38:37Z

Ksenija: Ksenija uploaded a new version of "File:Derinfo.png"

File:Use der syntax.png

2014-09-14T18:38:22Z

Ksenija: Ksenija uploaded a new version of "File:Use der syntax.png"

File:Huygens clock.png

2014-09-14T18:36:58Z

Ksenija: Ksenija uploaded a new version of "File:Huygens clock.png"

File:Plot events3.png

2014-09-14T18:36:20Z

Ksenija: Ksenija uploaded a new version of "File:Plot events3.png"

File:Plot events2.png

2014-09-14T18:36:05Z

Ksenija: Ksenija uploaded a new version of "File:Plot events2.png"

File:Plot events1.png

2014-09-14T18:35:10Z

Ksenija: Ksenija uploaded a new version of "File:Plot events1.png"

File:Vessel with siphon.png

2014-09-14T18:34:20Z

Ksenija: Ksenija uploaded a new version of "File:Vessel with siphon.png"

File:Hysteron.png

2014-09-14T18:33:28Z

Ksenija: Ksenija uploaded a new version of "File:Hysteron.png"

Publications

2014-07-07T15:13:19Z

Ksenija:

'''Note: not all the PDF links are present in the following text. If a link is missing, please using the original [http://www.ascend4.org/ascend_bibliography.htm ASCEND bibliography] page to find the document.'''

In date order, divided according to category.

== Doctoral theses ==

V Rico-Ramirez, 1998. ''Representation, Analysis and Solution of Conditional Models in an Equation-Based Environment'', Ph.D. Thesis, Carnegie Mellon University. [http://ascend.cheme.cmu.edu/ftp/pdfThesis/victhesis.pdf PDF]

K A Abbott, 1996. ''Very Large Scale Modeling'' Ph.D. Thesis, Carnegie Mellon University, Pittsburgh, PA. [http://ascend.cheme.cmu.edu/ftp/pdfThesis/kirkthesis.pdf PDF]

B A Allan, 1997. ''A More Reusable Modeling System'', Ph.D. Thesis, Carnegie Mellon University, Pittsburgh, PA. [http://ascend.cheme.cmu.edu/ftp/pdfThesis/benthesis.pdf PDF]

B T Safrit, 1996. ''Synthesis of Azeotropic Batch Distillation Separation Systems'', Ph.D. Thesis, Carnegie Mellon University, Pittsburgh, PA. [http://ascend.cheme.cmu.edu/ftp/pdfThesis/safritthesis.pdf PDF]

M Thomas, 1996. ''Tool and Information Management in Engineering Design'', Ph.D. Thesis, Carnegie Mellon University, Pittsburgh, PA. [http://ascend.cheme.cmu.edu/ftp/pdfThesis/markthesis.pdf PDF]

R S Huss, 1995. ''Collocation Methods For Flexible Distillation Design'', Ph.D. Thesis, Carnegie Mellon University, Pittsburgh, PA. [http://ascend.cheme.cmu.edu/ftp/pdfThesis/bobthesis.pdf PDF]

J Zaher, 1995. ''Conditional Modeling'', Ph.D. Thesis, Carnegie Mellon University, Pittsburgh, PA. [http://ascend.cheme.cmu.edu/ftp/pdfThesis/joethesis.pdf PDF]

P Piela, 1989. ''ASCEND: An Object-Oriented Computer Environment for Modeling and Analysis'', Ph.D. Thesis, Carnegie Mellon University, Pittsburgh, PA. [http://ascend.cheme.cmu.edu/ftp/htbibliography/b9.html Abstract]

== Journal papers ==

A W Westerberg, 2003, ''A retrospective on design and process synthesis'', Computers & Chemical Engineering, '''28''' (4), 447-458. {{doi|10.1016/j.compchemeng.2003.09.029}}.

V Rico-Ramírez, B A Allan and A W Westerberg, 1999, ''Conditional Modeling. 1. Requirements for an Equation-Based Environment'', Ind. Eng. Chem. Res., 38 (2), pp 519–530. {{doi|10.1021/ie9800593}}

V Rico-Ramírez and A W Westerberg, 1999, ''Conditional Modeling. 2. Solving Using Complementarity and Boundary-Crossing Formulations'', Ind. Eng. Chem. Res., 38 (2), pp 531–553. {{doi|10.1021/ie9800602}}

Benjamin A Allan and Arthur W Westerberg, 1999, ''Anonymous Class in Declarative Process Modeling'', Ind. Eng. Chem. Res., '''38''' (3), pp 692–704. {{doi|10.1021/ie980297y}}

Kirk A Abbott, Benjamin A Allan, and Arthur W Westerberg, 1997, ''Global preordering for Newton equations using model hierarchy'', AIChE Journal '''43''' (12), {{doi|10.1002/aic.690431207}}

Robert S Huss and Arthur W Westerberg, 1996, ''Collocation Methods for Distillation Design. 2. Applications for Distillation'', Ind. Eng. Chem. Res., '''35''' (5), pp 1611–1623, {{doi|10.1021/ie9503508}}

B T Safrit, A W Westerberg, U Diwekar, O M Wahnschafft, 1995, ''Extending Continuous Conventional and Extractive Distillation Feasibility Insights to Batch Distillation'', Ind. Eng. Chem. Res., 34, 3257-3264. {{doi|10.1021/ie00037a012}}

B Katzenberg and P Piela, 1993, ''Work Language Analysis and the Naming Problem'', Communication of the ACM, Vol. 36, No. 6, 86-92. {{doi|10.1145/153571.163286}}

P Piela, R McKelvey, and A Westerberg, 1993, ''An Introduction to the ASCEND Modeling System: Its Language and Interactive Environment'', ; Journal of Management Information Systems ,Vol. 9, No.3, 91-121. [http://ascend.cheme.cmu.edu/ftp/htbibliography/b6.html Abstract], [http://www.jstor.org/stable/40398044 Full article]

P Piela, B Katzenberg, and R McKelvey, 1992, ''Integrating the User into Research on Engineering Design Systems ASCEND: An Object-Oriented Computer Environment for Modeling'', Research in Engineering Design, Vol. 3, 211-221. {{doi|10.1007/BF01580843}}

Oliver J Smith IV and Arthur W Westerberg, 1991, ''The optimal design of pressure swing adsorption systems'', Chemical Engineering Science, '''48'' (12), {{doi|10.1016/0009-2509(91)85001-E}}

P Piela, T Epperly, K Westerberg and A Westerberg, 1991, ''ASCEND: An Object-Oriented Computer Environment for Modeling and Analysis: The Modeling Language.'', Computers and Chemical Engineering, Vol. 15, No. 1, 53-72, . {{doi|10.1016/0098-1354(91)87006-U}}

R F Woodbury, 1990, ''Variations in Solids: A Declarative Treatment'', Computers and Graphics, Special Issue on Features and Geometric Reasoning, Vol. 14, No. 2, 173-188. {{doi|10.1016/0097-8493(90)90030-2}}

S Kuru and A W Westerberg, 1985, ''A Newton-Raphson based strategy for exploiting latency in dynamic simulation'', Computers & Chemical Engineering, Vol 9., No. 2, {{doi|10.1016/0098-1354(85)85007-9}} or [[Media:Kuru1985.pdf|PDF]]

A W Westerberg and D R Benajamin, 1985, ''Thoughts on a Future Equation-Oriented Flowsheeting System'', Computers and Chemical Engineering, Vol. 9, No. 5, 517-526. {{doi|10.1016/0098-1354(85)80026-0}}

Michael H Locke and Arthur W Westerberg, 1983, ''The ascend-II system—a flowsheeting application of a successive quadratic programming methodology'', Computers & Chemical Engineering '''7''' (5) 615-630. {{doi|10.1016/0098-1354(83)80007-6}}

A W Westerberg and S W Director, 1978, ''A modified least squares next term algorithm for solving sparse n × n sets of nonlinear equations'', Computers and Chemical Engineering, Vol. 2, pp 77-81. {{doi|10.1016/0098-1354(78)80011-8}} (this paper describes the basis of our [[QRSlv]] solver)

== Conference presentations and papers ==
A Sadrieh, P A Bahri, ''Application of Graphic Processing Unit in Model Predictive Control'', Computer Aided Chemical Engineering, Elsevier, 2011, Volume 29, 21st European Symposium on Computer Aided Process Engineering, Pages 492-496, ISSN 1570-7946, ISBN 9780444538956, {{doi|10.1016/B978-0-444-53711-9.50099-7}} or try [http://www.sciencedirect.com/science/article/pii/B9780444537119500997 this link].

J Pye, K Lovegrove and G Burgess. ''[http://solar-thermal.anu.edu.au/files/2010/10/pye-2010-solarpaces-combined-cycle.pdf Combined-cycle solarised gas turbine with steam, organic and CO2 bottoming cycles]'', Proceedings of SolarPACES 2010, Perpignan, France, Sept 2010.

J Coventry and J Pye, 2009, ''[http://solar-thermal.anu.edu.au/files/2010/02/SolarPACES2009_CoventryPye.pdf Coupling supercritical and superheated direct steam generation with thermal energy storage]''. SolarPACES 2009, Berlin, 15-18 Sept.

H G Silva and R L R Salcedo, 2006, ''Modeling and Optimization of Chemical Processes: ASCEND IV and Stochastic Optimizers'', Proceedings of [http://www.actapress.com/Content_Of_Proceeding.aspx?ProceedingID=381 Modelling and Simulation], Montreal, Canada.

Allan, B.A. and A. W. Westerberg, ''Compiling and Solving 100,000 equations on a PC in (3) minutes'' Carnegie Mellon University, Pittsburgh, PA, 1998. Presented April 1998 at INFORMS Montreal conference. [http://ascend.cheme.cmu.edu/ftp/pdfPapersRptsSlides/informs98.pdf PDF].

Bhargava, H., Krishnan, R., and Piela, P., ''Formalizing the Semantics of ASCEND'', Proceeding of the 27th Hawaii International Conference on the System Sciences, IEEE Computer Society Press, 1994.
[http://ascend.cheme.cmu.edu/ftp/htbibliography/b1.html Abstract], {{doi|http://dx.doi.org/10.1109/HICSS.1994.323312}}.

Krishnan, R., Piela, P., and Westerberg, A., ''Reusing Mathematical Models in ASCEND'', Recent Developments in Decision Support Systems, C. W. Holsapple and A. B. Whinston (eds.). Proceedings of the NATO ASI on Decision Support Systems. Springer-Verlag, in cooperation with NATO Scientific Affairs Division, 1993. [http://ascend.cheme.cmu.edu/ftp/htbibliography/b5.html Abstract], read [http://books.google.com.au/books?id=UefgA0jltcgC&pg=PA275 via Google Books].

Westerberg, A.W., Abbott, K.A., and Allan, B.A., ''Plans for ASCEND IV: Our Next Generation Equational-Based Modeling Environment'', AspenWorld 94, Boston Massachusetts, November 1994. [http://ascend.cheme.cmu.edu/ftp/htbibliography/b10.html abstract], [[Media:westerberg-1994-aspenworld.pdf|paper]], [[Media:AspenWorldASCENDSlides.pdf|slides]].

Westerberg, A., Piela, P., McKelvey, R., and Epperly, T., ''The ASCEND Modeling Environment and Its Implications'', Proceedings of the 4th International Symposium on Process Systems Engineering, I: Design, I.2.1 - I.2.12, 1991. [http://ascend.cheme.cmu.edu/ftp/htbibliography/ba.html Abstract].

Westerberg, A., Piela, P., Subrahmanian, E., Podnar, G., and Elm, W., ''A Future Computer Environment for Preliminary Design'', Proceedings of the Third International Conference on Foundations of Computer Aided Process Design, J. Siirola, I. Grossmann, and G. Stephanopoulos (eds.), 1989. [http://ascend.cheme.cmu.edu/ftp/htbibliography/bb.html Abstract]

== Technical reports ==

L Cisternas, N Luza, E Gálvez, 2007. ''[[Media:HydroSim-reporte-tecnico.pdf|HydroSim: Una Librería para Simulación, Modelación y Optimización en Hidrometalurgia]]'' (in spanish; english title: ''HydroSim: a library for simulation, modelling and optimisation in hydrometallurgy''), Technical report, Universidad de Antifagasta, Design of Product and Process Group and CICITEM. '''PDF'''

J L Perry and B A Allan, 1996. ''[http://repository.cmu.edu/compsci/2143/ Design and Use of Dynamic Modeling in ASCEND IV]'', Carnegie Mellon University, EDRC Technical Report 06-224-96. '''PDF''' ([http://ascend.cheme.cmu.edu/ftp/pdfPapersRptsSlides/AscendIVP.pdf alternative link])

P Piela and A W Westerberg, 1994. ''[http://ascend.cheme.cmu.edu/ftp/pdfPapersRptsSlides/processModeling.pdf Equation-Based Process Modeling]'', Carnegie Mellon University, EDRC Technical Report. '''PDF'''.

A W Westerberg, 1994. ''[http://ascend.cheme.cmu.edu/ftp/pdfPapersRptsSlides/ascendIntro.pdf ASCEND Modeling Language and Environment Notes]'', CAPD Short Course Notes, Carnegie Mellon University.

T Epperly, 1989. [[Media:epperlyPaperICES05-29-89.pdf|''Implementation of an ASCEND Interpreter'']], Carnegie Mellon University, EDRC Technical Report 05-29-89, '''PDF'''.

B A Allan and A W Westerberg, 1997.''[http://ascend.cheme.cmu.edu/ftp/pdfPapersRptsSlides/aiche97-201f.pdf Reusability and Scalability in Modeling] (Previously titled "Reusability in Modeling and Modeling Support Environments")'', Carnegie Mellon University, ICES Technical Report, November 1997. Pittsburgh, PA. (AIChE 1997 Meeting Paper 201f.) (see also [http://ascend.cheme.cmu.edu/ftp/pdfPapersRptsSlides/slides-201f.pdf presentation slides])

B A Allan and A W Westerberg, 1997. ''The ASCEND IV Language Syntax and Semantics'', Carnegie Mellon University, EDRC Technical Report. (this report is now edited and incorporated into our [[:Category:Documentation|documentation]]).

B A Allan, V Rico-Ramirez, M Thomas and K Tyner, 1996. ''ASCEND IV: A portable Mathematical Modeling Environment'', Carnegie Mellon University, ICES Technical Report, October. [http://ascend.cheme.cmu.edu/ftp/pdfHelp/ascend-help-BOOK-3.pdf PDF].

K Dee and A Westerberg, 1922. ''CEPHDA: Chemical Engineering Process Hierarchical Design with ASCEND'', Carnegie Mellon University, EDRC Technical Report 06-140-92 [http://ascend.cheme.cmu.edu/ftp/htbibliography/b2.html Abstract].

V Rico-Ramirez, B A Allan and A W Westerberg, 1998. ''Conditional Modeling in an Equation-Based Environment'', Carnegie Mellon University, ICES Technical Report 06-242-98. [http://ascend.cheme.cmu.edu/ftp/pdfPapersRptsSlides/tech_modeling.pdf PDF].

V Rico-Ramirez and A W Westerberg, 1998. ''Complementarity Formulation for the Representation of Algebraic Systems Containing Conditional Equations'', Carnegie Mellon University, ICES Technical Report 06-243-98. [http://ascend.cheme.cmu.edu/ftp/pdfPapersRptsSlides/tech_complementarity.pdf PDF].

K Westerberg, 1989. ''Development of Software for Solving Systems of Nonlinear Equations'', Carnegie Mellon University, EDRC Technical Report 05-36-89, Pittsburgh, PA. [http://ascend.cheme.cmu.edu/ftp/pdfPapersRptsSlides/nonlinear_rpt.pdf PDF]

K Westerberg, 1989. ''Development of Software for Solving Systems of Linear Equations'', Carnegie Mellon University, EDRC Technical Report 05-35-89. [http://ascend.cheme.cmu.edu/ftp/pdfPapersRptsSlides/linear_rpt.pdf PDF]

J Zaher, 1991. ''Developing Reusable Libraries in the ASCEND environment'', Carnegie Mellon University, EDRC Technical Report 06-108-91. [http://ascend.cheme.cmu.edu/ftp/htbibliography/bf.html PDF].

== Undergraduate theses ==

K Bestuzheva, 2014. ''Computer modelling of successive collisions of a particle with moving boundaries in an external magnetic field'' (in Russian ''Компьютерное моделирование последовательных отражений частицы от движущихся поверхностей при наличии внешнего магнитного поля''). Final thesis for the degree of a specialist, State Management University, Moscow, Russia. [[:File:Bestuzheva_final_thesis.pdf | PDF]]

R Hogan, 2012. ''Efficiency of supercritical carbon dioxide for the working fluid in the power cycle for high-temperature concentrated solar power''. BE thesis, Australian National University. Available on request.

K Sio, 2010. ''Cost Optimisation of DSG Thermal Storage System''. BE thesis, Australian National University. Available on request.

2013-09-20T14:18:56Z

Ksenija:

File:Plot events2.png

2013-09-20T14:18:07Z

Ksenija:

File:Plot events1.png

2013-09-20T14:17:10Z

Ksenija:

User:Ksenija

2013-09-20T13:26:25Z

Ksenija:

'''Ksenija Bestuzheva'''
 Development branch: {{srcbranchdir|ksenija|}}, {{srcbranchdir|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:

<source lang=a4c>
EVENT (x > 0)
...
END EVENT;
</source>

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:

<source lang=a4c>
a IS_A submodel1;
...
EVENT(...)
CASE TRUE:
USE a;
x0 = x;
END EVENT;
</source>

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

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

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

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

{{srcbranch|ksenija2|models/ksenija/bball_event3.a4c}} - another version of the bouncing ball model.

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

[[File:hysteron.png|thumb|none|A sketch of the cilinder]]

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

{{srcbranch|ksenija2|models/ksenija/simultaneous.a4c}} - a model with two boundaries crossed simultaneously.

{{srcbranch|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".

[[File:vessel_with_siphon.png|thumb|none|A sketch of the vessel]]

{{srcbranch|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, 02 - September, 08.'''

* Added a model of the Huygens pendulum clock: {{srcbranch|ksenija2|models/ksenija/clock.a4c}}

[[File:huygens_clock.png|thumb|none|The phase plane of the Huygens clock model]]

* 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: {{srcbranch|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.

[[File:Plot_event.png|thumb|none|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.

==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 [https://en.wikipedia.org/wiki/Delta-sigma_modulation 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:
<source lang=a4c>
EVENT (condition)
(* Some equations solved or assignments done when the condition becomes true *)
END EVENT
</source>

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:
<source lang=a4c>
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;
</source>

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

# '''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.
# '''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.
# '''Change the independent statement.''' Add a special flag ‘independent’ to the solver_var.
# '''Add DerInfo to RealAtomInstance.''' I will need to think about the implementation.

'''June, 10.'''

# '''Finish changing the functions from typedef.c''' so that they would use the new names.
# '''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.
# '''Test''' everything that has been done on this stage.
# '''Change the instantiator''' in order to process new names.

'''July, 13.'''

# '''Change integrator_find_indep_var''': with the use of the diffvars code it could find the independent variable in a much simpler way.
# '''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.
# 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:

<source lang=a4c>
DERIVATIVE OF x,y WITH t;
</source>

or

<source lang=a4c>
DERIVATIVE OF x,y;
INDEPENDENT t;
</source>

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: {{srcbranch|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 ({{srcbranch|ksenija|ascend/system/k12_analyze.h}} and {{srcbranch|ksenija|ascend/system/k12_diffvars.h}}). The global variable g_use_dersyntax defined in {{srcbranch|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.

[[File:Use_der_syntax.png|200px|thumb|center|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.

[[File:Derinfo.png|200px|thumb|center|The DerInfo window]]

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

<source lang=py>
x = M.x
y = M.y
print "Der(x,y) = ",float(der(x,y))
</source>

where <tt>M</tt> is a simulation Python object (see {{srcbranch|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.

-- [[User:Jpye|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.

--[[User:Ksenija|Ksenija]] 00:21, 11 February 2012 (EST)

[[Category:ASCEND Contributors]]
[[Category:GSOC2012]]

Ksenija:

'''Ksenija Bestuzheva'''
 Development branch: {{srcbranchdir|ksenija|}}, {{srcbranchdir|ksenija2|}}

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, control theory. I have experience with C, C++, Java programming through my studies.

==GSOC2013==

===Example models.===

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

{{srcbranch|ksenija2|models/ksenija/thermostat.a4c}} - a model of a thermostat.

{{srcbranch|ksenija2|models/ksenija/bball_event.a4c}}, {{srcbranch|ksenija2|models/ksenija/bball_event2.a4c}} - a bouncing ball. When it hits the ground, it instantly changes its speed. Some problems concerning this models still remain. First, I don't know what to do with the system not being square. Adding and subtracting y is definitely not a good idea. Second, when the ball almost rests on the floor it is hard to process all boundary crossings correctly.

{{srcbranch|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 cilinder, the position of the piston is the input u, the position of the cilinder's right end is the output x.

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

{{srcbranch|ksenija2|models/ksenija/simultaneous.a4c}} - a model with two boundaries crossed simultaneously.

{{srcbranch|ksenija2|models/ksenija/vessel_with_siphon.a4c}} - a model of auto-oscillations of fluid level in the vessel with simultaneous inflow and outflow.

'''Models which are not yet written:'''

* A vessel with a sliding mode.
* Some other models?

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

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

The screenshot of the plot for the test model:

[[File:Plot_event.png]]

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

==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 [https://en.wikipedia.org/wiki/Delta-sigma_modulation 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 unambigious, flexible and clear.

====Proposed syntax.====

So, I think that the new statement will look something like:
<source lang=a4c>
EVENT (condition)
(* Some equations solved or assignments done when the condition becomes true *)
END EVENT
</source>

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:
<source lang=a4c>
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;
</source>

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

# '''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.
# '''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.
# '''Change the independent statement.''' Add a special flag ‘independent’ to the solver_var.
# '''Add DerInfo to RealAtomInstance.''' I will need to think about the implementation.

'''June, 10.'''

# '''Finish changing the functions from typedef.c''' so that they would use the new names.
# '''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.
# '''Test''' everything that has been done on this stage.
# '''Change the instantiator''' in order to process new names.

'''July, 13.'''

# '''Change integrator_find_indep_var''': with the use of the diffvars code it could find the independent variable in a much simpler way.
# '''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.
# 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:

<source lang=a4c>
DERIVATIVE OF x,y WITH t;
</source>

or

<source lang=a4c>
DERIVATIVE OF x,y;
INDEPENDENT t;
</source>

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: {{srcbranch|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 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 ({{srcbranch|ksenija|ascend/system/k12_analyze.h}} and {{srcbranch|ksenija|ascend/system/k12_diffvars.h}}). The global variable g_use_dersyntax defined in {{srcbranch|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.

[[File:Use_der_syntax.png|200px|thumb|center|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.

[[File:Derinfo.png|200px|thumb|center|The DerInfo window]]

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

<source lang=py>
x = M.x
y = M.y
print "Der(x,y) = ",float(der(x,y))
</source>

where <tt>M</tt> is a simulation Python object (see {{srcbranch|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.

-- [[User:Jpye|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.

--[[User:Ksenija|Ksenija]] 00:21, 11 February 2012 (EST)

[[Category:ASCEND Contributors]]
[[Category:GSOC2012]]

2013-07-08T07:29:28Z

Ksenija:

File:Plot.png

2013-07-08T07:21:12Z

Ksenija: uploaded a new version of "File:Plot.png"

User:Ksenija

2013-07-01T13:57:28Z

Ksenija:

'''Ksenija Bestuzheva'''
 Development branch: {{srcbranchdir|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, control theory. I have experience with C, C++, Java programming through my studies.

==GSOC2013==

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

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

==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: http://ascend4.org/Event_handling. Some more interesting examples may concern encoding analog signals into digital signals (see https://en.wikipedia.org/wiki/Delta-sigma_modulation). 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 unambigious, flexible and clear.

====Proposed syntax.====

So, I think that the new statement will look something like:
<source lang=a4c>
EVENT (condition)
(* Some equations solved or assignments done when the condition becomes true *)
END EVENT
</source>

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:
<source lang=a4c>
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;
</source>

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

# '''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.
# '''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.
# '''Change the independent statement.''' Add a special flag ‘independent’ to the solver_var.
# '''Add DerInfo to RealAtomInstance.''' I will need to think about the implementation.

'''June, 10.'''

# '''Finish changing the functions from typedef.c''' so that they would use the new names.
# '''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.
# '''Test''' everything that has been done on this stage.
# '''Change the instantiator''' in order to process new names.

'''July, 13.'''

# '''Change integrator_find_indep_var''': with the use of the diffvars code it could find the independent variable in a much simpler way.
# '''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.
# 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:

<source lang=a4c>
DERIVATIVE OF x,y WITH t;
</source>

or

<source lang=a4c>
DERIVATIVE OF x,y;
INDEPENDENT t;
</source>

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: {{srcbranch|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 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 ({{srcbranch|ksenija|ascend/system/k12_analyze.h}} and {{srcbranch|ksenija|ascend/system/k12_diffvars.h}}). The global variable g_use_dersyntax defined in {{srcbranch|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.

[[File:Use_der_syntax.png|200px|thumb|center|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.

[[File:Derinfo.png|200px|thumb|center|The DerInfo window]]

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

<source lang=py>
x = M.x
y = M.y
print "Der(x,y) = ",float(der(x,y))
</source>

where <tt>M</tt> is a simulation Python object (see {{srcbranch|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.

-- [[User:Jpye|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.

--[[User:Ksenija|Ksenija]] 00:21, 11 February 2012 (EST)

[[Category:ASCEND Contributors]]
[[Category:GSOC2012]]

User:Ksenija

2013-06-24T06:31:59Z

Ksenija:

'''Ksenija Bestuzheva'''
 Development branch: {{srcbranchdir|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, control theory. I have experience with C, C++, Java programming through my studies.

==GSOC2013==

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

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

==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: http://ascend4.org/Event_handling. Some more interesting examples may concern encoding analog signals into digital signals (see https://en.wikipedia.org/wiki/Delta-sigma_modulation). 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 unambigious, flexible and clear.

====Proposed syntax.====

So, I think that the new statement will look something like:
<source lang=a4c>
EVENT (condition)
(* Some equations solved or assignments done when the condition becomes true *)
END EVENT
</source>

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:
<source lang=a4c>
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;
</source>

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

# '''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.
# '''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.
# '''Change the independent statement.''' Add a special flag ‘independent’ to the solver_var.
# '''Add DerInfo to RealAtomInstance.''' I will need to think about the implementation.

'''June, 10.'''

# '''Finish changing the functions from typedef.c''' so that they would use the new names.
# '''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.
# '''Test''' everything that has been done on this stage.
# '''Change the instantiator''' in order to process new names.

'''July, 13.'''

# '''Change integrator_find_indep_var''': with the use of the diffvars code it could find the independent variable in a much simpler way.
# '''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.
# 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:

<source lang=a4c>
DERIVATIVE OF x,y WITH t;
</source>

or

<source lang=a4c>
DERIVATIVE OF x,y;
INDEPENDENT t;
</source>

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: {{srcbranch|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 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 ({{srcbranch|ksenija|ascend/system/k12_analyze.h}} and {{srcbranch|ksenija|ascend/system/k12_diffvars.h}}). The global variable g_use_dersyntax defined in {{srcbranch|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.

[[File:Use_der_syntax.png|200px|thumb|center|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.

[[File:Derinfo.png|200px|thumb|center|The DerInfo window]]

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

<source lang=py>
x = M.x
y = M.y
print "Der(x,y) = ",float(der(x,y))
</source>

where <tt>M</tt> is a simulation Python object (see {{srcbranch|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.

-- [[User:Jpye|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.

--[[User:Ksenija|Ksenija]] 00:21, 11 February 2012 (EST)

[[Category:ASCEND Contributors]]
[[Category:GSOC2012]]

User:Ksenija

2013-06-17T13:57:18Z

Ksenija:

'''Ksenija Bestuzheva'''
 Development branch: {{srcbranchdir|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, control theory. I have experience with C, C++, Java programming through my studies.

==GSOC2013==

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

==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: http://ascend4.org/Event_handling. Some more interesting examples may concern encoding analog signals into digital signals (see https://en.wikipedia.org/wiki/Delta-sigma_modulation). 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 unambigious, flexible and clear.

====Proposed syntax.====

So, I think that the new statement will look something like:
<source lang=a4c>
EVENT (condition)
(* Some equations solved or assignments done when the condition becomes true *)
END EVENT
</source>

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:
<source lang=a4c>
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;
</source>

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

# '''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.
# '''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.
# '''Change the independent statement.''' Add a special flag ‘independent’ to the solver_var.
# '''Add DerInfo to RealAtomInstance.''' I will need to think about the implementation.

'''June, 10.'''

# '''Finish changing the functions from typedef.c''' so that they would use the new names.
# '''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.
# '''Test''' everything that has been done on this stage.
# '''Change the instantiator''' in order to process new names.

'''July, 13.'''

# '''Change integrator_find_indep_var''': with the use of the diffvars code it could find the independent variable in a much simpler way.
# '''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.
# 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:

<source lang=a4c>
DERIVATIVE OF x,y WITH t;
</source>

or

<source lang=a4c>
DERIVATIVE OF x,y;
INDEPENDENT t;
</source>

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: {{srcbranch|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 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 ({{srcbranch|ksenija|ascend/system/k12_analyze.h}} and {{srcbranch|ksenija|ascend/system/k12_diffvars.h}}). The global variable g_use_dersyntax defined in {{srcbranch|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.

[[File:Use_der_syntax.png|200px|thumb|center|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.

[[File:Derinfo.png|200px|thumb|center|The DerInfo window]]

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

<source lang=py>
x = M.x
y = M.y
print "Der(x,y) = ",float(der(x,y))
</source>

where <tt>M</tt> is a simulation Python object (see {{srcbranch|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.

-- [[User:Jpye|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.

--[[User:Ksenija|Ksenija]] 00:21, 11 February 2012 (EST)

[[Category:ASCEND Contributors]]
[[Category:GSOC2012]]

File:Solow plot.png

2013-04-23T14:14:25Z

Ksenija:

File:Solow.png

2013-04-23T14:12:17Z

Ksenija: