bugGNU Octave - Bugs: bug #59063, Event Handling different between...

 
 

bug #59063: Event Handling different between Ode-Solvers (especially ode15s)

Submitter:  None
Submitted:  Fri 04 Sep 2020 11:13:21 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
Status:  Fixed Assigned to:  None
Originator Name:  Michael Erhard Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * dev
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Fri 13 Nov 2020 08:07:52 PM UTC, comment #14: 

I slightly re-grouped the conditions and pushed your patch here:
https://hg.savannah.gnu.org/hgweb/octave/rev/b1b09b88585a

Marking as ready for test.

Markus Mützel <mmuetzel>
Group administrator
Fri 13 Nov 2020 05:57:21 PM UTC, comment #13: 

Attached the tested patch.

(file #50271)

Hg200 <hg200>
Thu 12 Nov 2020 10:33:50 PM UTC, comment #12: 

I think the longer change is fine as well.  You should probably add a comment ahead of the code that explains the code a little bit because this is quite complicated now.

Rik <rik5>
Group administrator
Thu 12 Nov 2020 06:48:11 PM UTC, comment #11: 

Thank you for the feedback. Here are some thoughts on this:

I have tested a lot and can confirm that the mentioned code section causes the problem. The suggested fix in comment #10 solves the problem in principle. But I wonder to what extent "corner cases" have to be considered. I have tested the following cases with Matlab R2018b. "val" refers to the "value" - variable returned by the event callback function:

Case 1).
"direction = 0" and "val = [0, 0, 0, ...]"
--> Matlab continues with the integration (does not stop). My interpretation is that "val" is not monotic (neither increasing nor decreasing)

Case 2).
"direction = 0" and "val = [-2, -1, 0, ...]"
--> Matlab stops at (and includes) the value of i = 3. "val" is "increasing". A similar case as the function would "decrease".

Case 3).
"direction = +1" (locates only zeros where the event function is increasing)
and "val = [1, 1, 0, 0, 0, 1, 1, 1, 1, ...]".
--> R2018b outputs the values for the indices i = 1..5. The value i=6 (equal to 1) is rejected.

I do not know much about ODE, but I can imagine that somebody wants to control the integration e.g. by a state machine (e.g. if clause), which returns boolean values {0, 1} only. In this case it would matter if we stop the integration in case 3) or not. The fix in comment #10 does not do this.

The following patch also catches case #3, but does not discard the value of i=6.


if ((val(i) >= 0 && oldval(i) < 0 && dir(i) != -1) // increasing
   || (val(i) <= 0 && oldval(i) > 0 && dir(i) != 1) // decreasing
   || (val(i) > 0 && oldval(i) <= 0 && dir(i) != -1) // increasing
   || (val(i) < 0 && oldval(i) >= 0 && dir(i) != 1)) // decreasing
 {
    index.resize (index.numel () + 1);
    index (index.numel () - 1) = i;
 }


Please let me know if this too hypothetically.

Hg200 <hg200>
Wed 11 Nov 2020 11:30:13 PM UTC, comment #10: 

Yes, it seems like this simple change might be enough:


-            if ((val(i) > 0 && oldval(i) < 0 && dir(i) != -1) // increasing
-                || (val(i) < 0 && oldval(i) > 0 && dir(i) != 1)) // decreasing
+            if ((val(i) >= 0 && oldval(i) < 0 && dir(i) != -1) // increasing
+                || (val(i) <= 0 && oldval(i) > 0 && dir(i) != 1)) // decreasing


Attached as a diff.

(file #50254)

Rik <rik5>
Group administrator
Wed 11 Nov 2020 09:00:20 PM UTC, comment #9: 

I have taken a closer look at the time based event problem. A minimal test code showing the problem is attached. In the attached example the integration should stop after 1 second, but it does not. Debugging in _ode15_.c shows that during integration "IDA::event" does not return a "true" status when the event value in the callback function becomes zero (_ode15_.c around 653).

Further investigation of IDA::event shows that the "val" - values from the callback function are available and that they are complete. The problem is rather in line 789 in IDA::event:


if ((val(i) > 0 && oldval(i) < 0 && dir(i) != -1) // increasing
     || (val(i) < 0 && oldval(i) > 0 && dir(i) != 1)) // decreasing


If the event hits exactly the output timestamp, "val" or "oldval" becomes exactly zero. But then this code falls through, because a crossing through "zero" is not recognized. An equal check "=" check is missing for this. I made a "quick" build check and this seems to solve the problem. If i find some time i will prepare a patch.

(file #50253)

Hg200 <hg200>
Sat 07 Nov 2020 09:07:55 PM UTC, comment #8: 

Sorry. I wasn't aware this bug was about two different issues.
Re-opening since I can't tell if both are fixed now.

Markus Mützel <mmuetzel>
Group administrator
Sat 07 Nov 2020 03:49:06 PM UTC, comment #7: 

Can this bug be closed yet?  It seems to me that there were two issues.  The first was 1-based indexing which is now resolved.  The second was that time-based events aren't captured if the event time and the output step size are exactly the same.

Rik <rik5>
Group administrator
Fri 06 Nov 2020 07:06:37 PM UTC, comment #6: 

I changed the respective BISTs and pushed your patch here:
https://hg.savannah.gnu.org/hgweb/octave/rev/1aa69571a313

Closing as fixed.

Markus Mützel <mmuetzel>
Group administrator
Tue 08 Sep 2020 08:17:08 AM UTC, comment #5: 

Thank you for the information and the submitted patch,
but I'm sorry, I would call me an "end user" only using the windows installer, infrequently looking in the maintainers/news blog, hopefully waiting for new releases with new features. Using Octave in a Matlab Simscape like manner is greatful (the stiff solvers ode12s/ode15i with event handling introduced in v4.4.1 brought the functionality to a new level), although there are some issues related to solver setup between stiff and nonstiff problems.

With the knowledge of these issues, making a probable working code is feasible for me. I'll waiting for a new release.

Anonymous
Tue 08 Sep 2020 03:09:13 AM UTC, comment #4: 

I won't be able to make much progress as I don't have access to a new enough version of SUNDIALS to work with Octave.  I took a stab at a patch for the zero-based indexing of the "ie" variable which is attached.  Perhaps you can apply it, compile, and see if it fixes the issue for you.

(file #49763)

Rik <rik5>
Group administrator
Tue 08 Sep 2020 01:17:11 AM UTC, comment #3: 

A quick answer to your last question, ode23s will be a part of version 6.1 and it uses the same event handler function as ode23/ode45 so it should work properly.

Rik <rik5>
Group administrator
Mon 07 Sep 2020 08:59:04 AM UTC, comment #2: 

Thanks for your quick response,
please find below a short (maybe physical meaningless) example of an oscillator with 2 events.
The first event is a state event (oscillator hits the ground-rebounding) and the second event is a time event (every second change in velocity).

The state event is no problem for all solvers, the only "bad" thing  is the change of index count (which is not documentated).
The time event isn't captured by ode15s, if the output step size is exactly an integer of the event time (here 1 s and 1.0e-2 s), here only ode23/ode45 works. Only if I change from 1.0e-2 s to 1.0001e-2 s ode15s gets the time events).

I think, that I read something about availability of ode23s in version 6.1. Is there an equal behaviour present?

Please find below a file with the requested example.

(file #49752)

Anonymous
Fri 04 Sep 2020 05:17:32 PM UTC, comment #1: 

Very likely true.  I worked on ode23/45 to have EventFcn work similarly to Matlab, but ode15i/15s have a different architecture and were written by a different author.

Do you have a simple test case for items #1 and #2?  That makes debugging easier and also allows verification once a fix is proposed.


Rik <rik5>
Group administrator
Fri 04 Sep 2020 11:13:21 AM UTC, original submission:  

In cases using an event function, ode23/ode45 gives identical behaviour, whereas ode15s (ode15i not yet tested) behaves different.

1st)
Indices of triggered events in EventFcn are counted from 1,2,3... in ode23/ode45, whereas ode15s count from 0,1,..., if multiple constraints are used.

2nd)
ode15s is not capable of detecting discrete time events, when the constraint time is a multiple of the output step size. ode23/ode45 doesn't have any problems with this

An example can be the bouncing ball problem by changing the event constraint from state to a discrete time before ball hits the ground or by adding an additional constraint for time.

Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #50271:  59063.event_fix.cset added by hg200 (1KiB - application/octet-stream)
file #50254:  59063.diff added by rik5 (777B - text/x-patch)
file #50253:  odemwe.m added by hg200 (1KiB - text/x-objcsrc)
file #49763:  59063.ie_fix.cset added by rik5 (804B - application/octet-stream)
file #49752:  ODEwithEvents.m added by None (3KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by hg200 (Updated the item)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 17 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2020-11-17 rik5 StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2020-11-13 mmuetzel StatusConfirmed Ready For Test
    2020-11-13 hg200 Attached File- Added 59063.event_fix.cset, #50271
    2020-11-12 rik5 StatusNone Confirmed
    2020-11-11 rik5 Attached File- Added 59063.diff, #50254
    2020-11-11 hg200 Attached File- Added odemwe.m, #50253
    2020-11-07 mmuetzel StatusFixed None
        Open/ClosedClosed Open
    2020-11-06 mmuetzel StatusNeed Info Fixed
        Open/ClosedOpen Closed
    2020-09-08 rik5 Attached File- Added 59063.ie_fix.cset, #49763
    2020-09-07 None Attached File- Added ODEwithEvents.m, #49752
    2020-09-04 rik5 Item GroupIncorrect Result Matlab Compatibility
        StatusNone Need Info
        Release5.2.0 dev
        Operating SystemMicrosoft Windows Any

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code