Online Cycle Detection for Models with Mode-Dependent Input and Output Dependencies
OOnline Cycle Detection for Models withMode-Dependent Input and Output Dependencies
Heejong Park a , Arvind Easwaran a , Etienne Borde b a Nanyang Technological University, Singapore b LTCI, T´el´ecom Paris, Institut polytechnique de Paris, France
Abstract
In the fields of co-simulation and component-based modelling, designersimport models as building blocks to create a composite model that providesmore complex functionalities. Modelling tools perform instantaneous cycledetection (ICD) on the composite models having feedback loops to rejectthe models if the loops are mathematically unsound and to improve simula-tion performance. In this case, the analysis relies heavily on the availabilityof dependency information from the imported models. However, the cycledetection problem becomes harder when the model’s input to output de-pendencies are mode-dependent, i.e. changes for certain events generatedinternally or externally as inputs. The number of possible modes createdby composing such models increases significantly and unknown factors suchas environmental inputs make the offline (statical) ICD a difficult task. Inthis paper, an online ICD method is introduced to address this issue for themodels used in cyber-physical systems. The method utilises an oracle as acentral source of information that can answer whether the individual mod-els can make mode transition without creating instantaneous cycles. Theoracle utilises three types of data-structures created offline that are adap-tively chosen during online (runtime) depending on the frequency as wellas the number of models that make mode transitions. During the analysis,the models used online are stalled from running, resulting in the discrepancywith the physical system. The objective is to detect an absence of the instan-taneous cycle while minimising the stall time of the model simulation thatis induced from the analysis. The benchmark results show that our methodis an adequate alternative to the offline analysis methods and significantlyreduces the analysis time.
Keywords:
Instantaneous cycle, modelling, cyber-physical system,
Preprint submitted to Journal of Systems Architecture February 4, 2021 a r X i v : . [ c s . S E ] F e b imulation, causality loop
1. Introduction
Today, there exist a plethora of modelling tools for capturing the be-haviour of physical systems. Many of these tools provide a set of well-established libraries that can be reused for constructing more complex sys-tems. Reusability is a significant part of many engineering tasks where de-signers can rely on the correctness of individual components, which havebeen used and tested extensively by others. A library of well-defined compo-nents, however, does not always guarantee the correctness of a model madeof such components. For example, one of the useful techniques in modellingis a feedback loop where output events of a system are routed back as itsinputs and create cyclic dependencies. Nevertheless, creating feedback loopsin an undisciplined manner can result in instantaneous cycles . It refers tothe situation when the cause of input and output events for a model are in-terdependent with each other that can make composite models (models thatare made of other smaller models) mathematically unsound. This leads toa simulation that diverges from the operations of the physical system andincreases the simulation time of the composite model [1].The classical approach for checking the absence of instantaneous cyclesin a composite model is to reject any feedback loops that are not brokenby a unit delay. A delay is typically introduced in a model via componentssuch as an integrator and unit delay components whose output events do notdepend on their present or future input events at any time instance. Thus, amodelling tool performs instantaneous cycle detection (ICD) during designtime (offline) based on the model composition and statically rejects mod-els with feedback loops whose computational dependency cannot be solvedmathematically. Yet, in the environment where models are imported as third-party components [2], detecting cycles in a composite model can be difficult.This is mainly because the library vendors hide the internal implementationof the models due to Intellectual Property (IP) issues, prohibiting accurateanalysis of the design. Another approach consists in rejecting any modelcompositions that structurally appear as feedback loops on a top-level orautomatically inserting a delay component in every loop. Yet, such methodscan potentially reject models that are acyclic when the models are flattenedto the lowest level of the hierarchy [3].2CD becomes more difficult when the input to output dependencies ofmodels dynamically change based on their internal mode (like in hybrid au-tomata). There is a body of work [4, 5, 1] which employ causal dependency information of individual sub-model’s input and output ports to check in-stantaneous cycles in a composite model. Most often this information isprovided by the library vendors or the tools that generate models so that thetool that performs the assembly of components perform ICD. For a compositemodel consisting of mode-dependent sub-models, the analysis should verifythat all reachable modes exclude instantaneous cycles. However, staticallyfiltering out such erroneous modes in a composite model through a reacha-bility analysis is typically an undecidable problem [5] for white-box modelsand even infeasible for component-based models. These components are typ-ically grey-box models whose mathematical behaviours are hidden and onlypartial information, such as input to output dependency, is exposed to theirexternal environment.In this paper, we instead tackle the mode-dependent ICD problem viaan online method that provides an efficient way to reduce the stall time ofa component-based model simulation due to the analysis. The stall time isthe duration that the models are blocked from execution due to housekeep-ing works such as ICD. We target areas of cyber-physical systems includingbut not limited to: model driven control systems where the models continu-ously synchronise with the physical system such that the results of the modelsimulation influence operations of the physical system or vice-versa. There-fore, a large stall time causes an unwanted deviation of models from theirrepresented physical counterparts. To our knowledge, there is no literatureto this date that tries to address the overhead due to the online
ICD for mode-dependent input to output dependencies in grey-box component-based models. In particular, our method employs an oracle , which can pro-vide ‘yes’ or ‘no’ answers to individual models requesting a mode changeauthorisation. Therefore, our approach does not need to verify all reachablemodes combinations from the model composition but detects an instanta-neous cycle online whenever the models are about to change their modes.The key problem in our approach is to reduce the stall time required for theoracle to decide if the requested mode changes are acceptable therefore doesnot create discrepancies between the models and the physical system. Thedetailed contributions of this paper are:1. An online ICD method for Cyber-Physical Systems (CPS) consisting3 ecurecommunicationchannel Mode change requestModel 1 Model 2ApprovalOracle
Sync
Cyber Space PhysicalSpace
Figure 1: Oracle-based instantaneous cycle detection of component-based models whose input and output dependencies aremode-dependent.2. An adaptive method to reduce the stall time (time consumed by ICDalgorithm) of models that synchronise with physical systems in a real-time manner and thus can be used in various CPS applications.3. A set of benchmark results and a case study that shows the applicabilityof our method in the fields of industrial manufacturing systems.The rest of this paper is organised as follows: Section 2 defines the prob-lem that we address in this paper. Section 3 enumerates the related researchon ICD followed by the introduction of the workpiece sorting system as ourmotivating example in Section 4. Preliminary background is presented inSection 5. The methodology of our online ICD is introduced in Section 6.Sections 7 and 8 present a set of experimental results that evaluate the per-formance of the analysis in various settings. A set of use-case scenarios ofour technique is presented in Section 9 through existing real-world examples.Finally, conclusion and future work are given in Section 10.
2. Problem Definition
A graphical overview of the proposed oracle-based online ICD is shownin Fig. 1. We consider models as a set of components with their input andoutput ports interconnected with each other to exchange data during thesimulation. Furthermore, the models simulated in the cyber space synchro-nise operations with their physical counterparts in the physical space. In the4yber space, the arrows indicate a data dependency from output to input between the models 1 and 2 or vice-versa within the models. For example, togenerate output data from the output port b from the model 1, it requires aninput value from port a , which is generated from the output port e from themodel 2. When there is no arrow incident to an output port, for example e , itcan generate data without the need for the input port values to be resolved.The composition of these models shown in the figure is acyclic because thereare no paths along the input and output dependencies that create a cycle,i.e. creating an instantaneous cycle. In addition, we also consider the modelsthat can change the input to output dependencies along with their operatingmodes. For example, creating a dependency between the ports d and e upona mode change in model 2. In this case, an instantaneous cycle is created forthe path e → a → b → d → e . Our objective is to quickly detect the cre-ation of such a cycle during runtime due to mode changes in the models. Toachieve this, we assume the models make mode change requests (MCR) to thecentral oracle, which monitors the creation of the instantaneous cycles beforethe mode change can happen. Each MCR from a model only contains localdependency information and the oracle ensures this IP-sensitive informationis not shared among the models. The oracle checks whether the mode changeis valid and responds with an approval message back to the models within afinite time bound. We call this bound stall time ( t stall ) because the modelsmust stall and not proceed to the next simulation step until they receive aresponse from the oracle for their MCR requests. Stalling models do notrequire additional features in these models since their execution is controlledby an external entity, e.g. main loop, that coordinates progression of themodel’s time. In this work, we consider grey-box, component-based modelsthat expose their current mode’s input-to-output dependency information tothe environment. Here we define the problem as follows. Problem 1.
Develop an online ICD technique that minimises the stall timefor grey-box, component-based models whose input and output dependenciescan change during runtime depending on their internal modes.
When the oracle detects an instantaneous cycle from the MCR requests,it rejects the request from the model that creates the loop. Such rejectionshould not result in a fault in an on-going cyber-physical operation. There-fore, the corresponding model makes a transition to a safe mode , which isspecified by a designer beforehand and guaranteed not to create an instan-taneous cycle. These safe modes would require additional logic in models;5owever, we think this is a little addition to the models that support multi-mode features which we target in this paper. We also assume the modelsexpose a limited set of safe modes and their combinations with other mod-els are checked statically offline. The size of such combination is typicallymuch smaller than the entire combination of all possible modes. Models thatmade a transition to a safe mode can return to a normal operating mode inthe same way they make an MCR to the oracle. In this way, the analysiscan be performed independently by a trusted entity (oracle) and models areonly required to provide changes in the local input-to-output dependencyinformation for each mode change.
3. Related Work
The problem of detecting cyclic dependencies is a subset of causality anal-ysis problem, which can be found in formal modelling literature. Authorsin [5] introduced dependency algebra to formulate the causality problem fora network of actors whose input and output dependencies are fixed. Depen-dency algebra is implemented in Ptolemy II [6] to check the instantaneouscycle via solving the algebraic equation, which is derived from the composi-tion of actors. Their method checks if every simple cycle in the network ofactors is broken by the delay component (e.g. the integrator block). Findingall simple cycles, however, requires traversals of the whole actor network asmany times as there are cycles in the network. Ptolemy II also provides an op-tion to support the mode-dependent online ICD where only input-to-outputdependencies in active modes are considered for the analysis. However, thisoption still employs the method introduced in [5] and thus has the sameruntime complexity as solving algebraic equations for the actors with fixedmodes. On the other hand, we are introducing an adaptive technique to ef-ficiently detect the mode-dependent instantaneous cycles that is faster thanfinding simple cycles on the whole network.Authors in [7] introduced a causal analysis method for concurrent hybridautomata. Similar to our approach, their approach tries to detect causalitycycles online. They introduced a notion of compatibility where two automataare compatible if no output variables are shared in any mode combinations.Nevertheless, checking such compatibility is an expensive task when the num-ber of all possible mode combinations grows exponentially. Our method, onthe other hand, does not require checking of all possible mode combinationsamongst concurrently running automata, therefore more amenable in the6 j e c t o r E j e c t o r E j e c t o r Bin Bin BinSensor Sensor Sensor EntryExit C CC C
Figure 2: A workpiece sorting system online setting. Work on the structural analysis of multi-mode differentialalgebraic equation (DAE) systems is presented in [8]. The authors showedan example of modelling a clutch in a car between engaged and releasedmodes. The approach is based on the analysis on the equation (program)level which is hard to directly apply to our case that utilises input and outputdependencies on the component level.The problem of cycle detection especially in a dynamically changing graphdataset has been applied in a variety of applications. In [9], authors intro-duced a technique to find cycles in a large-scale graph that satisfy both lengthand some attribute constraints. The technique is deployed at Alibaba in ane-commerce system to monitor fraudulent activities upon cycle detection.Online cycle detection technique is also used in the pointer analysis in aprogram [10] as well as in a distributed deadlock detection algorithm [11].In this paper, on the other hand, we introduce the dynamic cycle detectionproblem in the domain of component-based modelling and present an efficienttechnique by classifying the types of mode changes within the models.We foresee our approach can be applied in a variety of applications suchas checking dependencies in a task mapping scenarios [12], dependency reduc-tion algorithms [13] and fault detection algorithms that have task executiondependencies as constraints [14].
4. Motivating Example
Consider an example of an industrial automation system shown in Fig. 2.It consists of a set of ejectors that extend to place workpieces into one ofthe bins in front of them. The linear conveyor belt can transfer workpiecesto the next workstation via an exit point if they are yet to be consideredas final products. Each machine in the system is controlled by its controller7 ushPullTokenOutEjEndEjStartTokenIn WpInWPPos
ControllerModel EjectorModel
Figure 3: An ejector system model consisting of a controller and a plant that runs in thecyber space (annotated with ‘C’) which makes them a closed-loop system. In addition, weassume the case where both the controller and the machine (plant) modelsrun concurrently in real-time in the cyber-space. Without further goinginto implementation details of each model, Fig. 3 illustrates the top-levelcomponents of the controller and the ejector, which are interconnected viainput (white rectangle) and output (black rectangle) ports. The semantics ofcommunication between these two models are depending on their model ofcomputation (MoC). In this paper, we assume these models are based on thesynchronous reactive (SR) MoC [15] with an extension of continuous time,where discrete events between these models are instantaneously propagatedvia input and output ports within the same simulation time step called a tick .Fig. 4 shows the implementation of the controller and plant of the workpiecesorting system in Fig. 3 in an SR language called Esterel [16].The controller model has two outputs
Push and
Pull that push and pullthe ejector to move the workpiece into the bin. Feedback signals
EjStart and
EjEnd provided by the ejector model indicate if the ejector has beencompletely retracted or extended, respectively. The input
WpPos is a con-tinuous variable that indicates the locations of the workpieces derived fromthe conveyor belt model (not shown in the figure).
WpIn is a sensory in-put to the controller that triggers it to issue a push signal to the ejector.
TokenIn and
TokenOut are connected between adjacent controllers to im-plement a ring-token, which is used to evenly distribute workpieces into thebin. In Fig. 4, the controller (lines 1-14) and ejector (lines 14-32) models areexecuted concurrently with each other indicated by the synchronous paral-lel operator || (line 14). In this program, we implement modes using theconditional present statement where each branch indicates a single mode.8 [ % Controller Model loop present TokenIn then present notFull then % Mode -1 ( Push ) await [ EjStart and WpIn ]; abort sustain Push when EjEnd ; await EjEnd ; abort sustain Pull when EjStart ; end present ; emit TokenOut ; end present ; pause end loop ] || [ % Plant Model var pos : float in loop present notFull then % Mode -1 present case Push do pos := integrate ( pos , 1) case Pull do pos := integrate ( pos , 0) end present ; % 1. Check the location of WP ( WpPos ) % 2. Emit EjEnd , EjStart based on pos else % Mode -2 present Push then emit EjStart % Cycle end present end present ; pause end loop end var ] Figure 4: Esterel implementation of the controller and plant models for the workpiecesorting system.
These models have several distinct modes to adapt to different operatingscenarios. The cyclic dependency problem occurs due to the feedback loopssuch as
EjStart , EjEnd and
WpIn as shown in this example. In many cases,it is not clear if these loops are instantaneous cycles in this top-level com-position. Indeed, it depends on the internal implementation of each model.Fig. 3 shows an example of input to output dependencies for these modelswhere e i denotes a dependency in a mode i . When the controller is in mode9, the output Push depends on both
EjStart and
WpIn because the con-troller can only trigger the operation when the ejector is fully retracted and the workpiece is detected via a sensor. This dependency is shown in Fig. 4at line 5 using the Esterel statements await that blocks the controller untilboth
EjStart and
WpIn become ‘present’ before outputting the signal
Push at line 6. The ejector model does not require
Push to generate
EjStart sincethe dynamics of the eject operation breaks this dependency, i.e. an integratorblock in the model breaks the dependency. This is shown at lines 18-23 inFig. 4 (some parts are omitted for brevity) that is performed solely on thevalue of pos . Therefore the model is acyclic in this mode .We can introduce an alternative mode for these models when the bins arefull so that a workpiece cannot be pushed by the ejectors. In this scenario,when the controller tries to push the workpiece via Push , EjEnd never be-comes true since the bin is fully occupied. Instead,
EjStart becomes falseimmediately until the controller stops issuing
Push . In this case, we cre-ated an instantaneous cycle since
EjStart is now immediately dependent on
Push , i.e. the path where edges of mode 2, noted e , are created in Fig. 3. InFig. 4, this dependency is shown at lines 25-27 where EjStart is only emittedif the signal
Push is present. The output of the Esterel compiler (v5.92) uponcompilation of this code indicates there exist a cycle between these signalsand the compiler cannot generate statically scheduled code. If we closelyinspect the code in Fig 4, such cycle cannot be created if the guarded signal notFull stays present while the control-flow of the controller model is withinlines 5-6, i.e. the plant model does not enter the lines 25-27 that creates acycle. However, such case could not be guaranteed at compile time and theprogram is rejected for code generation. Therefore, static analysis can stillreject models that are not cyclic based on the active modes. Furthermore,it would be even harder to detect cycles created due to the composition ofmodels with multiple modes, for example when the controller and plant mod-els are separately developed in component-based modelling environment asshown in Fig. 3. This example motivates the need for online technique thatdynamically checks the absence of instantaneous cycles.The controllers pass a token in a circular fashion to select which bin tostore an incoming workpiece. When a bin is full, the corresponding controller Note that dependencies between
Pull and
EjEnd are also similar but omitted in thefigure for the sake of conciseness e between TokenIn and
TokenOut is created to immediately pass a token to the next controller. Therefore aninstantaneous cycle can also be created in this case such that all controllerspass the token to each other indefinitely. The frequency of the mode changesfor the ejector models can cause a lag in the model simulation time due toadditional analysis times and result in discrepancies between the models andits physical counterparts.Offline analysis of the instantaneous cycle becomes more difficult when thesystem is larger and the total number of modes among the models is signif-icant, especially when the models exhibit concurrent behaviours. Moreover,the analysis has to introduce many pessimistic assumptions when the internalimplementations of the model’s behaviour are hidden due to IP restrictions.One possible way to mitigate this problem is to perform the analysis duringruntime on every mode change and provide safe mode transitions, which areknown to be correct, upon a ICD. The safe mode transition is required inthe applications where the models are synchronised with the physical systemas shown in Fig. 1, to prevent unwanted events, such as physical damagesin the system, that cannot be undone. One major advantage of this ap-proach is that there is no need to filter any erroneous modes that requiremany assumptions which may not be true. Nevertheless, this runtime analy-sis introduces an unwanted stall time that blocks the whole composite modelfrom execution, thus creating a discrepancy such as a delay in the operationbetween the model and the physical system. In this work, we are introducinga method to reduce such stall time.
5. Background
A composite model C in our modelling approach is a tuple (cid:104) M, I, O, S, T (cid:105) where M is a set of component-based models, I is a set of input ports, O is a set of output ports and S : O × I is a set of interface signals thatdescribe dependencies from output to input ports. A signal is a status anda value pair ( s st , s v ) where s st ∈ { , , ⊥} , s v ∈ R . A sequence of ticks T = { ( n, r ) | n ∈ N , r ∈ R , n ≥ ∧ r ≥ } is a shared time among allmodels within the same composite model. We extend the logical time n inthe traditional SR MoC with a real number r to describe continuous timemodels such as Ordinary Differential Equations (ODEs).11he time of a composite model C progresses when all the signal statuses s st are resolved from unknowns ⊥ to either 0 or 1. In the presence of theinstantaneous cycle, not all signal statuses may be resolved that blocks themodel simulation from progressing to the next tick. This paper tackles thisproblem via an online method to prevent such a deadlock in composite modelsbeing executed with the SR MoC where the ICD is performed at every tickboundary. Given a directed acyclic graph (DAG) G = (cid:104) V, E (cid:105) where V is a set ofvertices and E ⊆ V × V , transitive closure (TC) of G denoted as R + isthe smallest transitive relation on V that includes E as a subset. ∀ ( u, v ) ∈ V × V, ( u, v ) ∈ R + iff there exist a path from u to v in E . Thereforetransitive closure of G gives a reachability relation between all vertices in V .One can build a TC between all input and output vertices in G via matrixmultiplication of an adjacency matrix A in O ( n ω ) time where ω = 2 .
38 [17].Dynamic transitive closure (DTC) maintains a data-structure that can beupdated upon insertion and addition of edges in the DAG. There are threemain operations for a DTC: • insert(x,y) – To add a transitive relation in a DTC between two verticesfrom x to y . • delete(x,y) – To remove a transitive relation in a DTC between twovertices from x to y . • query(x,y) – To check if a vertex y is reachable from x .Complexities of these operations vary depending on the implementation ofthe data-structure. Typically, improving performance of one operation de-grades the others and vice-versa [18, 19].Transitive reduction (TR) of G denoted as R − is a minimal set of E whosetransitive closure is identical to the transitive closure of G . In other words, itgives a graph that has the minimal number of edges with the same reachabil-ity relation as that of the original DAG. It is known that the time complexityof computing TR is in the same class as that of computing TC [20]. As itwill be shown in the later section, we employ both TC and TR for reducingthe stall time of the ICD. 12 nputProvider Process Outputs Designer
CompositeModel
Initialanalysis
Process Start
CGTRTCInitialiseOracleModels MCR InstantaneousCycleDetectionCycle? Yes
No(Example)
CPS ToolSafe ModeTransition
Inputs OutputConsumer
OfflineOnline
Figure 5: Workflow diagram for the oracle-based instantaneous cycle detection
6. The Oracle-based Online Instantaneous Cycle Detection
To enable the analysis of grey-box component-based models, we assumeinput-to-output dependencies are embedded in each model by the modellingtools such as in Functional Mock-up Unit [2]. In addition, we also assume themodels can provide the updated dependency information upon the changeof their modes. Since IP-related issue is one of the main reasons that themodel vendors might hide the implementation of their model’s behaviour,our method utilises a single trusted entity so-called an oracle. It thus limitsthe amount of information that must be (i) exchanged among model vendors,or (ii) provided to the models integrator, for the ICD.
The overall workflow of the ICD is shown in Fig. 5. The process is dividedinto two phases: offline and online, which are indicated by the blue and yellowboxes, respectively. At the beginning of the offline process, a composite modelcomprised of a set of interconnected sub-components via input and outputports is given as an input by a designer. This model is processed during theinitial analysis phase, which generates three types of data structures. The13
CR Y/NOICDOracle
Figure 6: An effect of the stall time for model simulation first two are a composite graph G c and its transitive closure indicated by G tc . The other is a set of graphs { G tr , ..., G ktr } for some k where each G itr isobtained from the transitive reduction on the vertex-partitioned sub-graphsof G tc . A brief discussion on how to partition G tc based on MCR types ispresented in Section 6.3.The outputs of the offline phase are given as inputs to the oracle at thebeginning of the online phase as indicated by the yellow box in Fig. 5. Whenmodels need to change their modes, the information on their input to outputdependency is transmitted to the oracle via an MCR message. The oraclechooses one of G c , G tc and G itr to perform the ICD and answers either ‘yes’or ‘no’ back to the models. The time between the MCR and the answergenerated by the oracle is called stall time , during which all the models areblocked and cannot advance to the next simulation step. This is illustratedin Fig. 6 where the MCR call is made to the oracle at the end of the currentsynchronous tick n and the subsequent execution is blocked for the durationof t stall due to online ICD (OICD). The construction of G c = (cid:104) V c , E c (cid:105) from the composite model C involvescreating vertices V c = I ∪ O and an edge set E c = S ∪ E L where E L is a setof input to output edges for the initial modes for all models m ∈ M and S is a set of interface signals as explained in Section 5.1. S is fixed throughoutthe life-time of C whereas E L changes along with the internal modes of themodels within C . The presence of cycles in this graph is considered as thepresence of instantaneous cycles in the composite model: Definition 1 (Instantaneous cycle) . A composite model C contains the in-stantaneous cycle if its composite graph G c contains any cycle. G c using depth-first search (DFS) whenever the structure of the graph changesdue to the mode changes in the models. However, as shown in [5], thecomplexity of DFS for the whole G c would be costly especially with the largenumbers of nodes and edges in the graph. Therefore we construct a transitiveclosure G tc of G c during the offline phase shown in Fig. 5 that can be usedfor efficiently querying the existence of instantaneous cycle in the compositemodel upon a mode change.Formally, G tc = (cid:104) V c , E tc (cid:105) is a transitive closure of G c . In this work, weadopt the technique introduced in [21] to compute G tc , which maintains theadjacency matrix explicitly. The adjacency matrix A xy of G tc contains thenumber of paths for all ( x, y ) ∈ E tc ; the existence of a path between twovertices in G tc can be checked in O (1) time. Maintaining the number ofpaths in A allows us to dynamically update G tc upon addition or deletion of( x, y ) ∈ E c using the following algorithm [21]: ∀ u ∈ pred ( x ) , ∀ v ∈ succ ( y ) , A uv ← A uv ± A ux · A yv (1)where pred ( x ) is a set of all direct predecessors of x and succ ( y ) is a set of alldirect successors of y . The ± sign is + for insert and − for delete operation.The instantaneous cycle can be checked before ( x, y ) is inserted to E c (andequivalently E tc ) in O (1) by the following query: x ∈ succ ( y ) (2)For a sparse graph, the runtime complexity of Eq. (1) is O ( | V c | ). For asingle edge update, Eq. (2) gives the fastest stall time. In case when thereare multiple edges to be updated in a single MCR, the oracle has to performEq. (1) as many times as the number of edges that are inserted if it onlymaintains G tc . An alternative method would be inserting or removing alledges requested by the MCR in G c all at once and perform DFS for checkingan instantaneous cycle. Again this would not be amenable if the size of G c is large.To alleviate the time required to check instantaneous cycles in a large G c and to handle multiple edge updates in a single MCR, we perform transi-tive reduction on vertex-partitioned sub-graphs G itc = (cid:104) V ic , E itc (cid:105) of G tc where ∀
15s explained in Section 5.2. Since G itr has the same reachability relation as G itc (which is the same as that of G tc and G c for the vertex partition V ic ) butwith the minimal number of edges, the instantaneous cycle detection on G itr via DFS for some i would be faster than on the entire G c . The oracle adaptively chooses one of G c , G tc and G itr to perform the ICD.Formally, an MCR call made by a model is a sequence of edges E mcr = { e ln } ji =1 = { e l , . . . , e lj } to be updated in G c where l indicates insertion (1)and deletion (0) of an edge. A response to the model γ ∈ { , } from theoracle indicates either an acceptance (1) or rejection (0) for the correspondingMCR. Next, we define three types of MCR calls:1. MCR type a – When a single model makes an MCR with a change ina single input to output dependency, i.e. | E mcr | = 1.2. MCR type b – When an E mcr consists of changes in multiple inputto output dependencies where ∃ =1 G itr ∈ { G tr , ..., G ktr } (there exists a unique G itr in the set) such that ∀ ( x, y ) ∈ E mcr , x, y ∈ V ic .3. MCR type c – When an MCR consists of changes in multiple input tooutput dependencies where ∃ G itr , G jtr ∈ { G tr , ..., G ktr } , i (cid:54) = j, E mcr ∩ E itr (cid:54) = ∅ ∧ E mcr ∩ E jtr (cid:54) = ∅ .The decision on how to divide G tc into a set of sub-graphs for computing G itr depends on applications and the analysis on the frequency of edge updatesin the MCRs. For example, we can group vertices in the same G itr whenall v ∈ V ic appear frequently in the MCRs in the same synchronous tick.Another possibility is to consider each G itr belongs to a single model m ∈ M since all v ∈ V ic would most likely appear in the same synchronous tick. Inthis way, the oracle will receive more MCRs of type b , which are faster toanalyse than MCRs of type c . It should be noted that upon mode change,models only have to provide changes in their input to output dependenciesin E mcr to minimise the amount of data need to be transferred to the oracle.Furthermore, oracle does not need to consider all possible combinations ofmodes from sub-models since only dependency changes in the target modesare needed for the analysis.The pseudocode of the oracle that performs the online ICD is shown inFig. 7. The MCR type a is checked using Eq. (2) as shown at lines 2-3. Forthe MCR type b (line 4), we check if there exist a pre-computed G itr where16 : function DetectCycle ( G c , G tc , { G tr , ..., G ktr } , E mcr )2: if | E mcr | = 1 then (cid:46) MCR type-a3: Perform a cycle check using Eq. (2).4: else if ∃ =1 G itr s.t. ∀ ( x, y ) ∈ E mcr , x, y ∈ V ic then
5: Update G itr with E mcr (cid:46) MCR type-b6: Perform DFS on G itr for presence of any cycles.7: else (cid:46) MCR type-c8: Update G c with E mcr
9: Perform DFS on G c for presence of any cycles.10: end if
11: Send results to models12: if Instantaneous cycle is found then
13: Revert updates on G itr or G c else UpdateGraphs ( G c , G tc , { G tr , ..., G ktr } , E mcr )16: end if end function Figure 7: Pseudocode for the online ICD all the edge updates required by the MCR affect vertices in V ic are in E itr . Ifit is the case, a DFS is performed on G itr for ICD, which is shown at line 6.For the MCR type c (line 7), the edges in E mcr are across multiple G itr . Inthis case, the oracle performs DFS on G c as a fallback option as shown atline 9. Note that both G itr and G c are updated with E mcr before the analysisas shown at lines 5 and 8. Upon detection of an instantaneous cycle, theseupdates are reverted at line 13.After the ICD is done, the oracle updates all three graphs which is indi-cated by the function call UpdateGraphs at line 15 in Fig. 7. The pseu-docode of
UpdateGraphs is shown in Fig. 8 where the lines 9-11 imple-ments Eq. (1) that updates G tc . During the execution of UpdateGraphs ,the oracle may receive a new MCR and we assume it can preempt the execu-tion of
UpdateGraphs and perform ICD for the new MCR. However, theoracle cannot use G tc or G itr for the analysis for the new request because thesegraphs are not up-to-date with the current modes of the models. Instead,DFS on G c is performed as a fallback method, which algorithm is identical tothe lines 8-9 in Fig. 7. In this case, we also assume that the previous invoca-tion of UpdateGraphs is abandoned and it is invoked again after checkinginstantaneous cycle. The new invocation would include all edge updates fromthe previous MCRs that have not yet been processed by
UpdateGraphs .The size of E mcr in UpdateGraphs can grow unbounded if the inter-17 : function UpdateGraphs ( G c , G tc , { G tr , ..., G ktr } , E mcr )2: if | E mcr | > c then RecomputeTC ( G c )4: UpdateTR ( G c , G tc , { G tr , ..., G ktr } )5: return end if A ← adjacency matrix of G tc for ∀ i ∈{ ,..., | E mcr |} ( x, y ) i ∈ E mcr do for all u ∈ pred ( x ) do for all v ∈ succ ( y ) do A uv ← A uv ± A ux · A yv end for end for end for UpdateTR ( G c , G tc , { G tr , ..., G ktr } , E mcr )16: end function Figure 8: Pseudocode for updating G tc function UpdateTR ( G c , G tc , { G tr , ..., G ktr } , E mcr )2: for i ∈ { , ..., k } do A ← adjacency matrix for G ic ⊆ G c ∀ u,v ∈{ , ··· , | V ic |} B uv ← query tc ( u, v )5: C ← A · B ∀ u,v ∈{ , ··· , | V ic |} C uv ← if A uv > ∧ C uv = 0 else
07: Set adjacency matrix for G itr to C end for end function Figure 9: Pseudocode for updating all G itr ∈ { G tr , ..., G ktr } arrival time of MCRs is faster than the speed of the update process of UpdateGraphs . Therefore,
UpdateGraphs invokes
RecomputeTC ifthe size of accumulated E mcr is bigger than certain value c , which is shownat lines 2-6. In other words, if the time required to update G tc from E mcr be-comes longer than recomputing G tc from the scratch, the algorithm invokes RecomputeTC .All G itr are updated in UpdateTR whose pseudocode is shown in Fig. 9.The main part of this algorithm, based on the method introduced in [20],is the multiplication of the adjacency matrices A and B at line 5 for G ic and G itc where G ic = (cid:104) V ic , E ic (cid:105) , E ic = { ( u, v ) | ∀ u, v ∈ V ic , ( u, v ) ∈ E c } . Itshould be noted that G itc is not a transitive closure of G ic . If the adjacency18atrix of G itc is not maintained explicitly, the algorithm has to make querieson G tc to create B as shown at line 4. The result of the multiplication C isfurther processed at line 6 where each element at C ij is set to 1 if its previousvalue was 0 and if the element in A ij at the same index is greater than 0.Otherwise, the value at C ij is set to 0. The adjacency matrix for G itr is thenupdated with C (line 7). The complexity and correctness of the online ICD is presented in thissection. The update time ∆ of a single G itr by the oracle is∆ = δ q + δ tr (3)where δ q is the time required to make | V ic | reachability queries on G tc and δ tr is the actual time required to compute G itr . It is shown in [20] that thetime complexity to compute transitive reduction is in the same class as thetransitive closure computation. Therefore, the fastest known algorithm tocompute G itr is boolean matrix multiplication, and hence O ( δ tr ) = O ( | V ic | ω ).To compute δ q we need to make reachability queries on G tc as many as | V ic | .Yet if we maintain A of G tc explicitly, δ q = 0. Therefore ∆ = O ( | V ic | ω ). Then,runtime complexities of RecomputeTC and
UpdateTR are presented inthe following lemma:
Lemma 1.
The worst-case runtime complexity of the functions
Recom-puteTC and
UpdateTR is O ( | V c | ω ) where ω = 2 . is the greatest lowerbound for the matrix multiplication.Proof. The runtime complexity of transitive closure (
RecomputeTC ) isknown to be in the same class as that of matrix multiplication, which is O ( | V c | ω ). The runtime complexity of UpdateTR is O ( (cid:80) ki =1 | V ic | ω ) where k is the total number of vertex-partitioned sub-graphs G itc and | V ic | ω is fromthe matrix multiplication as shown at lines 2 and 5 in Fig. 9, respectively.From the multinomial and generalised binomial theorems, we can see | V c | ω > (cid:80) ki =1 | V ic | ω where | V c | = (cid:80) ki =1 | V ic | for any i > ω >
1. Therefore, thecomplexity of
UpdateTR is reduced to O ( | V c | ω ).It is obvious that, as shown at lines 2-10 in Fig. 7, the worst-case timebound for the stall time of the models when the oracle was in DetectCycle is dominated by the DFS on G c , which is O ( | V c | ), for the MCR type c . Since19he oracle is able to preempt the graph update process in UpdateGraphs upon receiving new MCRs from the models, the worst-case stall time is also O ( | V c | ). However, note the physical stall times are expected to be muchlower than this worst-case, especially if the oracle is able to frequently exploitthe benefits of MCR types a and b . Time complexity to check if the incomingMCR is type b (line 4 in Fig. 7) is O ( | E mcr | ) where | E mcr | ≤ | E c | since foreach ( x, y ) ∈ E mcr , the oracle can check in which vertex-partitioned V ic x and y are mapped in O (1).The following theorem proves the correctness of our ICD technique: Theorem 1.
The algorithm in Fig. 7 correctly detects instantaneous cyclesfor all three MCR types.Proof.
By definition, G tc has reachability relation between any vertices in V c and the relation is transitive. Then for the MCR type a , it is trivial tocheck if an edge ( x, y ) ∈ E mcr would create a cycle in G c via the query shownin (2) on G tc . Concretely, G tc contains the transitive relation of ( x, y ) in G c where x ∈ succ ( y ) iff there exist a path from y to x in G c . Therefore,one can verify if the addition of an edge ( x, y ) would create a cycle viachecking the existence of a path from y to x via x ∈ succ ( y ). For theMCR type b , since any vertex-partitioned sub-graph G itc has a set of edges E itc = { ( u, v ) | ∀ u, v ∈ V ic , ( u, v ) ∈ E tc } , G itc also has the same reachabilityrelation between any vertices in V ic as G tc and G c and is transitive. Bydefinition, transitive reduction of G itc denoted as G itr has the same reachabilityrelation as G itc and G c for the vertex set V ic . Therefore, for any given G c and G itr , if there exists a set of edges E cycle ⊆ E c that creates a cycle in G c , thenthere also exists the same set of edges that creates a cycle in G itr such that E cycle ⊆ E itr and vice-versa. For the MCR type c , it is trivial to see anyDFS algorithm can find the existence of a cycle in G c , and this concludes theproof.
7. Experiment: Finding the Saturation Point from the Online In-stantaneous Cycle Detection
In this section, we present a set of experiments using synthetically gener-ated examples where the oracle performs the online ICD based on incomingMCR types. All of the benchmarks were written in MATLAB version 2019band carried out on the machine with a Core i7-8650U CPU @ 1.9Ghz laptopwith 16GB RAM. 20 able 1: G c candidates for the experiment G c | V c | | M | | E | | M mcr | | E mcr | G c,
100 5 1274 3 6 G c,
200 10 4943 10 20 G c,
400 10 19991 10 20 G c,
600 10 44880 10 30 G c,
800 10 79900 10 30
The main objective of this experiment is to investigate how fast the mod-els can make MCR before the oracle always falls back to DFS on G c for ICD.We call this threshold a saturation point . Table 1 shows an overview of theexperimental candidates for G c . We have generated a random set of G c inMATLAB based on the parameters shown in the table. The size of ports(vertices) in G c is increased up to 800 to visualise the trend of the saturationpoints for different sizes of G c . As mentioned in Section 6.2, we assume each G itr is a transitive reduction of G itc for a model m i ∈ M . The ports in V c are evenly distributed among the models and the ratio of input to outputports in V ic for each model is set to 1. We use the notation G c,n for indi-cating a composite graph G c with | V c | = n . The number of models | M | isfixed to 10 for all candidates except G c, , which has 5. This is to increasethe complexity of the ICD for MCR type b along with the G c size, i.e. toincrease the size of each G itr . The edges are created based on the uniformlydistributed random numbers such that ∀ u,v ∈ V P (( u, v ) ∈ E ) = 0 . | M mcr | is the maximum number of models and | E mcr | is the number of edges to up-date per MCR, respectively. The types of MCR calls made by the models areevenly distributed in this experiment for each a , b and c and in total 100 callsare made altogether. Finally, we designed this experiment in a way that themodels only send MCR calls that do not create instantaneous cycles. Thisresults in the longest possible stall time for each MCR call since the oraclerequired to search the complete vertices and edges in the graph. The stall time distributions for each experiment are shown in Fig. 10 withvarying MCR call periods denoted as P mcr . The smallest P mcr on the y-axisin each figure indicates the saturation point for the oracle whereas the largest21 · − P m c r ( s ) (a) G c, · − P m c r ( s ) (b) G c, · − P m c r ( s ) (c) G c, · − P m c r ( s ) (d) G c, · − P m c r ( s ) (e) G c, Figure 10: Distributions of stall times for (a) G c, , (b) G c, , (c) G c, , (d) G c, and(e) G c, P mcr is when all of the MCRs are served without falling back to DFS on G c .We say the oracle is saturated when the oracle starts to perform ICD for allMCRs using the fallback method, i.e. all MCRs are treated as type c .For G c, shown in Fig. 10a, the oracle could serve MCRs without sat-uration up to P mcr = 0 .
001 whereas all of the MCRs are served withoutthe fallback method when P mcr > . G c, shown in Fig. 10b, theinterquartile range (IQR) of the box with P mcr = 0 .
005 is relatively shorterthan other boxes. This is due to the small variations between the measuredstall times since most of the MCR calls served by the oracle are done viathe fallback method (DFS) on G c . On the other hand, the IQR becomeslarger with longer P mcr because the oracle completes updating the graphs22 . · − P mcr F a ll b a c k ( % ) (a) G c, · − P mcr F a ll b a c k ( % ) (b) G c, .
14 0 .
16 0 . P mcr F a ll b a c k ( % ) (c) G c, . .
65 0 . . P mcr F a ll b a c k ( % ) (d) G c, . . P mcr F a ll b a c k ( % ) (e) G c, Figure 11: Percentage of MCRs served by the oracle via the DFS fallback method on G c for (a) G c, , (b) G c, , (c) G c, , (d) G c, and (e) G c, more often before the arrival of the next MCR, and therefore more MCRswith types a and b are served by the oracle without the fallback method.The same trend can be seen in the experiments for G c, , G c, and G c, ,which are shown in Figs. 10c, 10d and 10e, respectively. For G c, , P mcr hadto be increased up to 0.18 second to avoid saturation whereas for G c, and G c, this had to be further increased up to 0.75 and 2 seconds, respectively.Note that some of the outliers shown in the figures are due to the cold startof the benchmark program and the measurements are quickly settled downafter a few simulation steps.More details on the proportion of the MCR calls served by the oracle viathe fallback method are shown in Fig. 11. The P mcr values on the x-axis ineach figure correspond to the P mcr values on the y-axis in Fig. 10 of the samesubfigure label. The results indicate that the number of fallback methodstaken by the oracle decreases almost linearly with increasing P mcr values. Oneinteresting outcome is that the variations of the stall times are significantlyreduced when the amount of fallback methods drops to approximately 50%.This can be seen by comparing the results shown in Figures 11 and 10 wherethe lower quartile values are increased such that they are close to medianvalues. This indicates that at this point most of the MCR types a and b ,23 · − . . . . P mcr A cc u m u l a t e d s t a ll t i m e ( s ) (a) G c, . . . · − . . . P mcr A cc u m u l a t e d s t a ll t i m e ( s ) (b) G c, .
13 0 .
14 0 .
15 0 .
16 0 .
17 0 . . . P mcr A cc u m u l a t e d s t a ll t i m e ( s ) (c) G c, . .
63 0 .
65 0 .
68 0 . .
73 0 . P mcr A cc u m u l a t e d s t a ll t i m e ( s ) (d) G c, . . . . . P mcr A cc u m u l a t e d s t a ll t i m e ( s ) (e) G c, Figure 12: Accumulated stall times for (a) G c, , (b) G c, , (c) G c, , (d) G c, and(e) G c, for 100 MCR calls to the oracle which are the fastest to check, are served via the fallback methods.The total accumulated stall times for 100 MCR calls are shown in Fig. 12for different sizes of G c . The results indicate that total stall time is inverselyproportional to P mcr as higher the P mcr value, the number of MCRs servedvia the fallback method decreases. The decrease in the total stall times issmall with smaller sizes of G c . This is expected since the oracle can performDFS on G c fast enough so that the stall time does not differ much from theICD performed on G tc and G itr . On the other hand, the difference with larger G c is significant, for example, it is reduced by 70% for G c, whereas only0.075% for G c, , which are shown in Figs. 12a and 12e, respectively. Thedecreases in the total stall times for other G c are 37% for G c, , 61% for G c, and also 61% for G c, as shown in Figs. 12b, 12c and 12d.The previous experiment was carried out with 100 MCR calls with thetypes a , b and c that are evenly distributed. The effect of different proportionsof the MCR types on the stall time was investigated in the next experimentand the results are shown in Fig. 13. Each of the x, y and z axis representsthe proportion of the MCR types that are sent to the oracle during thebenchmark runs. The colours of the marks indicate the total accumulatedstall times in seconds. This benchmark was run with G c with the settings24 igure 13: Accumulated stall times for models with varying MCR types in 100 calls | V c | = 800, | E | = 79900, | M mcr | = 10 and | E mcr | = 20. It is shown thatthe MCR types a and b did not significantly increase the total stall timeswhereas the time increased up to 5.5 seconds when all of the MCR calls weretype c .Fig. 14 illustrates effects of various partition settings on G tc . The exper-iment was performed on the different sizes of partitions with a fixed MCRtype c to observe how the partitions affect stall times on various G c sizes.The first experiment on G c, in Fig. 14a shows conversions of MCR type c to type b as the number of partitions decreases. With a single large partition,about half of the MCRs were served by the oracle as type b , see the left graphin Fig. 14a. On the other hand, the percentage of the DFS fallback methodincreased with the large partition because the time required to update G itr and G tc would also increase as explained in Section 6.4. We can observe aconstant decrease in the accumulated stall time with increasing size (decreas-ing the number) of partitions with G c, . The same trends can be observedfor G c, and G c, as shown in Fig. 14b and Fig. 14c, respectively. How-ever, it is shown that the accumulated stall times started to increase withlarge partitions for G c, and G c, . This is due to the increased size of the G itr that makes the time required for doing DFS on it becomes longer thanthose smaller G itr in G c, . To decrease stall time, this experiment suggeststhat one needs to partition G tc in a way to convert MCRs of type c to type b as many as possible while avoiding to create too large partitions especially25 P e r ce n t ag e ( % ) Fallback type-c type-b
24 12 6 1012 A cc u m u l a t e d S t a ll t i m e ( s ) (a) G c,
83 40 20 10 1020406080 P e r ce n t ag e ( % ) Fallback type-c type-b
83 40 20 10 10123 A cc u m u l a t e d S t a ll t i m e ( s ) (b) G c,
130 65 34 17 8 4 1 P e r ce n t ag e ( % ) Fallback type-c type-b
130 65 34 17 8 4 1 A cc u m u l a t e d S t a ll t i m e ( s ) (c) G c, Figure 14: A set of experiments with different partition sizes for G tc when the models within these partitions make frequent MCR calls.To apply our approach in a real-time setting, we foresee the use of a tech-nique such as the probabilistic worst-case execution time (pWCET) analy-sis [22] to compute the worst-case stall time bound with certain degrees ofconfidence. In this technique, hardware or software components with hard-to-analyse timing behaviours are subject to statistical randomisation to reducetheir computational complexity. Then, the rare event theories such as ex-treme value theory (EVT) can estimate the probability of maximum of thetiming behaviours. In our case, such complexity arises from the distributionof MCR types on certain choices of partitions as well as their frequency that26ay or may not incur the DFS fallback method. Based on the observed be-haviours of these parameters, we believe pWCET is applicable to the onlineICD in the actual cyber-physical system.
8. Experiment: The Workpiece Sorting System Case Study
This section presents the workpiece sorting system case study that isintroduced earlier in Section 4. The functional requirement of this experimentis to correctly place the workpieces into the bins in the presence of the stalltime that is induced from the ICD. The system consists of twelve pairs ofcontroller and ejector plant models that sort incoming workpieces into one ofthe twelve bins. See Fig. 2, which shows the first three ejectors placed alongthe conveyor belt. Additionally, there is also a controller that controls thespeed of the conveyor belt. The operations performed by the machines uponthe arrival of each workpiece consist of a sequence of mode changes. Our aimis to validate such mode changes via our online ICD approach without failingto sort incoming workpieces into correct bins due to the additional stall timethat would create discrepancy between the model and the physical plant.The oracle is added to include ICD capability for the mode changes in thecontrollers as well as the plant models for the ejectors. The ejector requiresthree mode changes, Idle-to-Push, Push-to-Pull and Pull-to-Idle, that pushthe workpiece into the bin. The corresponding MCRs are sent to the oraclein sequence when the workpiece is detected via the infra-red sensor at theejector.In addition to the system’s operation described in Section 4, we haveextended the scenario with a constraint where the speed of the conveyor beltshould decrease below the threshold of 0.5 unit/s before the ejector is ableto push the workpiece into the bin. Furthermore, the real-time locations ofthe workpieces can only be detected via the infra-red sensors that are placedin line with the ejectors or at the entry point of the conveyor belt.The control of the conveyor belt’s speed is done via the optimisation prob-lem where the controller maximises the speed of the conveyor belt for a finitetime-horizon. The stall time induced by the ICD lags the model’s estima-tion on the physical location of the workpiece. Therefore the conveyor belt’sspeed might not decrease below the threshold before the physical workpiecereaches the ejector, if the stall time is too large. In this case, we consider theejector incorrectly pushes the workpiece and hence fails. The lagged modelscan be re-synchronised with the physical system when the workpiece location27 racle ........ Connectedto the nextejectorsubsystem........ Clear signal forresetting the 'Full' stateToken to be received from the last EjectorSubSystem
EjectorDevice
Clear_inWp_inToken_in Wp_outToken_outControllerModeEjectorMode
EjectorSubSystem1[Z] CB
In1In2Delay Out1wpdetect
Conveyor Belt1UntypedWPGen uwpc t_o
Oracle t_oDelayOutt_oDelayInOrig
In1 Out1
TokenGenerator
Figure 15: Simulink top-level view of the workpiece sorting system is detected by one of the infra-red sensors. The models are developed usingSimulink where the controller for the conveyor belt is implemented usingthe MATLAB script whereas the SimEvent library is used for the conveyor’splant model. Models for the ejector device are implemented using the State-flow diagram as well as the Simulink’s continuous time blocks for capturingthe ejector’s physical movement. We have developed a Simulink-to-compositemodel converter from which G c , G tc and { G tr , ..., G ktr } are generated.An overview of the workpiece sorting system is shown in Fig. 15. Theoracle annotated with a blue rectangle has concatenated input vectors ofMCRs mcr , . . . , mcr n where n is a total number of ejectors which is 12. Anadditional boolean typed input signal wpd indicating the exit of the workpiecefrom the conveyor belt is used for logging purpose. The output of the oracleis the analysis time for an MCR, which is fed to the input of the ConveyorBelt model for simulating the lagged delay. The block EjectorSubSystem1has the controller and plant models for the ejector and the output red linesare connected to the adjacent ejector model, which is omitted due to space . . Simulation time (s) C o n v e y o r s p ee d ( un i t / s ) Figure 16: Change in the conveyor belt speed when the instantaneous cycle detectionperformed only on G c (red line), on G tc (blue line) and G itr and when the stall time iszero (green line) . . Simulation time (s) S t a ll t i m e ( s ) Figure 17: Change in the accumulated stall times due to the instantaneous cycle detection limitations. In this example, the current state of the controller and ejectormodels are used for constructing an MCR call. A token of the last ejectoris returned to the EjectorSubSystem1 via TokenGenerator that implementstoken ring topology for deciding ejection order.To generate G c , G tc and G tr from the Simulink model shown in Fig. 15,we traverse models based on their inter-connections using MATLAB’s built-in Simulink programming model editing APIs and consider every Simulinkmodel libraries as an individual model, i.e. we have (cid:83) ki =1 V ic = V c where k is the total number of model libraries. In addition, we assume every modelexcept the controller and ejector has fixed dependencies between every inputand output ports. Each of the controller and ejector models have four andthree modes, respectively, and we create dynamic input to output dependen-cies based on these modes as described in Section. 4. The generated G c has | V c | = 1409, | E c | = 1441 and | M | = 706.The results of the simulations are shown in Fig. 16. The plot consistsof three simulation runs, indicated by the red, blue and green lines, untilthe completion of the first three ejector operations. The vertical lines arethe instances when, in terms of simulation time, one of the infra-red sensorsdetected the workpiece. The red line that starts from 2 unit/s indicates the29hanges in the conveyor belt speed when the ICD was performed using DFSon G c such as in [5]. The result shows the conveyor speed was correctlydecreased below 0.5 unit/s for the first workpiece since no MCR was issuedprior to this workpiece. However, the estimations of the next two workpiecesarriving at the ejector were incorrectly computed due to the large stall timeinduced from the ICD on G c . As a result, the model simulation significantlylagged behind the physical system’s operation and this is indicated by thevertical red line at around 34.11 s that appeared when the conveyor beltspeed was above 0.5 unit/s threshold. On the other hand, the blue lineshows the result of our approach where the conveyor belt speed was correctlyreduced below the threshold before the workpiece reached the ejector. Wehave partitioned the models in a way that the oracle could utilise G tc or G itr for any MCRs. The ideal case is shown as the green line, when thestall time was zero, for comparison purposes. The workpiece detection times(vertical lines) are different depending on the analysis methods as shown fromthe second ejector operation. This is because the conveyor belt speeds werereduced at different simulation times due to the varying stall times dependingon the ICD methods.The changes in the accumulated stall times for the results presented inFig. 16 are shown in Fig. 17. It is clearly shown that the stall time increasedsignificantly up to 0.434 s when the ICD was done only on G c whereas it isonly up to 0.067 s when our oracle-based approach was used. The modelsare re-synchronised with the physical system when the workpiece is detectedby the sensor. This is indicated by the stall time that drops at around 17.5 sand 57 s. At the same time, the ejector needs to perform the push operationthat again induced additional stall time, and therefore the adjusted stall timedid not drop to zero as shown in the figure.
9. Further Discussion on the Use-Case Scenario
With the development of Functional Mock-up Interface (FMI) standard [1],many modelling tools start to support exporting and importing third-partymodels encapsulated in a standard container called the Functional Mock-upUnit (FMU). FMI standard defines how IO dependencies can be embeddedin FMU. Nevertheless, the standard only supports fixed dependencies for allpossible modes in FMU. For example, [23] demonstrates the use of FMU thatencapsulates thermal models of the more-electric aircraft. This model, whenflattened and consisting of only basic atomic blocks, has 1668 ports ( | V c | )30nd 1492 signal connections ( | E c | ). The model has a delay block in front ofthe FMU in order to break potential instantaneous cycles. Such delay blocksmight not be needed or even result in undesired behaviour when the modelsare executed in the SR semantics and need to be synchronised in the exactSR tick (Section 5.1). Our technique can improve the dependency analysisfor the models that employ FMU especially with the upcoming FMI version(3.0) that supports co-simulation of models with discrete events.Several hybrid system models are presented in [24] with discussions onthe creation of instantaneous cycles. For example, three-point masses movingon a flat surface that collides with each other (discrete reset) during whichinstantaneous cycles can occur in the model. Another example called thefull wave rectifier is presented in [24] where an instantaneous cycle is createdwhen the load to the system is replaced with a pure resistive load. Again,the existence of mode switches in this model makes the static analysis on thecomponent-based models difficult. Alternatively, an online technique can al-leviate this problem and our method targets reducing potential discrepanciesbetween the model and physical system due to the ICD. The three-pointmass example consists of approximately 70 ports and 35 signal connectionswhile 21 ports and 20 signal connections for the rectifier example.In [7], an example of fault tolerant micro-chemical plant was introducedwhere a chemical substance produced from multiple reactants has to be mixedbefore it reaches the final output tank. The fault scenarios that impair mul-tiple pumps and heat exchangers require reconfiguration of mixers so thesubstance can be correctly rerouted to the output tank. In this case, mixerschange their input to output port dependencies that can result in instanta-neous cycle. Therefore, our online approach is also applicable in this sce-nario when the number of combined modes for the system is significantlylarge (10 ) that static analysis infeasible as indicated by the authors. Fur-thermore, when the models are used online to quickly adapt to faults andprovide possible reroutes, reducing stall time from our ICD technique wouldbe beneficial. While the authors presented a simplified version of the plantthat has 2 reactants, 20 valves, 4 pumps and 2 heat exchangers, which resultin approximately 54 ports and 43 IO connections, the size of real system istypically much larger. 31
0. Conclusions and Future Work
This paper addressed the problem of online ICD for the models withmode-dependent input and output dependencies that change during run-time. The method utilises an oracle that determines whether the models canchange their internal modes at each synchronous tick boundary. The methodcan be applied to various modelling formalisms such as the synchronous re-active model of computation if there is a notion of synchronisation point forcommunication with the oracle. The oracle maintains three types of data-structures: a composite graph and its transitive closure as well as transitivereduction, which are adaptively chosen to reduce the stall time of the modelsimulation induced by the online ICD. The experimental results showed theoracle can provide an answer to the mode update queries within an accept-able stall time frame. We also demonstrated the industrial use-case of ourtechnique via the workpiece sorting system. For future work, we would liketo further investigate the impacts of various data-structures for the dynamictransitive closure with varying query and update complexities.
Acknowledgement
This work was supported by Delta-NTU Corporate Lab for Cyber-PhysicalSystems with funding support from Delta Electronics Inc. and the NationalResearch Foundation (NRF) Singapore under the Corp Lab@University Scheme.
References [1] D. Broman, C. Brooks, L. Greenberg, E. A. Lee, M. Masin, S. Tripakis,M. Wetter, Determinate Composition of FMUs for Co-Simulation, in:Proceedings of the Eleventh ACM International Conference on Embed-ded Software, IEEE Press, 2013, p. 2.[2] C. Gomes, C. Thule, D. Broman, P. G. Larsen, H. Vangheluwe, Co-Simulation: A Survey, ACM Computing Surveys (CSUR) 51 (3) (2018)49.[3] R. Lublinerman, C. Szegedy, S. Tripakis, Modular Code Generationfrom Synchronous Block Diagrams: Modularity vs. Code Size, in: ACMSIGPLAN Notices, Vol. 44, ACM, 2009, pp. 78–89.324] E. Matsikoudis, E. A. Lee, The Fixed-Point Theory of Strictly CausalFunctions, Theoretical Computer Science 574 (2015) 39–77.[5] Y. Zhou, E. A. Lee, Causality Interfaces for Actor Networks, ACMTransactions on Embedded Computing Systems (TECS) 7 (3) (2008)29.[6] E. A. Lee, Constructive Models of Discrete and Continuous PhysicalPhenomena, IEEE Access 2 (2014) 797–821.[7] M. W. Hofbaur, F. Wotawa, A Causal Analysis Method for ConcurrentHybrid Automata, in: PROCEEDINGS OF THE NATIONAL CON-FERENCE ON ARTIFICIAL INTELLIGENCE, Vol. 21, Menlo Park,CA; Cambridge, MA; London; AAAI Press; MIT Press; 1999, 2006, p.840.[8] A. Benveniste, B. Caillaud, H. Elmqvist, K. Ghorbal, M. Otter,M. Pouzet, Structural Analysis of Multi-Mode DAE Systems, in: Pro-ceedings of the 20th International Conference on Hybrid Systems: Com-putation and Control, 2017, pp. 253–263.[9] X. Qiu, W. Cen, Z. Qian, Y. Peng, Y. Zhang, X. Lin, J. Zhou, Real-TimeConstrained Cycle Detection in Large Dynamic Graphs, Proceedings ofthe VLDB Endowment 11 (12) (2018) 1876–1888.[10] D. J. Pearce, P. H. Kelly, C. Hankin, Online Cycle Detection and Differ-ence Propagation: Applications to Pointer Analysis, Software QualityJournal 12 (4) (2004) 311–337.[11] D. Lee, M. Kim, A sistributed scheme for dynamic deadlock detectionand resolution, Information Sciences 64 (1-2) (1992) 149–164.[12] C.-H. Huang, Hda: Hierarchical and dependency-aware task mappingfor network-on-chip based embedded systems, Journal of Systems Ar-chitecture 108 (2020) 101740.[13] Y. Kim, S. Choi, J. Jeong, Y. H. Song, Data dependency reduction forhigh-performance fpga implementation of deflate compression algorithm,Journal of Systems Architecture 98 (2019) 41–52.3314] W. Jiang, L. Wen, J. Zhan, K. Jiang, Design optimization ofconfidentiality-critical cyber physical systems with fault detection, Jour-nal of Systems Architecture 107 (2020) 101739.[15] E. A. Lee, H. Zheng, Leveraging Synchronous Language Principles forHeterogeneous Modeling and Design of Embedded Systems, in: Pro-ceedings of the 7th ACM & IEEE international conference on Embeddedsoftware, ACM, 2007, pp. 114–123.[16] G. Berry, G. Gonthier, The ESTEREL Synchronous Programming Lan-guage: Design, Semantics, Implementation, Science of computer pro-gramming 19 (2) (1992) 87–152.[17] F. Le Gall, Powers of Tensors and Fast Matrix Multiplication, in: Pro-ceedings of the 39th international symposium on symbolic and algebraiccomputation, 2014, pp. 296–303.[18] B. M. Kapron, V. King, B. Mountjoy, Dynamic Graph Connectivityin Polylogarithmic Worst Case Time, in: Proceedings of the twenty-fourth annual ACM-SIAM symposium on Discrete algorithms, Societyfor Industrial and Applied Mathematics, 2013, pp. 1131–1142.[19] C. Demetrescu, G. F. Italiano, Trade-Offs for Fully Dynamic TransitiveClosure on DAGs: Breaking Through the O ( n ) barrier, Journal of theACM (JACM) 52 (2) (2005) 147–156.[20] A. V. Aho, M. R. Garey, J. D. Ullman, The Transitive Reduction of aDirected Graph, SIAM Journal on Computing 1 (2) (1972) 131–137.[21] V. King, G. Sagert, A Fully Dynamic Algorithm for Maintaining theTransitive Closure, Journal of Computer and System Sciences 65 (1)(2002) 150–167.[22] F. J. Cazorla, E. Qui˜nones, T. Vardanega, L. Cucu, B. Triquet,G. Bernat, E. Berger, J. Abella, F. Wartel, M. Houston, et al., Proartis:Probabilistically Analyzable Real-Time Systems, ACM Transactions onEmbedded Computing Systems (TECS) 12 (2s) (2013) 1–26.[23] G. Dudgeon, More-electric-aircraft-in-simscape (June 2020).URL https://github.com/mathworks/More-Electric-Aircraft-in-Simscapehttps://github.com/mathworks/More-Electric-Aircraft-in-Simscape