Skip to content

Commit

Permalink
[NB] Collect composite time events (#9234)
Browse files Browse the repository at this point in the history
Events of the form
  sample(t0, dt) and (f(x) > 0)
do not need to introduce zero crossings and can be treated
as time events.
  • Loading branch information
phannebohm committed Aug 1, 2022
1 parent c855dd5 commit d78c600
Show file tree
Hide file tree
Showing 4 changed files with 220 additions and 63 deletions.
111 changes: 109 additions & 2 deletions OMCompiler/Compiler/NBackEnd/Modules/2_Pre/NBEvents.mo
Expand Up @@ -399,7 +399,7 @@ public
end createSample;

function createSampleTraverse
"used only for StateEvent traversel to encapsulate sample operators"
"used only for StateEvent traversal to encapsulate sample operators"
input output Expression exp "has to be LBINARY() with comparing operator or a sample CALL()";
input output Bucket bucket "bucket containing the events";
algorithm
Expand All @@ -414,6 +414,111 @@ public
end match;
end createSampleTraverse;

function createComposite
"Find special events of the form: sample(t0, dt) and (f(x) > 0)
These events can only occur at the sample times. At that time the additional condition
is checked only once, no state event necessary!
NOTE: This does not work for SIMPLE_TIME, e.g. (time > 0.2) and (f(x) > 0)"
input output Expression condition;
input output Bucket bucket;
output Boolean failed = false "returns true if composite event list could not be created";
protected
Pointer<Variable> aux_var;
ComponentRef aux_cref;
algorithm
(condition, bucket, failed) := match condition
local
Expression exp, exp2;
Call call;

// base case: sample is the left operand to AND
case Expression.LBINARY(exp1 = exp as Expression.CALL(call = call), operator = Operator.OPERATOR(op = NFOperator.Op.AND))
guard BackendUtil.isOnlyTimeDependent(exp)
algorithm
(call, exp2, bucket, failed) := checkDirectComposite(call, condition.exp2, bucket);
if not failed then
exp.call := call;
condition.exp1 := exp;
if not referenceEq(exp2, condition.exp2) then
condition.exp2 := exp2;
end if;
end if;
then (condition, bucket, failed);

// base case: sample is the right operand to AND
case Expression.LBINARY(exp2 = exp as Expression.CALL(call = call), operator = Operator.OPERATOR(op = NFOperator.Op.AND))
guard BackendUtil.isOnlyTimeDependent(exp)
algorithm
(call, exp2, bucket, failed) := checkDirectComposite(call, condition.exp1, bucket);
if not failed then
exp.call := call;
condition.exp2 := exp;
if not referenceEq(exp2, condition.exp1) then
condition.exp1 := exp2;
end if;
end if;
then (condition, bucket, failed);

// recursion: sample might be nested (all parent operators have to be AND)
// e.g. (sample(t0, dt) and f1(x)) and f2(x)
case Expression.LBINARY(operator = Operator.OPERATOR(op = NFOperator.Op.AND))
algorithm
(exp, bucket, failed) := createComposite(condition.exp1, bucket);
if not failed then
condition.exp1 := exp;
(exp, bucket, failed) := createComposite(condition.exp2, bucket);
if not failed then
// TODO what if there is more than one sample()?
condition.exp2 := exp;
end if;
failed := false; // we know we have a composite time event in the first half
else
(exp, bucket, failed) := createComposite(condition.exp2, bucket);
if not failed then
condition.exp2 := exp;
end if;
end if;
then (condition, bucket, failed);

else (condition, bucket, true);
end match;

if not failed then
if TimeEventTree.hasKey(bucket.timeEventTree, condition) then
// time event already exists, just get the identifier
aux_var := TimeEventTree.get(bucket.timeEventTree, condition);
aux_cref := BVariable.getVarName(aux_var);
else
// make a new auxiliary variable representing the state
(aux_var, aux_cref) := BVariable.makeEventVar(NBVariable.TIME_EVENT_STR, bucket.auxiliaryTimeEventIndex);
bucket.auxiliaryTimeEventIndex := bucket.auxiliaryTimeEventIndex + 1;
// add the new event to the tree
bucket.timeEventTree := TimeEventTree.add(bucket.timeEventTree, condition, aux_var);
// also return the expression which replaces the zero crossing
end if;
condition := Expression.fromCref(aux_cref);
end if;
end createComposite;

function checkDirectComposite
"Checks if call is a sample call and if it is creates the appropriate events.
Also checks the rest exp for composite events, not sure if this is necessary."
input output Call call "sample call";
input output Expression exp;
input output Bucket bucket;
output Boolean failed;
protected
Boolean failed2;
algorithm
(call, bucket, failed) := createSample(call, bucket);
if not failed then
(exp, bucket, failed2) := createComposite(exp, bucket);
if not failed2 then
// TODO what if there is more than one sample()? Can we simplify this?
end if;
end if;
end checkDirectComposite;

function getIndex
input TimeEvent timeEvent;
output Integer index;
Expand Down Expand Up @@ -763,9 +868,11 @@ protected
protected
Boolean failed = true;
algorithm
// try to create time event
// try to create time event or composite time event
if BackendUtil.isOnlyTimeDependent(condition) then
(condition, bucket, failed) := TimeEvent.create(condition, bucket);
else
(condition, bucket, failed) := TimeEvent.createComposite(condition, bucket);
end if;

// if it failed create state event
Expand Down
@@ -1,6 +1,7 @@
TEST = ../../../../rtest -v

TESTFILES = \
compositeEvent.mos \
eventSystem.mos \
hybridSys.mos \

Expand Down
@@ -0,0 +1,49 @@
// name: compositeEvent
// keywords: NewBackend
// status: correct

loadString("
model CompositeEvent
Real x(start=0, fixed=true);
Integer n(start=0, fixed=true);
equation
der(x) = 1;
when x < 1 and (x > 0.5 and sample(0.125,0.25)) then
n = pre(n) + 1;
end when;
end CompositeEvent;
"); getErrorString();
simulate(CompositeEvent, simflags="-lv LOG_EVENTS"); getErrorString();
val(n,0.125);
val(n,0.375);
val(n,0.625);
val(n,0.875);

// Result:
// true
// ""
// true
// ""
// record SimulationResult
// resultFile = "/mnt/home/ph/ws/models/compositeEvent/CompositeEvent_res.mat",
// simulationOptions = "startTime = 0.0, stopTime = 1.0, numberOfIntervals = 500, tolerance = 1e-06, method = 'dassl', fileNamePrefix = 'CompositeEvent', options = '', outputFormat = 'mat', variableFilter = '.*', cflags = '', simflags = ''",
// messages = "LOG_EVENTS | info | status of relations at time=0
// LOG_EVENTS | info | status of zero crossings at time=0
// LOG_SUCCESS | info | The initialization finished successfully without homotopy method.
// LOG_EVENTS | info | time event at time=0.125
// | | | | | [1] sample(0.125, 0.25)
// LOG_EVENTS | info | time event at time=0.375
// | | | | | [1] sample(0.125, 0.25)
// LOG_EVENTS | info | time event at time=0.625
// | | | | | [1] sample(0.125, 0.25)
// LOG_EVENTS | info | time event at time=0.875
// | | | | | [1] sample(0.125, 0.25)
// LOG_SUCCESS | info | The simulation finished successfully.
// "
// end SimulationResult;
// ""
// 0.0
// 0.0
// 1.0
// 2.0
// endResult
122 changes: 61 additions & 61 deletions testsuite/simulation/modelica/NBackend/event_handling/hybridSys.mo
@@ -1,71 +1,71 @@
function func1
input Real x;
output Real func1_out;
input Real x;
output Real func1_out;
algorithm
func1_out := x;
func1_out := x;
end func1;

function func2
input Real x1 ;
input Real x2 ;
output Real func2_out;
input Real x1;
input Real x2;
output Real func2_out;
algorithm
func2_out := x1 + x2;
func2_out := x1 + x2;
end func2;

model hybridSys
parameter Integer Niter=4;
//Variablesofthediscreteeventmodel
Boolean phase_Start(start=true);
Boolean phase_Loop1(start=false);
Boolean phase_Loop2(start=false);
Boolean phase_Loop3(start=false);
Boolean phase_End(start=false);
Real x_Start(start=0);
Real x_Loop1(start=0);
Real x_Loop2(start=0);
Real x_Loop3(start=0);
Real x_End(start=0);
Boolean startCondition(start=false);
Boolean loopCondition1(start=false);
Boolean loopCondition2(start=false);
Boolean loopCondition3(start=false);
Boolean endCondition(start=false);
//Variablesofthecontinuous-timemodel
Real x1(start=10),x2;
parameter Integer Niter=4;
// Variables of the discrete event model
Boolean phase_Start(start=true);
Boolean phase_Loop1(start=false);
Boolean phase_Loop2(start=false);
Boolean phase_Loop3(start=false);
Boolean phase_End(start=false);
Real x_Start(start=0);
Real x_Loop1(start=0);
Real x_Loop2(start=0);
Real x_Loop3(start=0);
Real x_End(start=0);
Boolean startCondition(start=false);
Boolean loopCondition1(start=false);
Boolean loopCondition2(start=false);
Boolean loopCondition3(start=false);
Boolean endCondition(start=false);
// Variables of the continuous-time model
Real x1(start=10),x2;
equation
//---------------------
//Continuous-timemodel
der(x1)=-func1(x1);
//Nodiscrete-to-continuousinteraction
//x2=func1(x1);
//Discrete-to-continuousinteraction
x2=func2(x1,x_End);
//--------------------
//Discrete-eventmodel
startCondition=time>1;
loopCondition1=pre(x_Loop1)<Niter+1;
loopCondition2=pre(x_Loop2)<Niter+1;
loopCondition3=pre(x_Loop3)<Niter;
endCondition= not loopCondition3;
phase_Start=pre(phase_Start) and not startCondition;
phase_Loop1=pre(phase_Start) and startCondition or pre(phase_Loop3) and loopCondition3 or pre(phase_Loop1) and not loopCondition1;
phase_Loop2=pre(phase_Loop1) and loopCondition1 or pre(phase_Loop2) and not loopCondition2;
phase_Loop3=pre(phase_Loop2) and loopCondition2 or pre(phase_Loop3) and not (loopCondition3 or endCondition);
phase_End=pre(phase_Loop3) and endCondition;
when phase_Start then
x_Start=pre(x_Start)+1;
end when ;
when phase_Loop1 then
x_Loop1=pre(x_Loop1)+1;
end when ;
when phase_Loop2 then
x_Loop2=pre(x_Loop2)+1;
end when ;
when phase_Loop3 then
x_Loop3=pre(x_Loop3)+1;
end when ;
when phase_End then
x_End=pre(x_End)+1;
end when ;
end hybridSys;
//---------------------
// Continuous-time model
der(x1) = -func1(x1);
// No discrete-to-continuous interaction
//x2 = func1(x1);
// Discrete-to-continuous interaction
x2 = func2(x1,x_End);
//--------------------
// Discrete-event model
startCondition = time>1;
loopCondition1 = pre(x_Loop1)<Niter+1;
loopCondition2 = pre(x_Loop2)<Niter+1;
loopCondition3 = pre(x_Loop3)<Niter;
endCondition = not loopCondition3;
phase_Start = pre(phase_Start) and not startCondition;
phase_Loop1 = pre(phase_Start) and startCondition or pre(phase_Loop3) and loopCondition3 or pre(phase_Loop1) and not loopCondition1;
phase_Loop2 = pre(phase_Loop1) and loopCondition1 or pre(phase_Loop2) and not loopCondition2;
phase_Loop3 = pre(phase_Loop2) and loopCondition2 or pre(phase_Loop3) and not (loopCondition3 or endCondition);
phase_End = pre(phase_Loop3) and endCondition;
when phase_Start then
x_Start = pre(x_Start)+1;
end when;
when phase_Loop1 then
x_Loop1 = pre(x_Loop1)+1;
end when;
when phase_Loop2 then
x_Loop2 = pre(x_Loop2)+1;
end when;
when phase_Loop3 then
x_Loop3 = pre(x_Loop3)+1;
end when;
when phase_End then
x_End = pre(x_End)+1;
end when;
end hybridSys;

0 comments on commit d78c600

Please sign in to comment.