Flexible Nets: a modeling formalism for dynamic systems with uncertain parameters

The modeling of dynamic systems is frequently hampered by a limited knowledge of the system to be modeled and by the difficulty of acquiring accurate data. This often results in a number of uncertain system parameters that are hard to incorporate into a mathematical model. Thus, there is a need for modeling formalisms that can accommodate all available data, even if uncertain, in order to employ them and build useful models. This paper shows how the Flexible Nets (FNs) formalism can be exploited to handle uncertain parameters while offering attractive analysis possibilities. FNs are composed of two nets, an event net and an intensity net, that model the relation between the state and the processes of the system. While the event net captures how the state of the system is updated by the processes in the system, the intensity net models how the speed of such processes is determined by the state of the system. Uncertain parameters are accounted for by sets of inequalities associated with both the event net and the intensity net. FNs are not only demonstrated to be a valuable formalism to cope with system uncertainties, but also to be capable of modeling different system features, such as resource allocation and control actions, in a facile manner.


Introduction
The development of appropriate models is crucial for the design, analysis and control of dynamic systems. The usefulness of a model depends on both its capacity to capture the relevant features of the system, and its capacity for mathematical analysis. These capacities of the model largely rely on the adopted modeling formalism, i.e. on the set of modeling principles and rules that are used to build the model. The task of modeling is often hindered by the lack of detailed system information.
This paper exploits the particular features of Flexible Nets (FNs), a modeling formalism introduced in Júlvez et al. (2018) to study Wilson disease, to model and analyze dynamic systems with uncertain parameters, to account for partially observable systems and, to compute the control actions that optimize a given control objective. Roughly speaking, a dynamic system can be seen as a set of state variables that are modified by means of processes (by process, we mean any event, operation or activity whose occurrence has the potential to change the state of the system). These two basic entities, state and process, are mutually related: on the one hand, the execution of the processes determines how the state changes; on the other hand, the state determines the speed of the processes. These relationships between state and processes are clear, for instance, in a chemical system where the state is given by the amount of molecules and the processes are the reactions taking place in the system. On the one hand, the occurrence of reactions produces a change in the amount of molecules that satisfies the stoichiometry of these reactions. On the other hand, the rate of the reactions depends on the amount of molecules. FNs capture these relationships between state and processes by means of two different nets: the event net and the intensity net.
Let us introduce some of the basic features of FNs by means of a simple chemical reaction network. Assume that the reaction network is composed of the following two reactions: R 1 : ∅ → A R 2 : A → nB + 4nC Reaction R 1 models the production of compound A (each occurrence of R 1 increases the concentration of A, which is denoted [A], one unit), and reaction R 2 models how A is decomposed into compounds B and C. The amounts which [B] and [C] are increased by the occurrence of R 2 depend on n. Let us assume that n is uncertain, but known to be in the interval [20,22]. That is, each occurrence of R 2 decreases [A] one unit, increases [B] n units, and increases [C] 4n units where n ∈ [20,22]. Let us further assume that the initial concentrations of [B] and [C] are 0, and the initial concentration of [A] is known to be in the interval [9.9, 1.1].
The FN in Fig. 1 models the described reaction network. Namely, each chemical compound is a associated with a circle (which will be called a place), and each reaction is associated with a rectangle (which will be called a transition). The stoichiometry of the reactions is modeled by the equations associated with the dots labelled v 1 and v 2 (which will be called event handlers). The uncertain stoichiometry of R 2 is accounted for by the Fig. 1 FN modeling a simple chemical reaction network inequalities associated with v 2 , i.e. a=v, 20v≤b≤22v and c=4b; and the uncertain initial concentration [A] is captured by the inequalities associated with A, i.e. 9.9≤m 0 [A]≤1.1.
The event net of an FN is a graph with three different types of vertices: places, transitions and event handlers. While places and transitions are used to model the state and the processes of the system, respectively, event handlers are used to determine the quantities by which the state is changed when given processes occur. Each place is associated with a state variable, and the value of that variable at a given instant is called marking, or number of tokens, in that place. Similarly, each transition is associated with a process, and the number of times the process has taken place is the number of actions in the transition. Each event handler connects a set of transitions to a set of places, and is associated with a set of linear inequalities that relates actions to marking. Given a number of actions in the connected transitions, any solution of the set of linear inequalities can be used to update the number of tokens in the connected places. Thus, the amount by which the marking changes is allowed to be nondeterministic. This feature of event nets allows the model to account for the different system evolutions that can arise as a consequence of uncertainty in the system. In Fig. 1, the event net is composed of the places, transitions, event handlers and arcs and edges in green.
Let us further assume that the rate of reaction R 1 is uncertain, but constrained to the interval [4.5, 5.5], i.e. the number of reactions that occur per time unit is in [4.5, 5.5], and the rate R 2 satisfies 0.9[A] ≤ rate(R 2 ) ≤ 1.1[A], i.e. it is proportional to [A] with an uncertainty of 10%. These reaction rates are modeled in the FN by the inequality associated with R 1 and the inequality associated with the dot labelled s 1 (which will be denoted intensity handler).
Similarly to event nets, an intensity net is a graph with three different types of vertices: places, transitions and intensity handlers. Places and transitions have the same role as in event nets. Each intensity handler connects a set of places to a set of transitions, and is associated with a set of linear inequalities that relates the number of tokens with the speed of the transitions, i.e. of the processes modeled by the transitions. The intensity, or speed, of a transition determines the rate at which actions are created in the transition. As in event nets, any solution of the inequalities can be used to determine the speed of transition, thus, linear inequalities can be used to model uncertainties in the dynamics of the system. In Fig. 1, the intensity net is composed by the places, transitions, intensity handlers and arcs and edges in blue.
An FN is the result of combining an event net and an intensity net, i.e. an FN is a graph with four types of vertices: places, transitions, event handlers and intensity handlers. While the intensity net establishes the speeds at which actions are generated as a function of the marking, the event net specifies how the generated actions are executed and how a new marking is computed. Thus, although FNs are inspired by Petri nets (Murata 1989), their structure is different, since in addition to places and transitions, FNs have handlers which are connected to places and transitions by arcs and edges. FNs offer both high modeling power and appealing analysis possibilities that aim to make use of all the information provided by the available uncertain parameters. Namely, FNs can accommodate uncertainties in the initial marking, in the marking change produced by the firing of the transitions, in the default speeds of transitions, and in the speed of transitions produced by the marking.
A number of different formalisms can be found in the literature that can, to some extend, incorporate uncertain parameters in their models. Depending on the domain of their state variables, these formalisms can be roughly classified as those whose variables are integer numbers, and those whose variables are real numbers (hybrid approaches combine both types of variables).
In the discrete domain, extensions of the most popular modeling formalisms exist that can handle uncertain parameters. For instance, in the Petri nets arena, uncertain knowledge of the marking is particularly well handled by fuzzy (Looney 1988) and possibilistic Petri nets (Cardoso et al. 1999), uncertainty in the firing of transitions can be accounted for by labeled (Cabasino et al. 2010) and interpreted (Ramirez-Trevino et al. 2003) Petri nets, and uncertain firing times can be modeled by time Petri nets (Merlin and Faber 1976) and stochastic Petri nets (Ajmone Marsan et al. 1995). Other discrete formalisms that include extensions to account for uncertain parameters are probabilistic Boolean networks (Shmulevich et al. 2002), which are an extension of Boolean networks (Wang et al. 2012), and influence diagrams (Detwarasiti and Shachter 2005), which are generalizations of Bayesian networks (Needham et al. 2007) that can solve decision problems under uncertainty. Stochastic extensions of timed automata (Alur and Dill 1994) also exist in which delays and discrete choices are made randomly (Bertrand et al. 2014). In a similar vein, stochastic extensions of process algebras (Priami et al. 2001;Ciocchetta and Hillston 2008) have been proposed to describe components with uncertain behaviour (Clark et al. 2007).
A major difference of the above mentioned approaches with respect to FNs is that the state variables of FNs are real numbers. This implies that genuinely discrete systems cannot be modeled by FNs. Nevertheless, the use of real variables facilitates, in general, the use of more efficient analysis methods since the state explosion problem inherent to large discrete systems is avoided, and linear programming techniques can be applied. Moreover, large discrete populations can be approximated reasonably well in many cases by means of real variables (Bortolussi et al. 2013). With respect to the existing stochastic approaches and extensions, it should be said that they offer the possibility to perform useful statistical analyses, which usually require information about the probability distributions of the system, and often involve a significant computational cost. In contrast, FNs do not require information about probability distributions, just about the intervals in which the uncertain parameters lay. This results in efficient analysis techniques based on linear programming.
The most popular modeling approaches in the continuous domain are based on differential equations (Braun et al. 1983;Tyson et al. 2011;van den Berg et al. 2008). In particular, the relaxation of the integrality constraint in a discrete formalism usually leads to models that are governed by differential equations. For instance, the evolution of continuous Petri nets (Silva et al. 2011), which can be seen as a relaxation of Petri nets (Murata 1989), is determined by a set of ordinary differential equations. Another popular modeling formalism that can graphically represent systems in different domains and that can be easily converted to state space representation is bond graphs (Borutzky 2010). However, it should be emphasised that a potential difficulty in the design of models based on differential equations is that exhaustive and accurate information about the system dynamics is required, i.e. uncertain parameters cannot be easily handled. Note that the time trajectory of a system modeled by ordinary differential equations is continuous and deterministic. On the other hand, constraint-based models, which are popular in systems biology (Varma and Palsson 1994;Orth et al. 2011), can incorporate uncertain dynamic information but their analysis capabilities are limited to the steady state.
An important feature of FNs is that they can bridge the gap between deterministic models and constraint-based models by allowing the incorporation of uncertain parameters. Thus, on the one hand and similarly to continuous Petri nets which are deterministic (Jiménez et al. 2004), FNs can model positive linear systems; on the other hand and similarly to constraint-based models (Varma and Palsson 1994), FNs can model systems with uncertain initial state and uncertain process speeds. In contrast to Petri nets, the timing of transitions and the marking changes in places are explicitly separated in FNs: the timing is handled by the intensity handlers, and the marking changes by the event handlers. This can lead to a clearer and more concise graphical representation of the system. Moreover, efficient computational methods exist to analyse the transient state of FNs.
The rest of the paper is organized as follows: Section 2 introduces event nets and shows how a partially observable system can be modeled. Intensity nets are presented in Section 3. The combination of these nets leads to FNs, which are defined in Section 4. Section 5 shows how FNs can handle systems with uncertain parameters and analyze them. Section 6 concludes the paper.

Event nets
In the following, the reader is assumed to be familiar with Petri nets (see Murata 1989 for a gentle introduction).

Definition and state equations
This section introduces event nets, which can be denoted as T V P nets, i.e. actions in transitions T produce and consume tokens in places P through event handlers V . Event handlers connect places and transitions, and determine the marking changes according to the actions in the transitions. In contrast to Petri nets, the net elements that produce changes in the marking are the event handlers, and such changes are allowed to be nondeterministic. The places, depicted as circles, model the different types of components or elements in the system, e.g. resources, products, items, etc. The transitions, depicted as rectangles, model the different types of operations, activities or processes in the system. Such operations require time to be performed and have the potential to change, i.e. produce and consume, the amount of components, i.e. the marking. The event handlers, depicted as dots, model the different ways in which the transitions can change the marking.

Definition 1 (Event net) An event net is a tuple
The vertices of the net are connected by the edges in E V . Each pair of vertices can be connected by at most one edge. The set E V is partitioned into two sets E P V and E T V , where E P V is a set of directed edges connecting places to event handlers and vice versa, and E T V is a set of undirected edges connecting transitions and event handlers. For simplicity, directed edges are referred as arcs, and undirected edges as edges. More formally: -Every e ∈ E P V is either an arc e = (p i , v k ) from a place p i to a handler v k , or an arc e = (v k , p i ) from a handler v k to a place p i .
-Every e ∈ E T V is an edge e = {t j , v k } connecting a transition t j and a handler v k . Notice that direct connections among places and transitions are not allowed. The following notation is used: Example 1 The event net in Fig. 2a has three places, P = {p 1 , p 2 , p 3 }, one transition, T = {t 1 }, and two event handlers The set of all arcs and edges is As examples for the introduced notation, the set of output handlers of p 1 is p v 1 = {v 1 , v 2 }, the set of transitions connected to v 2 is t v 2 = {t 1 }, and the set of input handlers of p 3 is v p 3 = {v 2 }.
In an event net, each place contains a number of tokens (or marking), and each transition contains a number of actions that represent the potential of the system to carry out the associated process. In contrast to tokens, actions require time to be produced (the production rate of actions is determined by the intensity net, see Section 3). The state of an event net accounts not only for the marking and the number of actions, but also for the marking changes and the execution of actions: Definition 2 (State) The state of an event net N V is given by the tuple (σ, a T , a E , m, m), where: is the number of actions of t j executed by v k .
Thus, every net element (except for the event handlers) is associated with at least one nonnegative real variable. Since actions need time to be produced, at the initial state it holds σ = 0, a T = 0 and a E = 0. Notice that dealing with real state variables instead of discrete ones allows the model to incorporate real quantities and to approximate large discrete quantities as in continuous Petri nets (Silva et al. 2011). In any case, the state variables can be constrained to the nonnegative integers if required, see Section 2.2.
Each event handler v k ∈ V is associated with a set of linear inequalities that relate the number of actions executed in the connected transitions to the marking changes in the connected places. The coefficients of such a set of inequalities can be expressed by two matrices (A k , B k ) of real numbers and the same number of rows that are associated with each handler v k ∈ V . The number of actions executed by v k , a f , and the produced marking changes, m f , is given by A k m f ≤B k a f . The columns of , where 'a', 'b' and 'x' are used to label arcs and edges. More precisely, in v 1 :a=b, 'a' denotes the number of tokens in p 1 consumed by v 1 , and 'b' denotes the number of tokens in p 2 produced by v 1 . In the inequalities associated with v 2 , 'a' and 'b' denote the number of tokens in p 1 consumed by v 2 and the number of tokens in p 3 in produced by v 2 respectively, and 'x' denotes the number of actions in t 1 executed by v 2 . The equality v 1 :a=b means that the number of tokens in p 1 consumed by v 1 is equal to the number of tokens in p 2 produced by v 1 . In other words, for every token in p 1 consumed by v 1 , a token is produced in p 2 by v 1 , this can be interpreted as tokens moving from p 1 to p 2 through v 1 .
The equation a=x associated with v 2 means that the number of tokens in p 1 consumed by v 2 is equal to the number of actions in t 1 executed by v 2 , e.g. if one action is executed then one token is consumed. Moreover, the inequality x≤b≤2x means that the execution of one action in t 1 by v 2 , i.e. x = 1, produces a nondeterministic quantity b ∈ [1, 2] of tokens in p 3 (each execution of an action can produce a different amount b ∈ [1, 2] of tokens in p 3 ).
Notice that v 1 is not connected to any transition. This means that no process is required to move a token from p 1 to p 2 . For the sake of mathematical notation, it can be assumed that v 1 is connected to a fake transition, t f ake , that has no effect on the model. The matrices A 1 , B 1 , A 2 and B 2 that capture the inequalities associated with the event handlers are: where the indices of the columns of A 1 are ordered as (p 1 , v 1 ), (v 1 , p 2 ); the index of the column of B 1 is {t f ake , v 1 }; the indices of the columns of A 2 are ordered as (p 1 , v 2 ), (v 2 , p 3 ); and the index of B 2 is {t 1 , v 2 }. Thus, the number of actions executed and marking changes produced by v 2 are related by: Matrices A 1 and A 2 (B 1 and B 2 ) can be arranged diagonally to obtain A(B): Notice that the number of actions produced in a transition (by the intensity net), σ , is equal to the number of actions that have been executed by the connected event handlers, a E , plus the number of actions still available, a T , hence, it holds: Similarly, the number of tokens in a place p i is equal to the initial number of tokens, which is denoted m 0 [p i ], minus the number of tokens consumed plus the number of tokens produced by the connected event handlers: The event net establishes how the state evolves as event handlers are enabled and fire.
indexed by the arcs of v k exist such that: Inequality (3) guarantees that enough actions are available, (4) makes use of the matrices A k and B k to relate the number of executed actions to the marking changes, (5) guarantees that enough tokens are available in the input places to be consumed, (6) guarantees that the overall state change is not null. Notice that the inequalities (4) allow the modeling of uncertainty in the marking changes produced by the execution of actions.
Definition 4 (Firing) An event handler v k enabled at (σ, a T , a E , m, m) can fire. The firing of v k leads instantaneously to a new state (σ, a T , a E , m , m ) where only the variables associated with edges, arcs, places and transitions connected to v k are updated as follows: where a f and m f satisfy (3), (4), (5) and (6).
Notice that an enabled handler is not forced to fire, and that the state reached by the firing of an event handler is allowed to be nondeterministic (see inequality (4)). Moreover, the firing does not force the execution of a minimum number of actions nor the consumption or production of a minimum number of tokens. In fact, the equations in Definition 4 are trivially satisfied with a f = 0 and m f = 0. Thus, such equations also hold for every non-enabled handler with a f = 0 and m f = 0.
The overall change in the state produced by several firings is the result of adding the changes produced by each firing. This leads to a set of equations that are satisfied by the states that can be reached from the initial state.

Proposition 1 (State equations) Let the state of an event net
σ actions are available and no event handler has fired. Every state (σ, a T , a E , m, m) reachable from (σ, σ, 0, 0, m 0 ) belongs to SE N V (σ, m 0 ) where: where Y σ and Z m are determined by the net structure: Proof Let us show that equations (7) necessarily hold for every state (σ, a T , a E , m, m) (7) states that the number of actions produced, σ , is equal to the number of actions executed, a E , plus the number of actions available, a T . This is equivalent to (1) expressed in matrix form, and thus necessarily holds. Equation A m ≤ Ba E in (7) is the matrix form of (4) and accounts for all the actions executed and all the marking changes produced by the firing of all the event handlers. That is, it captures all the cumulative marking changes, m, produced by all the executed actions, a E , and, hence, must necessarily hold. Finally, equation m = m 0 + Z m m in (7) which updates the number of tokens, m, in all the places according to the cumulative marking changes, m, is the matrix form of (2) and, hence, must also hold.
Roughly speaking, the role of matrix Y σ is to distribute the actions in transitions among the handlers connected to them, see (1). The role of Z m is to collect and add the marking changes produced by the firings, see (2). Fig. 2a assuming again that v 1 is connected to a fake transition is: Let us assume that the marking of the net in Fig. 2a and no event handler has fired. This corresponds to the state ( The graph in Fig. 2b shows the potential evolutions of the net under the assumption that event handlers fire in discrete amounts, i.e. all markings and actions are integers. The components of the vectors of states in the graph correspond to the variables The arcs are labeled with the event handler that is fired. Remark that while the firing of v 1 produces a deterministic change (one token consumed from p 1 and one token produced in p 2 ), the state change produced by the firing of v 2 is nondeterministic (either one or two tokens can be produced in p 3 ).
Notice that equations (7) account for the cumulative effect, and not the sequence, of the firings. In particular, the availability of tokens and actions consumed by the sequence of firings is not checked. This can lead to spurious solutions (Silva et al. 1998) in the state equations. Hence, equations (7) represent a necessary condition for the reachability of (σ, a T , a E , m, m).
In order to account for linear relationships among the values of the initial marking, m 0 is assumed to be a vector constrained as: where J m and K m are real matrices of appropriate size. Note that the inequalities (8) (7) can be easily modified to take into account the relationships expressed by (8): Although event handlers are not forced to fire, it is useful in some cases to consider only those states in which all the actions of given transitions have been executed. Let T F ⊆ T be the set of transitions whose actions must have been executed, i.e. the number of available actions of t j ∈ T F must be 0. In order to constrain (9) to such a set of states, the following equation can be added:

Partial observability
In an event net, the marking change produced by the execution of an action is allowed to be nondeterministic. This is the case when there are several event handlers connected to a transition, or if the connected event handler has appropriate inequalities associated with it. The nondeterministic effect of the execution of actions can be used to develop nondeterministic models and, in particular, to model partially observable systems. The state equations (9) account for all the states that can be reached after the firing of event handlers. In the context of partial observability, these equations characterize the set of states that are consistent with a given observation of a partially observable system. The following example shows how a partially observable system can be modeled by an event net. For the sake of this example, let us assume that the event handlers of the net in Fig. 3a only fire in discrete amounts and that only the sequence of transitions whose actions are executed is observable. For clarity, the equations of event handlers that make all their labels equal are omitted, e.g. the equations of v 1 are a=b and a=x, and they are therefore omitted (this same omission is made in the nets of the rest of the paper).
Thus, the observation of t 1 corresponds to the firing of v 1 ; the firing of v 2 is unobservable (silent event) because it is not connected to a transition; the observation of t 2 corresponds either to the firing of v 3 or v 4 because both handlers use the actions in t 2 Fig. 3 a Event net modeling a partially observable system. b Potential evolutions of the marking with discrete firings (in other words, t 2 models events that cannot be distinguished by an observer); the observation of t 3 corresponds to the firing of v 5 .
The graph in Fig. 3b shows the potential evolutions of the net with m 0 = (1 0 0 0 0 0) and σ = (2 1 1), the components of the vectors in the graph correspond to the marking of places (p 1 p 2 p 3 p 4 p 5 p 6 ). Initially, only v 1 can fire and, hence, the execution of actions in t 1 is the only event that can be observed. The observation of t 1 means that a token was consumed from p 1 and a token was produced in p 2 . Since p 2 is the input place of v 2 whose firing cannot be observed, a token in p 2 can remain in p 2 or move silently to p 3 . Thus, the set of markings consistent with the observation of t 1 is {(0 1 0 0 0 0), (0 0 1 0 0 0)}, see second row of the graph. Assume that t 2 is now observed, i.e. either v 3 or v 4 has fired. If the marking was (0 0 1 0 0 0), then v 3 (the only enabled handler) has fired and that leads to marking (0 0 0 1 0 0). Otherwise, v 4 has fired and that leads either to (0 0 0 0 1 0) or (0 0 0 0 2 0) because one firing of v 4 can produce one or two tokens in p 5 . Thus, the set of markings consistent with the observation of the sequence of events t 1 and t 2 is {(0 0 0 0 1 0), (0 0 0 0 2 0), (0 0 0 1 0 0)}, see third row of the graph. Assume that t 3 is now observed, i.e. v 5 has fired. Notice that the firing of v 5 requires two tokens in p 5 , thus, after the observation of t 1 and t 2 , (0 0 0 0 2 0) is the only consistent marking at which v 5 can fire. Consequently, the only sequence of markings consistent with the observation of t 1 , t 2 and t 3 is (1 0 0 0 0 0), (0 1 0 0 0 0), (0 0 0 0 2 0) and (0 0 0 0 0 1).
Notice that the state equations (7) can be used straightforwardly to compute the set of consistent markings with a given observation. For instance, for the observation of t 1 , t 2 and t 3 discussed above, the marking m = (0 0 0 0 0 1) is the only solution of: where σ = (2 1 1) and a T = (1 0 0) are the number of actions available at the beginning and at the end respectively. That is, each transition was observed once, i.e. σ −a T = (1 1 1).

Definition and state equations
This section introduces intensity nets, which can be denoted as P ST nets, i.e. tokens in places P produce and consume intensities in transitions T through intensity handlers S. The intensity of a transition t j is the speed at which actions are produced in t j . In other words, the number of actions produced at t is given by the integral over time of the intensity of t j . Intensity nets and event nets operate in a similar fashion. In fact, the changes in the intensities are produced by tokens in the intensity net in the same way that changes in the marking are produced by actions in the event net. The vertices of the net are connected by the edges in E S . Each pair of vertices can be connected by at most one edge. The set E S is partitioned into two sets E T S and E P S , where E T S is a set of directed edges (or simply arcs) connecting transitions to intensity handlers and vice versa, and E P S is a set of undirected edges (or simply edges) connecting places and intensity handlers. Thus, although both event handlers and intensity handlers are represented as dots, they can be easily distinguished by the arcs and edges that connect them to transitions and places. More formally: -Every e ∈ E T S is either an arc e = (t j , s l ) from a transition t j to a handler s l , or an arc e = (s l , t j ) from a handler s l to a transition t j .
-Every e ∈ E P S is an edge e = {p i , s l } connecting a place p i and a handler s l . As in the event net, connections among places and transitions are not allowed. The following notation is used: t s l denotes the input transitions of s l , i.e. t s l = {t j |(t j , s l ) ∈ E T S } -s t l denotes the output transitions of s l , i.e. s t l = {t j |(s l , t j ) ∈ E T S } -s t j denotes the input handlers of t j , i.e. s t j = {s l |(s l , t j ) ∈ E T S } -t s j denotes the output handlers of t j , i.e. t s j = {s l |(t j , s l ) ∈ E T S } -p s l denotes the places connected to s l , i.e. p s l = {p i |{p i , s l } ∈ E P S } -p s i denotes the handlers connected to p i , i.e. p s i = {s l |{p i , s l } ∈ E P S } As in the event net, the places in the intensity net contain tokens. These tokens can be used by the intensity handlers to produce intensities. A token is active if it is being used by an intensity handler, otherwise it is idle. While idle tokens are associated with places, active tokens are associated with edges. An intensity handler determines how much intensity is produced in its arcs as a function of the number of active tokens in its edges.
The state of the net is given by the variables associated with the net elements. Formally: Definition 6 (State) The state of an intensity net N S is given by the tuple (m, μ P , μ E , λ, λ), where:

a vector indexed by P where m[p i ] is the number of tokens in
is the number of active tokens of p i being used by s l ,

≥0 is a vector indexed by E T S where λ[(t j , s l )] is a decrease of intensity in t j produced by s l , and λ[(s l , t j )] is an increase of intensity in t j produced by s l
The number of tokens in a place p i is equal to the number of its idle tokens plus the number of its active tokens: An intensity handler is said to be working when it is producing intensities. When an intensity handler s l ∈ S starts working, the number of idle tokens in p s l decreases, the number of active tokens in its edges increases (such tokens start being used by the handler), and intensities are produced in its arcs. Conversely, when an intensity handler s l stops working, the number of idle tokens in p s l increases (i.e. they are released by the handler), the number of active tokens becomes 0, and no intensities are produced in its arcs. Thus (in contrast to the firing of event handlers whose firing cannot be reversed once it has occurred) intensity handlers are allowed to start and stop working (or to increase and decrease their working rates), thus allocating tokens as active tokens and releasing them as idle tokens over time.
The relation between the number of active tokens and the intensities produced are given by a set of inequalities associated with each intensity handler s l ∈ S. The coefficients of these inequalities can be captured by two matrices (C l , D l ) of real numbers and same number of rows. The columns of C l are indexed by the arcs connecting s l to transitions. The columns of D l are indexed by the edges connecting places to s l . Matrix C(D) is obtained by arranging all the matrices C l (D l ) diagonally.
Each transition t j is assigned a default (or nominal) intensity λ 0 [t j ]. Thus, the intensity λ[t j ] in a transition t j is equal to λ 0 [t j ] plus the positive changes in intensity minus the negative changes in intensity: The number of active tokens, μ w ∈ R | p s l | ≥0 indexed by the edges of s l , being used by an intensity handler s l , and the intensities, λ w ∈ R | t s l |+|s t l | ≥0 indexed by the arcs of s l , produced by s l are related by the matrices C l and D l as follows: If 1μ w + 1 λ w > 0, then s l is said to be working. Similarly to event handlers, intensity handlers are not forced to work. When a number of intensity handlers work simultaneously, they share the tokens in places and collaborate in the production of intensities. In a similar way to (4), the inequalities (13) allow the modeling of uncertainty in the intensity changes produced by the active tokens.
Similarly to (7), the state equations of the intensity net determine the potential states of the net for a given marking m and default intensities λ 0 : N S be (m, m, 0, 0, λ 0 ), i.e. m idle tokens are available and no intensity handler is working. Every state (m, μ P , μ E , λ, λ) reachable from (m, m, 0, 0, λ 0 ) belongs to SE N S (m, λ 0 ) where:

Proposition 2 (State equations) Let the state of an intensity net
where Y m and Z λ are matrices determined by the net structure: -Y m is a matrix with rows indexed by P , columns indexed by E P S , and such that Y m [p i , {p i , s l }] = 1 ∀{p i , s l } ∈ E P S and the rest of the elements in Y m are 0, -Z λ is a matrix with rows indexed by T , columns indexed by E T S , and such that Z λ [t j , (t j , s l )] = −1 ∀(t j , s l ) ∈ E T S , Z λ [t j , (s l , t j )] = 1 ∀(s l , t j ) ∈ E T S and the rest of the elements in Z λ are 0, and μ P , μ E , λ and λ are nonnegative variables.
As in (8), inequalities can be considered to model linear relationships among the values of the default intensities, λ 0 . Let us assume that λ 0 is a vector constrained as: where J λ and K λ are real matrices of appropriate size. Similarly to (8), the inequalities (15) can be used to model the uncertain default intensity of a transition, or to establish linear constraints among default intensities of transitions. Equations (14) can be easily modified to take into account the relationships expressed by (15): Although intensity handlers are not forced to work, it is useful in some cases to consider only those states in which all the tokens of given places are active. Let P F ⊆ P be the set of places whose tokens must be active, i.e. the number of idle tokens of p i ∈ P F must be 0. In order to constrain (16) to such a set of states, the following equation can be added: Figure 4 shows some of the modeling capabilities of intensity nets. The default intensity, λ 0 [t], of a transition, t, can be written next to the transition, see Fig. 4b (default intensities equal to 0 are omitted in the figures). As in the event nets, labels are associated with arcs and edges to represent amounts of produced/consumed intensities and number of active tokens. The intensity net in Fig. 4a has one place p 1 , one transition t 1 and one intensity handler s 1 . The inequality associated with s 1 establishes that the intensity λ[(s 1 , t 1 )] produced in the arc (s 1 , t 1 ) by s 1 must satisfy 2μ

Modeling capabilities
is the number of active tokens in {p 1 , s 1 } (notice that given that m[p 1 ]=2, the number of active tokens is upper bounded by 2). The actual value of λ[(s 1 , t 1 )] is selected nondeterministically in this interval. Since the default intensity of t 1 is 0 and (s 1 , t 1 ) is the only arc connected to t 1 , it holds λ[t 1 ]= λ [(s 1 , t 1 )] what establishes the rate at which actions will be produced in t 1 .
The intensity handler s 1 in Fig. 4b makes use of the active tokens in p 1 to decrease the intensity in t 1 and increase the intensity in t 2 . This can be seen as an intensity transfer from one transition to the other. According to the equations associated with s 1 , which can be rewritten as μ E [{p 1 , s 1 }] = λ[(t 1 , s 1 )] = λ[(s 1 , t 2 )], the amount of this transfer is The net in Fig. 4c shows how the intensity of one transition, t 1 , can be used to produce intensity in other transitions, t 2 and t 3 . For this net, the state equations (14) become Given that λ 0 [t 1 ] = 6 and λ 0 [t 2 ] = λ 0 [t 3 ] = 0, these state equations reduce to λ[t 1 ]+λ[t 2 ]+λ[t 3 ] = 6 what summarizes the potential intensities of the net.
The net in Fig. 4d models a choice in place p 1 , i.e. each token in p 1 can be used either to produce an intensity within the interval [1, 2] in t 1 , or synchronize with a token in p 2 to produce an intensity within the interval [3, 5] in t 2 .

Flexible nets
This section introduces Flexible Nets (FNs), which can be denoted as P H T nets, i.e. places P and transitions T are connected by event and intensity handlers. Roughly, an FN consists of an event net and an intensity net that have the same set of places and the same set of transitions. In an FN, the event net determines the way actions produce marking changes, and the intensity net determines the way tokens produce intensity changes. The inequalities associated with handlers allow the modeler to cover a range of relationships between "actions and tokens" and "tokens and intensities". Thus, handlers can be seen as a flexible layer between places and transitions that offers the possibility to model uncertainties in both the way actions produce marking changes, and the way tokens produce intensity changes.
The FN in Fig. 5 is composed of the event net in Fig. 2a and the intensity net in Fig. 4a. While the event net determines the marking changes produced by the firing of event handlers, the intensity net establishes the rate at which actions are created in t 1 . Notice that the firing of v 2 implies the execution of actions in t 1 , i.e. actions need to be produced in t 1 so that v 2 can fire. On the other hand, v 1 is not connected to any transitions and, thus, it can fire when there is a positive marking in p 1 . It should be noted that this is not equivalent to an immediate transition in Petri nets, since the firing of t 2 is not forced to happen as soon as the marking of p 1 is positive, its firing can occur at any time at which the marking of p 1 is positive.
In order to compute the number of actions produced in transitions, the number of actions produced in the intensity arcs will be computed first. Let σ (τ ) denote the number of actions produced in the intensity arcs until time τ ( σ [e](τ ) with e ∈ E T S denotes the number of actions produced in e). The value of σ (τ ) is defined as the integral of λ over time: The overall number of actions, σ [t j ](τ ), produced in a transition t j can be computed by integrating λ[t j ], or equivalently, by making use of Z λ , see (14), and σ (τ ): In addition to the state variables of the event and intensity net, σ is included in the tuple of variables defining the state of the FN.

Definition 8 (State)
The state x of an FN is given by the tuple x = (m, μ P , μ E , λ, λ, σ, σ, a T , a E , m). All the state variables are time dependent. For the sake of clarity, the time dependency will be omitted when it is clear from the context, e.g. m(τ ) is shortened to m. At time 0 it holds σ = 0, σ = 0, a T = 0, a E = 0, m = 0, i.e. the initial state can be written as: (m, μ P , μ E , λ, λ, 0, 0, 0, 0, 0).
By making use of SE N V (σ, J m , K m ) in (9), SE N S (m, J λ , K λ ) in (16), (18) and (19), it is possible to write a set of equations that any potential state at time τ must satisfy.
Proposition 3 (State equations) Let N be an FN with initial marking m 0 satisfying J m m 0 ≤ K m , and default intensities λ 0 satisfying J λ λ 0 ≤ K λ . Every state (m, μ P , μ E , λ, λ, σ, σ, a T , a E , m) reachable at time τ belongs to SE N (τ, J m , K m , J λ , K λ ) where: where every variable is nonnegative.
This way, an FN is a continuous time model where time, denoted as τ , is the independent variable and all the state variables are nonnegative reals.
Equations (20) can be interpreted as follows: at a given time τ , some of the produced actions (σ ) are available (a T ), and the rest (a E ) were executed before τ . The executed actions produced marking changes ( m) which resulted in the marking m in places at τ . Some of the tokens in m are active (μ E ) and the rest are idle (μ P ). Active tokens produce intensity changes ( λ) which result in overall intensities (λ) in transitions at τ . The integral of the intensity changes and overall intensities over time after τ will produce more actions (σ ), i.e. σ is produced as time elapses. This behavior repeats over time: when a new marking is reached, intensities are updated, which can lead to the production and execution of new actions, which consequently results in a new marking.
As in the event and intensity nets, constraints (10) and (17) can be added to (20) to force the execution of actions and the activity of places.

Exploiting uncertainty
The state equations (20) account for all the potential states of the net at time τ . In order to facilitate the analysis of FNs, a set of necessary reachability conditions was developed (Júlvez et al. 2018). These conditions consist of linear and quadratic inequalities that all the solutions of (20) must satisfy during the interval [0, τ ]. In order to obtain a time trajectory of the state, i.e. values of the state at different time instants τ 1 , τ 2 , τ 3 , . . ., two methods are considered: 1. The first method consists of developing a unique set of necessary reachability conditions that combines the reachability conditions of each interval [0, τ 1 ], [τ 1 , τ 2 ], [τ 2 , τ 3 ], . . ., in such a way that all the states that satisfy the constraints at the end of a given interval are taken as potential initial states for the next interval (see Júlvez et al. 2018 for details). Once this set of constraints is obtained, a particular trajectory of the FN can be computed by adding an objective function to such a set of constraints, and by solving the resulting programming problem.
It should be noted that solving convex quadratic programming problems is required by both methods above. Given that the computational complexity required to solve such problems is polynomial, the proposed computational methods can be applied to large FNs. The following subsections present some of the modeling, analysis and control capabilities of FNs by modeling a linear system with uncertain parameters, a resource allocation system and a system with control actions.

Linear system with uncertain parameters
The FN in Fig. 6 models a linear system with uncertain dynamics. More precisely, if we assume that all the tokens are forced to be active and all the actions to be executed, the rate at which the marking changes can be expressed as: where q and h are uncertain parameters but known to be in the intervals: q ∈ [1.0, 1.5], h ∈ [0.9, 1.1]. Thus, any potential time trajectory of the system will satisfy (21) with values of q and h within the given intervals. The uncertain parameter q is modeled by the default intensity of t 1 , and h is modeled by the inequalities of s 3 . Notice that the FN in Fig. 6 combines transitions with constant speed, e.g. t 4 , transitions whose speed is proportional to the marking of a place, e.g. t 2 , and uncertain parameters. The trajectory in Fig. 7a is obtained by the objective function "min m[p 3 ]", i.e. the goal is to minimize the marking of p 3 at the end of each interval. In the plots,λ[t j ] denotes the average intensity of t j during each interval. For such an objective, the solution of the programming problem sets the uncertain parameters to λ 0 [t 1 ]=1.5 and 0.9a=r (which results in λ[t 3 ]=0.9m[p 1 ]). This setting minimizes the flow directed from p 1 to the branch composed of v 3 and v 4 . As expected, the intensity of t 1 and t 2 (t 3 and t 4 ) is the same at steady state.
The trajectory in Fig. 7b is obtained by the objective function "max m[p 3 ]". For such an objective, the solution of the programming problem sets λ 0 [t 1 ]=1.0 and r=1.1a (which results in λ[t 3 ] = 1.1m[p 1 ]). This setting maximizes the flow directed from p 1 to the branch composed of v 3 and v 4 .

Resource allocation
The FNs in Fig. 8 models a dynamic system in which shared resources can be allocated to different production lines. Such a net shows how the tokens of a given place can activate different processes (those places have several intensity edges) and can cooperate with active tokens of other places. Namely, there are two types of resources, p a and p b , and three production lines, t 1 , t 2 and t 3 . The production line associated with t 1 (t 2 ) uses the raw material modeled by the tokens in p 1 (p 3 ) and produces items modeled by the tokens in p 2 (p 4 ). The production line associated with t 3 produces tokens in p 5 and it is assumed that it requires no raw material (or equivalently, this raw material is inexhaustible). In order to operate, the production line associated with t 1 (t 3 ) requires the allocation of resources of type p a (p b ). The speed of these production lines, t 1 and t 3 , is proportional to the number of tokens allocated to them. The operation of production line t 2 requires the cooperation of both resources, p a and p b , i.e. tokens of both resources must synchronize in equal amounts to make t 2 work. The speed of t 2 is equal to the number of tokens of p a (or p b ) allocated to this production line.
While the event handlers, v 1 , v 2 and v 3 determine the relationship between the input and output material of the production lines, the intensity handlers specify the speed of these lines according to the number of active tokens assigned to each line. In particular, the intensity edges {p a , s 1 } and {p a , s 2 } model the fact that the active tokens of p a can be used either by s 1 or s 2 . In a similar way, the fact that the active tokens of p b can be used either by s 2 or s 3 is modeled by the intensity edges {p b , s 2 } and {p b , s 3 }. This way, s 2 is the intensity handler responsible for the synchronization of resources for t 2 . The actions of all transitions are forced to be executed in order to model the active tokens to make the production lines work. Note that the graphical representation of the system by an FN is reasonably clear and compact. If the system were modeled by a classical Petri net, each transition would have to be split into several transitions that would each model the acquisition and release of the resources and the speed of the production lines. This would lead to a more complex graphical notation and, potentially, to more involved analysis methods.
Let the initial marking of the net be m there are two copies of resource type p a and one copy of resource type p b . Assume that the goal is to compute how the resources must be allocated over time so that the objective functionm[p 2 ]+0.5m[p 4 ]+0.25m[p 5 ], wherem[p i ] denotes the average marking of p i , is maximized. In words, this objective function implies that the goal is to maximize the production of all items giving priority to the products of type p 2 , then p 4 and finally p 5 .
This resource allocation problem can be solved by a single programming problem (see first method in Section 5) that makes use of the reachability constraints in Júlvez et al. (2018) and the mentioned objective function. More precisely, in order to obtain time trajectories, 90 intervals, each of 0.05 time units, will be considered.
The time trajectories of the marking and the allocated resources are shown in Fig. 9a and b respectively. Four time periods with different resource allocations (or operation modes) can be distinguished in these figures. The first period, from time 0 to 1.25, allocates the two tokens of p a to s 1 . This gives a high yield in the production of the items in p 2 which has the highest priority. Given that the two tokens of p a are used by s 1 during this first period, the token in p b cannot be used by s 2 , and hence it is used by s 3 to produce the items in p 5 , which has the lowest priority. During the second time period, from time 1.25 to 1.75, one token of p a is used by s 1 , and the other token of p a is synchronized by s 2 with the token of p b to operate t 2 and produce the items in p 4 , which has medium priority. As a result, the speed of t 1 (t 2 )(t 3 ) is 1(1)(0) during the second time period. At time 1.75, the marking of p 1 becomes 0, and hence, the active token of p a allocated to s 1 is released and becomes idle. Thus, during the third period, from time 1.75 to 3.25, only t 2 is working. At time 3.25, the marking of p 3 becomes 0, and hence the two tokens of p a become idle. During the fourth period, from time 3.25 onward, only the token of p b is active, and is employed by s 3 to operate t 3 .

Control actions
Control actions can be modeled in FNs by means of default intensities. This section demonstrates the ability of FNs to model and solve a control problem in which the control action is dynamically constrained. Figure 10a depicts a net with three places and four transitions. All the tokens are forced to be active, and all the actions are forced to be executed. The initial marking is m 0 [p 1 ]=m 0 [p 2 ]=0 and m 0 [p 3 ]=9. The default intensities of t 1 , t 2 and t 3 are 0. The default intensity of t 4 , λ 0 [t 4 ], models the only control action that can be applied to the system, and is constrained to the interval [0, 1.5]. Given that the equations associated with s 4 are s 4 :y=x; z=2x, each intensity unit in t 4 increases the intensity in t 1 , λ[t 1 ], by one unit and decreases the intensity in t 2 , λ[t 2 ], by two units. This way, the same control action is used for the intensities of t 1 and t 2 . Thus, the intensities in transitions satisfy  The evolution of the system can be described by the following differential equations: ] = 9 holds, then, the control objective is to drive the system to a marking that is as close as possible to the target marking (1, 4, 4). Figure 10b and c show the trajectories obtained by MPC with a sample time of 0.1 time units and a prediction horizon of one step. Initially, the value of λ 0 [t 4 ] is low as it is constrained by m[p 2 ], which initially is 0. Then, λ 0 [t 4 ] increases so that m[p 2 ] increases and m[p 1 ] decreases. At time 1.0, λ 0 [t 4 ] hits the constraint 1.5 where it is kept constant for 0.4 time units. Then, λ 0 [t 4 ] decreases in order to approach further the target marking. At steady state, the average intensities of t 1 , t 2 and t 3 are the same and equal to 2.73, and the value of the control action is λ 0 [t 4 ] = 0.81. The steady state marking reached is (1. 93, 4.35, 2.73). It is important to note that the target marking (1, 4, 4) cannot be an achievable steady state marking with the proposed single control action.
All the trajectories in this paper have been obtained by the tool fnyzer (https:// bitbucket.org/Julvez/fnyzer.git). This tool makes use of the modeling language Pyomo (Hart et al. 2017;Hart et al. 2011) and solvers, such as Gurobi (Gurobi Optimization 2015) and CPLEX (IBM ILOG CPLEX Optimizer 2010), to solve the programming problems associated with the FNs. The CPU time (Intel i7, 2.00 GHz, 8 GiB, Ubuntu 14.04 LTS) to solve one step of the MPC approach for the FNs in Figs. 6 and 10 was 1.81s and 5.93s respectively. The CPU time to solve the only programming problem associated with Fig. 8 was 1.27s.

Conclusions
FNs consist of two nets, an event net and an intensity net, that make an explicit distinction between the parts of the system involved in updating the marking in places, i.e. the event net, and the parts of the system involved in the determination of the speeds of transitions, i.e. the intensity net. Both the event and the intensity net are tripartite graphs in which places and transitions are connected by event and intensity handlers, respectively. This way, handlers act as an intermediate layer between places and transitions, which results in a significant modeling power. For instance, a transition in an event net can consume tokens from different sets of places, and a place in an intensity net can regulate the speed of different transitions. The tripartite net structure of event and intensity nets has demonstrated to be useful to model partial observability and resource allocation. Different types of system uncertainties can be accommodated by FNs through sets of linear inequalities associated with places, transitions, event handlers and intensity handlers. Namely, these inequalities allow the modeling of uncertainty in: a) the initial marking (8); b) the default intensities (15); c) the marking change produced by the execution of actions (4); and d) the intensity change produced by the active tokens (13). FNs account for the potential system trajectories arising as a result of uncertainties by means of a set of constraints that represent necessary reachability conditions. The combination of these constraints with an objective function can be used to obtain a system trajectory that optimizes a given criterion. This approach was successfully used to compute trajectory bounds, for instance, in the presented linear system with uncertain parameters, or to obtain a control law in a system whose control action is modeled by a transition with uncertain default intensity.