Array-Aware Matching: Taming the Complexity of Large-Scale Simulation Models

1 INTRODUCTION

In modern engineering, dynamic modelling and simulation are ubiquitous [16, 64]. Besides providing “virtual prototypes” [52] to streamline plant and control design activities [76, 77], “Digital Twins” – significantly based on simulation [9] – are nowadays the backbone of advanced controls [6, 55], predictive, condition-based and autonomous maintenance [10, 24, 43], anomaly detection, forecast and mitigation [37, 49, 72], continuous integration [54], lifelong asset management [47] and many other applications, see for example the survey in [25]. As a result, unprecedented challenges need facing for rapidly creating, modifying and running simulation models of steadily growing size and complexity [30].

Focusing on simulation models made of Differential and Algebraic Equations (DAE), the scenario just sketched has boosted the adoption of declarative, Equation-Based (EB) modelling languages as opposite to procedural, imperative programming ones [79]. The key feature of EB languages is the ability of separating the activities of writing a model and of producing its solution algorithms. This ability stems from the fact that the fundamental statement in EB languages is the equation. Contrary to the assignment, where an l-value receives the result of computing the expression on the right hand side, an equation just prescribes that the expressions on the left and the right hand side must be made equal – within convenient tolerances – at every point in time when the solution of a model is computed during its simulation.

Said otherwise, while assignments directly compose the algorithm to compute the model solution, equations just impose constraints to that solution, therefore saying nothing about the solver (in general, numeric) that will be used to compute it. To synthetically express this separation between describing the model and computing its solution, EB models are called declarative.

In synthesis, then, EB languages relieve the analyst from the task of turning equations into imperative code to perform their numerical integration, significantly helping to tame the complexity and rapid evolution of modern simulation problems [23]. It is the task of a translator to automatically turn a declarative model into an equivalent code in some imperative programming language [34], which is then fed to a compiler.

This translator-compiler workflow was devised at the outset of EB languages, with the aim of decoupling the generation of imperative code (translation) from its optimised compilation into machine code—a task for which e.g., C compilers are very well suited. However, today’s modelling and simulation problems exhibit new characteristics, that require to re-discuss the above translation workflow. A prominent such characteristic, on which we focus in this paper, is the presence of large (and possibly multi-dimensional) array variables and equations. This feature is distinctive of “large-scale” models. Think for example of a 3D thermal model for a solid body with fine-grained spatial discretisation: the model will contain energy dynamic balance equations for each of the many subvolumes into which the solid will be partitioned, and these equations will account for thermal exchanges with the adjacent volumes. It is quite natural to write such a model compactly in EB form by defining suitable array variables and equations, in the latter case by means of looping constructs (examples follow starting from Section 4).

In such cases, as we will show, a trade-off is easily observed. Writing the model directly as imperative code is far more complex, error prone and hard to maintain than adopting an EB declarative framework and obtaining the imperative code by automatic translation. But on the other hand, the code obtained by automatic EB-to-imperative translation is significantly less efficient than the one manually written as imperative.

We argue that the origin of this inefficiency mainly resides in the way EB translators manage array variables and equations. Current production-grade EB translators just treat each component of an array variable or equation as an individual scalar one, which results in a loss of structural information that imperative language compilers cannot efficiently recover [7]. As such, we argue that to improve both the translation time and the efficiency of the produced imperative code, it is necessary to make the translation process “array-aware”.

In this paper we offer a contribution to this end, aiming both for an efficient translation and an efficient imperative code. In detail,

To the best of our knowledge, we are the first to introduce an idea of matching optimality tied to the efficiency of the obtained imperative code, as well as to propose a heuristics that aims for that optimality besides for a fast matching process.

Organisation of the paper. The rest of the paper is organised as follows. In Section 2 we briefly introduce some definitions and the architecture of state-of-the-art approaches. Then, we discuss the history of EB model translators and other array-aware approaches in Section 3. We delineate the scope and purpose of our contribution in Section 4. In Section 5, we define array-aware matching and its optimality metric, while we prove its NP-completeness in Section 6. Finally, in Section 7, we show an approximate algorithm for array-aware matching, and in Section 8 we draw some conclusions and highlight future research directions.

2 BACKGROUND ON MODEL TRANSLATION

In this section we outline the foundations of automatic EB-to-imperative translation and compilation. To avoid confusion between scalar and array problems, we first provide a few definitions.

It follows that, in EB modelling languages, array variables are akin to the concept of ordinary multidimensional arrays, commonly found in many programming languages. Also, in EB modelling languages, an array equation is obtained by encasing a scalar one in one or more nested looping constructs, that define the indices of the contained scalar variables within the arrays of which they are part, as well as the ranges for the said indices.

It is important to highlight that the mentioned looping constructs – differently from those of imperative programming languages – allows the user to predicate on array equations, such as (1) \(\begin{equation} x[i]=y[i]\; \forall i \in [1,3]; \end{equation}\) to be intended as an abbreviation of (2) \(\begin{equation} x[1]=y[1];\;x[2]=y[2];\;x[3]=y[3]; \end{equation}\)

This noted, the translation-compilation process – as per the current state of the art in both research and production tools, see e.g., [23, 33] and [1, 2, 4, 5], respectively – can be divided into the following steps.

After the translation and compilation process is complete, the obtained executable code is run to produce the simulation output, in the form of a table with the value of the model variables as a function of simulation time. This process is often part of a graphical modelling environment for rapid prototyping. Once the modeller decides to perform a simulation, all the translation, compilation and execution steps are on the critical path toward getting the simulation results. Consequently, the shortening of compilation and simulation times is especially important for this kind of use-cases.

3 RELATED WORK

In this section we spend some words to relate our proposal to the research scenario on EB modelling in general, and to neighbouring research on EB-to-imperative translation in particular.

In the landscape of simulation languages, EB ones appeared and gained visibility in the 80s/90s of the last century; notable examples are Omola [11], Dymola [28] and gPROMS [15]. A taxonomy would stray from the scope of this paper; the reader interested in the historical panorama can refer e.g., to [22]. Worth noticing, however, is the common ancestor for the boost of the declarative approach stemming from studies such as [13] and [29], where the idea that simulation languages had to abandon the imperative setting was set forth and preliminarily exploited. Research on the matter thus focused at first on the model manipulation [41, 50] required by going declarative, and after a long systematisation process, this resulted in the birth of the Modelica language [51] to which we refer herein (though all the ideas we propose are general to the EB context).

The engineering use of EB languages and tools sustainably spread out in various domains, ranging from the chemical [12] and process industry [35] to power generation [21] and transmission [71], mechatronics [75] and robotics [39], automotive [42], and vehicles at large [26], aerospace [53], buildings [68] and more, including control design [20] and diagnostics [17]; the papers in the necessarily limited list above also contain interesting bibliographies for the reader willing to investigate further.

Together with testifying the success of the EB approach, however, the expansion just mentioned also shed light on some relevant limitations of the existing EB tools [62] – not of EB languages by themselves, it is worth stressing – especially when dealing with large-size models [19, 46]. This was the motivation for a first wave of tool optimisation, having as a major point the introduction of sparse solvers, a well-treated and long-lasting matter in domain-specific tools – see e.g., [32, 45, 69, 70] – but a source of challenges in the inherently multi-domain EB one [57, 78]. Examples of this research – with specific reference to Modelica given our scope – are [18, 63].

The possibility of solving large models fast enough to widen the EB applicability perimeter evidenced however a second type of tool limitation, concerning the translation rather than the solution of such models [66, 67]. The matter became critical in recent years, together with the emergence of problems that require model-based prototyping [48] and can scale up to the order of \(10^5\) equations. When such models become part of the inherently iterative engineering process, the time spent in translating and compiling them can be comparable to that spent in running simulations, if not even dominant [14]. For the sake of clarity it is worth noticing that the million equation barrier was already approached in the past [74] and in some domains nowadays well trespassed [61] by simulation tools, but these tools are not of the EB type, and most notably, do not separate model description and solution—which is a primary goal of the EB approach.

As a result, EB tools are nowadays undergoing a second wave of optimisation, directed to efficient translation. Open-source parsers [59] for EB languages such as Modelica are available providing some degree of array preservation, thereby enabling the research community to experiment with making the translation pipeline array-aware. In this relatively new effort, a primary objective is to achieve an \(\mathrm{O}(1)\) scaling of the translation time with the size of the model arrays, that as already noted are the main cause for the inefficiency of scalar-only model manipulation. In this context, the nearest neighbouring work to our research is the paper by Zimmermann et al. [80], who introduce the concept of “set-based graph” as a means to re-state the matching problem (originally scalar) in such a way to achieve an \(\mathrm{O}(1)\) translation, together with proposing algorithms for other manipulation steps related to matching, such as the management of strongly connected components and scheduling.

The main difference of our research with respect to [80] is a twofold instead of a single goal. More precisely, we do not aim just for an efficient translation, but also for an efficient simulation code. If the efficiency of the simulation code is taken into account, the number of looping constructs in the EB model that are preserved in the imperative one comes to matter a lot. The set-based approach of [80] is not designed to take this aspect into account. Aiming at loop preservation straightforwardly entails the introduction of an idea of optimality. This moves the focus from array-aware matching to optimal array-aware matching, and owing to its NP completeness, to the need for heuristics.

Other works have addressed the array-aware matching problem, such as [65] which correctly noticed that preserving looping constructs can positively impact both translation and simulation time. They also present a prototype translator that is limited to handling systems without algebraic loops and resorts to flattening the model completely in cases where their array-aware matching algorithm does not produce a solution. The paper [56] addresses array-aware index reduction, and the corresponding Modia implementation [3] also includes an array-aware matching algorithm that however is very simple and cannot split variable and equation nodes, thus requiring the modeller to pre-process the input model so as to make sure that every array variable can be matched to exactly one array equation. Compared to the two previously quoted works, our paper introduces the concept of optimal array-aware matching and proves that achieving such optimality is an NP complete problem, as well as presenting a more complete array-aware matching algorithm.

4 RESEARCH MOTIVATION

In this section, we discuss the inefficiencies that arise from a non array-aware model translation. To ground the discussion on an example, consider the model of a thermally insulated metal wire with prescribed temperatures at its ends. The evolution of this system is ruled by the one-dimensional Fourier equation. Carrying out a uniform spatial discretisation with the finite-volume approach results in the following system of differential equations: (3) \(\begin{equation} c \dot{T}_i = {\left\lbrace \begin{array}{ll} g (2 T_{{\it left}} - 3 T_i + T_{i+1}) & i = 1 \\ g (T_{i-1} - 2 T_i + T_{i+1}) & \forall \; i \in [2, N-1] \\ g (T_{i-1} - 3 T_i + 2 T_{{\it right}}) & i = N \end{array}\right.} \end{equation}\) where \(N\gt 3\) is the number of finite volumes, c is the thermal capacity of a volume, g the thermal conductance between the centres of two adjacent volumes, \(T_i\) is the temperature of volume i, with \(T_1\) being the leftmost and \(T_N\) being the rightmost volume, and finally \(T_{{\it left}}\) and \(T_{{\it right}}\) are the prescribed side temperatures.

When expressed in Modelica, the wire model reads as follows.

Listing 1.

View Figure

The remarkable similarity between the Modelica model and original system of equations is apparent; the der operator is used for expressing the derivative with time.

It is evident that the availability of array equations make EB models assume a compact and easily readable form. However, present state-of-the-art translators provide looping constructs only as a convenience for the modeller, and do not take advantage of this structural information to improve the translation efficiency, nor that of the generated simulation code. As a consequence, some inefficiencies arise that would be completely unexpected in the world of imperative programming languages.

A first such inefficiency is that the amount of both time and memory needed for the translation scale superlinearly with the size of array equations, rather than exhibiting the expected \(\mathrm{O}(1)\) scaling typical of imperative languages. In fact, the compilation of imperative programming languages treats looping constructs explicitly. The instructions inside a loop are represented in the compiler data structures just once, irrespective of the loop iteration count, and all subsequent compilation steps are performed on this compact data structure.

On the contrary, state-of-the-art EB translation algorithms are designed to only work in terms of scalar equations (see Flattening in Section 2), which are stored in the translator data structures individually. These scalar equations are then passed on to subsequent translation steps, some of which scale superlinearly. It is not uncommon for models exceeding the \(10^5\) equations mark to require translation times in the order of hours, and working memory in the order of hundreds of gigabytes [14].

A second and consequent inefficiency is that the size of the produced imperative code scales linearly with the size of array variables and equations. Since the original EB model is scalarised, the simulation consists of procedures containing repetitive code instead of looping constructs, often amounting to several gigabytes for large-scale models. Although imperative compilers are very efficient and can achieve \(\mathrm{O}(1)\) scaling in compilation with respect to array sizes, this is only possible if they are given a source code with loops, not long lists of repetitive statements.

The last inefficiency regards the performance of the machine code produced by the compiler when fed with the automatically generated imperative code. Modern computer architectures are built upon assumptions such as the locality principle for both data and code, which are the theoretical foundations for caches [27] and other microarchitectural optimisations. However, the code produced by current-generation model translators is not able to exploit these optimisations. First, the execution of large blocks of straight-line code requires frequent instruction cache invalidation. Furthermore, the loss of structural information about arrays leads to the generation of code that exhibits irregular data access patterns. Therefore, repetitive machine code is not only large, but runs significantly slower than equivalent hand-written code, in some cases 100 times or more [8].

Both translators and compilers can in principle infer looping constructs and improve access locality, but their ability to perform this kind of optimisation is limited [44]. Additionally, such inferences require to repetitively scan the list of flattened statements, making \(\mathrm{O}(1)\) scaling apparently impossible. We thus argue that a far better approach would be to preserve array-awareness throughout the translation process, rather than try to regain it a posteriori.

Summing up, all the presented inefficiencies share their root cause in the fact that existing model translation algorithms are not array-aware. In this paper, we begin an effort to fill this gap. Since the flattening step is trivial to extend for array-awareness – it suffices to not unroll looping constructs – we focus on the second and first truly key step, i.e., the matching problem.

5 GRAPH REPRESENTATION AND PROBLEM STATEMENT

In this section, we introduce the notation to describe array-aware algorithms, and present the formal statement of the array-aware matching, the algorithm we focus on in this paper.

We denote with \(G=(N,D)\) a generic graph having set of nodes N and set of arcs D. Any set of scalar equations E in the scalar variables V can be represented with a bipartite graph in which

In this context, the set of arcs D represents dependencies between variable and equation nodes.

To include array variables and equations – including multidimensional ones – we extend the above notation as follows. We indicate with \(\underline{v}_i\) the generic multidimensional array variable, whose scalar components in the physical problems we target are invariantly real numbers. We denote by \(\zeta _i\) the dimensionality of \(\underline{v}_i\). We denote with \(K=\lbrace k_1, \ldots , k_\delta , \ldots , k_{\zeta _i}\rbrace\) the sequence of positive integers needed to reference a scalar component within the array variable \(\underline{v}_i\). Each generic index \(k_\delta\) ranges from 1 to the size of the corresponding dimension \(\delta\), with \(1 \le \delta \le \zeta _i\). For compactness, we synthetically write \(v_{i,K}\) to indicate the scalar component in \(\underline{v}_i\) of indices \(\lbrace k_1 \dots k_{\zeta _i}\rbrace\).

Consistently with the above notation, we define the following.

A model containing both scalar and array variables and equations can thus be represented with any of the two equivalent bipartite graphs defined as follows.

The above definitions imply that the array graph is homomorphic to the flattened one. However, for the two graphs to be equivalent, each arc in the array one needs to carry information about which components of the connected array variable and equation it refers to — a matter that does not pertain to the scalar case. Therefore, the arcs of the array graph must be endowed with the information needed to reconstruct the arcs of the flattened graph. To formalise this, we introduce the concept of local multidimensional incidence matrix. Given an array equation \(\underline{e}_i\) and an array variable \(\underline{v}_j\) of dimensionality \(\zeta _i\) and \(\zeta _j\) respectively, let \(K_i=\lbrace k_{i,1},\ldots k_{i,{\zeta _i}}\rbrace\) be the sequence of indices for \(\underline{e}_i\), and \(L_j=\lbrace l_{j,1},\ldots l_{j,{\zeta _j}}\rbrace\) be the sequence of indices for \(\underline{v}_j\). The multidimensional local incidence matrix \(u_{i,j}\) has dimensionality \(\zeta _{ij}=\zeta _i+\zeta _j\), and its sequence of indices \(Q_{ij}\) is the concatenation of \(K_i\) and \(L_j\). Its generic element is 1 iff the scalar variable \(v_{j,K}\) appears in the scalar equation \(e_{i,L}\), else it is 0. We name U the set of local multidimensional incidence matrices. An example of a Modelica model and the corresponding array graph can be found in Figure 1.

Fig. 1.

View Figure

The notation we just provided is capable of representing arbitrary multidimensional incidence matrices. However, it should be noted that those produced by EB models coming from equations of physics are significantly structured, and present patterns that arise out of the looping constructs and expressions used for accessing array variables within array equations. It follows that, although from a theoretical standpoint it is convenient to reason in terms of multidimensional incidence matrices, an industry-grade implementation would rely on an efficient pattern-based data structure to achieve \(\mathrm{O}(1)\) scaling. This important matter is discussed in Section 7.4.

5.1 The Array-aware Matching Problem

The array-aware matching is an operation that takes as input an array graph \(\underline{G} = (\underline{V} \cup \underline{E}, \underline{D})\), where we recall that \(\underline{D}\) is the set of dependencies between array variables \(\underline{V}\) and array equations \(\underline{E}\). Additionally, every dependence \(\underline{D}\) has an associated local multidimensional incidence matrix \(u_{i,j}\). The array-aware matching produces as output an array graph \(\underline{G}^{\prime } = (\underline{V} \cup \underline{E}, \underline{D}^{\prime })\) with the following properties.

(1)	\(\underline{D}^{\prime } \subseteq \underline{D}\)
(2)	Each arc \((\underline{v}_{i}, \underline{e}_{j}) \in \underline{D}^{\prime }\) has an associated local multidimensional incidence matrix \(m_{i,j}\), where \(m_{i,j} (\underline{v}_{i,K}, \underline{e}_{j,L}) \in \lbrace 0,1\rbrace\). We call M the set of \(m_{i,j}\).
(3)	\(\forall m_{i,j} \in M\), \(m_{i,j}\) has the same size and dimensionality of \(u_{i,j}\in U\)
(4)	\(\nexists (\underline{v}_{i}, \underline{e}_{j}) \in \underline{D}^{\prime } : m_{i,j}\) is a zero matrix
(5)	\(u_{i,j} (\underline{v}_{i,K}, \underline{e}_{j,L}) = 0 \Rightarrow m_{i,j} (\underline{v}_{i,K}, \underline{e}_{j,L}) = 0\)
(6)	\(\forall \underline{v}_{i,K} \in \underline{V} : \exists ! m_{i,j} (\underline{v}_{i,K}, \underline{e}_{j,L}) = 1\)

Property (1) states that no new arcs are added to the graph, (2–3) denote as “matching matrices” m the local incidence matrices of the output graph, (4) states that – consistently – arcs fully unmatched are removed, (5) ensures that a scalar equation is matched only with a scalar variable it contains, and (6) that each scalar equation is matched to one and only one scalar variable.

Definition 5.1

(Optimal Array-aware Matching).

Let us define the expansion function \(f_e(m)\), which maps a local incidence matrix m to the number of looping constructs required to implement the matching it expresses. \(f_e(m)\) has a value of zero if and only if m is a zero matrix, otherwise \(f_e(m) \ge 1\). We say that a matching \((\underline{G}^{\prime }, M)\) is optimal if it minimizes the following: (4) \(\begin{equation} \Omega (\underline{G}, M) = \sum _{m \in M}{f_e(m)} \end{equation}\)

For the purpose of this work, it is not necessary to fully specify the expansion function, because different model translators may be able to “efficiently” handle only a subset of all possible matching matrices. In this discussion, the metric of “efficiency” is therefore the ability to represent the matching described by \(m_{i,j}\) as a single looping construct consisting of the scalar equations of the node adjacent to the corresponding arc \(d_{i,j}\). Without loss of generality we can consider the ideal case in which every matching matrix maintains its correspondence with exactly one looping construct. In this case, the \(f_e\) function takes this form: (5) \(\begin{equation} f_e(m) = \left\lbrace \begin{array}{ll} 0 & \mbox{iff $m$ is a zero matrix}\\ 1 & \mbox{otherwise}\\ \end{array} \right. \end{equation}\)

In this case, the optimality metric \(\Omega\) is equal to the number of arcs where at least one equation is matched. Intuitively, the more arcs that have a matching, the higher the number of times the same array equation has been matched with multiple variables, and for each variable matched with the same array equation a separate looping construct must be introduced. Therefore, minimizing \(\Omega\) means minimizing the number of looping constructs.

5.1.1 Just Enough, but Not Too Much.

It is important to stress that the optimality metric \(\Omega\) is but a proxy for identifying the matching that best improves both translation and simulation time. With the definition we gave of array equation in this paper (Definition 2.3 in Section 2) where individual scalar equations in an array equation differ only by the array indices, reducing the number of looping constructs does indeed result in more efficient code. Other works however, such as [56] and the corresponding Modia implementation [3] allow if statements in array equations thereby allowing to merge in a single array equation also scalar equations that are structurally different. For example, the wire model written with such an extended definition of array equations would read as

View Figure

Although it would appear that such an array equation could improve the optimality metric by requiring a single looping construct that handles all the wire finite volumes instead of requiring one looping construct plus two scalar equations for the first and last volume, this would not result in efficient code. Indeed, if such an array equation were brought as-is till code generation, the imperative code would need if statements in the loop to handle the differences in the equation structure, and the cost at run-time of the if statements would need to be paid multiplied by every loop iteration and additionally multiplied by every simulated time step. This is a strong point in favour of our definition of array equations, that explicitly disallows grouping structurally different scalar equations in a single array equation. Moreover, to achieve the best efficiency, a translator would need to perform symbolic manipulations to transform inefficient code such as the one above by moving the structurally dissimilar equations outside of the loop before the matching step, should the modeller decide to write the model in that form.

6 COMPLEXITY OF OPTIMAL ARRAY-AWARE MATCHING

In this section we show that for what concerns the matching problem, preserving the array structure of both equations and variables results in NP-completeness. This theoretical result motivates the need for heuristic algorithms, that will be presented in Section 7.

Proof.

We prove Theorem 6.1 by reducing the max-2-sat problem to optimal array-aware matching. Max-2-sat was proven NP-complete by Garey, Johnson and Stockmeyer [36], and reads as follows: given a Boolean formula in conjunctive normal form, where each clause contains at most two literals, find a literal assignment such that the maximum number of clauses is satisfied.

First of all, we establish a procedure that allows to represent a max-2-sat problem in terms of array-aware matching. To this end we introduce three formal modifications to the way max-2-sat is expressed in order to simplify the reduction process.

(1)	Every clause in the form \((a)\) is rewritten as \((a \vee a)\). It is evident that this rewriting does not change the value of that clause for any literal assignment.
(2)	Instead of considering the input to max-2-sat as a conjunction of clauses, we consider it as an ordered list. This is legitimate, as maximising the number of satisfied clauses does not require knowing if the entire formula is satisfied.
(3)	We assume that, when traversing the list of clauses, the first encountered occurrence of each literal is not negated. Should this be false, one would simply have to replace that literal with another one defined as its complement.

The above said, to encode a max-2-sat instance into an array-aware matching one, we start by recalling that \(\begin{equation*} (a \vee b) = (a \wedge \lnot {b}) \vee (\lnot {a} \wedge b) \vee (a \wedge b) \end{equation*}\) and we replace each OR clause in the list with the above defined equivalent triple of AND ones. Observe that only one of the three AND clauses can be true, hence maximising the number of satisfied AND clauses in the reformulated list is equivalent to the original problem.

Let us start with two empty sets of nodes \(A = \emptyset\), \(B = \emptyset\). Then, we build an intermediate flattened bipartite graph \(G = (A \cup B, L)\) by traversing the list of AND clauses. By construction, we will retain the invariant \(A \cap B = \emptyset\). For each clause c we operate as follows.

First, we consider the first literal \(\ell\) in the clause. If this is the first occurrence of \(\ell\), we add to G the subgraph shown in Figure 2(a). In doing so, \(\alpha\) nodes should be considered belonging to the A set while \(\beta\) nodes should be considered part of the B set. We name \(e_{\ell ,1}\) the first edge associated to literal \(\ell\). We name \(\alpha _{\ell ,1}\) and \(\beta _{\ell ,2}\), respectively the start and the end node for that literal.

Fig. 2.

View Figure

If \(\ell\) was already encountered and appears here in non-negated form, we add to G the subgraph shown in Figure 2(b), and we connect \(\alpha _{\ell ,i}\) to the end node of the literal. The start node of the literal does not change, while its end node becomes \(\beta _{\ell ,i}\).

If \(\ell\) was already encountered and appears here in negated form, we add to G the subgraph shown in Figure 2(c), and we connect \(\alpha _{\ell ,i}\) to the end node of the literal. The start node of the literal does not change, while its end node becomes \(\beta _{\ell ,i}\).

We repeat the subgraph insertion process for the second literal in the clause.

For each clause c, we define a set of four nodes – taken from the two subgraphs just added – as the union of nodes highlighted in red in Figure 2. We name this set clause nodes of c and we name it \(N_c\).

We repeat the above for all clauses. When the end of the list is reached, we connect the end node of each literal with its start one. Doing so, we create for each literal a simple cycle \(\alpha _{\ell ,1}\), \(\beta _{\ell ,1}, \ldots , \alpha _{\ell ,n_\ell }\), \(\alpha _{\ell ,1}\) of even cardinality \(n_\ell\), ordered as just indicated. We number the edges of each cycle as \(l_{\ell ,1}\) through \(l_{\ell ,n_\ell }\). This concludes the construction of graph G.

Graph G is a flattened graph and as such does not have the form required by array-aware matching, thus we need to construct a different graph \(\underline{G}=(\underline{E} \cup \underline{V}, \underline{D})\) homomorphic to G, and expressed in terms of array equations and variables. Therefore, \(\underline{G}\) will satisfy the following properties:

•	For each node \(\alpha _{\ell ,i} \in A\) not part of a clause node \(N_c\), there exists a variable \(\underline{v}_{n} \in \underline{V}\) with dimensionality 1 and size 1.
•	For each node \(\beta _{\ell ,j} \in B\) not part of a clause node \(N_c\), there exists an equation \(\underline{e}_{m} \in \underline{E}\) with dimensionality 1 and size 1.
•	For each clause node \(N_c\) there exists a variable \(\underline{v}_{n} \in \underline{V}\) and an equation \(\underline{e}_{m} \in \underline{E}\), both with dimensionality 1 and size 2.
•	The local incidence matrix \(u_{n,m}\) is a square identity matrix in order to represent the original relationship found in the scalar graph.
•	No other variables nor equations exist in \(\underline{G}\).
•	Dependencies arcs \(\underline{D}\) are constructed so as to make \(\underline{G}\) homomorphic to G.

In other words, all the \(\alpha\) nodes are considered variables, all the \(\beta\) nodes are considered equations, and each labelled set of four nodes \(N_c\) forms an array equation of size 2, and a corresponding array variable of size 2, while unlabelled nodes translate to scalars.

To carry on, we now need the following lemma.

Lemma 6.2

(Complexity of the Construction of G)

The graphs G and \(\underline{G}\) can be constructed in polynomial time.□

Proof.

The construction process which defines G is linear with the number of clauses and the creation of \(\underline{G}\) can be done by simply enumerating the nodes and edges of G.□

Back to the main proof, once we obtained the matched graph by executing array-aware matching on \(\underline{G}\), we assign to each literal \(\ell\) the boolean value true if arc \(l_{\ell ,1}\) has been selected for the matching, false otherwise. This last step can also be performed in polynomial time.

Given the definition of optimal array-aware matching of Definition 5.1, an algorithm capable of solving the problem will maximise the number of \((l_{\ell _1,i},l_{\ell _2,j})\) arc pairs between nodes belonging to the clause node sets \(N_c\). In fact, each pair contributes only a unitary weight to the optimality metric \(\Omega\), while two non-paired arcs will contribute a weight of two.

By construction of the bipartite graph G, each node, be it an equation or a variable, has exactly two adjacent edges. Thus, the matching choice is binary. Since we have built a simple cycle for each literal, there are only two matching solutions for each cycle, one selecting \(l_{\ell ,1}\) and all the odd numbered arcs, the other one selecting all the even numbered ones.

Additionally, due to how we constructed each cycle, the arc connecting the nodes belonging to a clause node set \(N_c\) for cases (a) and (b) of Figure 2 is odd, while in case (c) it is even. Thus, if the arc in the red box of Figure 2(a) is selected, then all the arcs in boxes of subgraph type (b) will also be selected, the arcs in boxes of subgraph type (c) will not be selected, and vice versa.

It follows that an arc pair inside a labelled node set can be selected if and only if the literal assignment – as read from the graph – satisfies the corresponding clause. As such, the objective functions of max-2-sat and array-aware matching, given the proposed graph construction and interpretation, are equivalent. This implies that array-aware matching is NP-complete.

Example. For the convenience of the reader, we complement the formal proof with an example of how the intermediate bipartite graph can be built from the following list of AND clauses: (6) \(\begin{equation} \lbrace a \wedge b,\; \lnot {a} \wedge c,\; c \wedge d\rbrace \end{equation}\)

The list of clauses is scanned left to right. The first clause encountered is \((a \wedge b)\), and both a and b are literals being encountered for the first time. As a result, in the graph we add two structures of type (a) shown in Figure 2. The start node of a is \(\alpha _{a,1}\), the end node of a is \(\beta _{a,2}\). In the same way, the start node of b is \(\alpha _{b,1}\), the end node of b is \(\beta _{b,2}\).

Now, the second clause \((\lnot {a} \wedge c)\) is processed. The literal a was already encountered, and appears in negated form; as a result, a structure of type (c) from Figure 2 is inserted, and the \(\beta _{a,2}\) end node is connected with \(\alpha _{a,3}\). The end node of a is changed to the newly inserted \(\beta _{a,3}\). Instead, the literal c is newly encountered, and therefore we insert a structure of type (a) from Figure 2. The start node of c is \(\alpha _{c,1}\), and the end node of c is \(\beta _{c,2}\).

The process has now arrived at the third and last clause \((c \wedge d)\). The literal c makes a reappearance in positive form, and therefore a structure of type (b) is inserted. The end node \(\beta _{c,2}\) is connected to the new node \(\alpha _{c,3}\), and \(\beta _{c,3}\) becomes the new end node of c. Finally, the literal d is encountered for the first time, and a structure of type (a) is inserted. At this point the list of clauses is exhausted.

The last step for constructing the scalar graph consists in connecting the end node of each literal with its corresponding start node. For literal a, we connect \(\alpha _{a,1}\) with \(\beta _{a,3}\). For literal b, we connect \(\alpha _{b,1}\) with \(\beta _{b,2}\). For literal c, we connect \(\alpha _{c,1}\) with \(\beta _{c,3}\). Finally, for literal d, we connect \(\alpha _{d,1}\) with \(\beta _{d,2}\).

Figure 3(a) shows the bipartite graph generated through the steps we have just outlined. Additionally, the figure outlines in red each clause node, and highlights the first arc of each cycle, which is used to determine the value of each literal from the matching.

Fig. 3.

View Figure

Figure 3(b) highlights the optimal matching of the same graph, alongside with the value of each clause and each literal. Since arc \(l_{a,1}\) has not been used in the matching, a is assigned value 0. Instead, literals b, c and d are assigned a value of 1, since the arcs \(l_{b,1}\), \(l_{c,1}\), and \(l_{d,1}\) are all used in the matching. This matching maximises the number of array equations and variables matched, and therefore also minimizes \(\Omega\): two arrays are matched, out of the three arrays in the input graph. Each array corresponds to a clause, and the matched arrays represent clauses whose value is one. In fact, the literal assignment makes the second and third clause true, and the first clause false.

7 ARRAY-AWARE MATCHING ALGORITHM

The NP-completeness proof of optimal matching highlights how, in order to efficiently handle large-scale problems, there is a need to introduce suitable heuristics. In this section, we propose an algorithm of reduced complexity.¹

Our proposal is a partially heuristic two-step procedure. The first step identifies obligatory matching choices and removes them from the problem, thereby reducing its size. We call this first step, presented in Section 7.2, the simplification step. No heuristics are involved in it.

For the second step, that we call the matching step, we propose an extension of the Hopcroft-Karp [40] algorithm to array graphs, which allows to exploit local incidence matrices in such a way to preserve the existing looping constructs. The matching step, described in Section 7.3, terminates at the first solution found, whence its heuristic nature.

Before presenting the simplification and matching algorithms, it is necessary to define some operations on local incidence matrices that are used throughout the said algorithms. This is done in Section 7.1, but requires an important preliminary remark. When dealing with matching problems that contain multidimensional arrays, incidence matrices become multidimensional as well. To lighten the treatise, in this work we nonetheless stick to talking about “rows” and “columns”. This notation does not cause any generality loss, however. All the proposed algorithms can work in the case of multidimensional incidence matrices by simply interpreting “rows” and “columns” as sets of dimensions, referring respectively to equations and variables. We also talk about “row” and “column vectors”, that generalize to matrices having only one of the two sets of dimensions.

7.2 Simplification Algorithm

The simplification step is dedicated to performing all the obligatory matches. Its importance in the handling of real-world problems becomes evident when considering that, as explained in [23, Chapter 7, and more specifically Section 7.2], when matching differential equations, the state variables of the system are considered to be known and need not be matched. Thus, the variables to be matched are either non-state variables, or derivatives of state variables. In modelling the evolution of physical systems over time, it is very common to write equations in the form \(\begin{equation*} \dot{x_i} = \mathrm{f}(x,u) \end{equation*}\) where x and u are respectively the sets of state variables and inputs. These equations introduce equation nodes in the bipartite matching graph with a single arc connecting them to the corresponding derivative, and that thus can only be matched with one variable. For example, this applies to all the three arcs in Figure 1. The commonplace presence of obligatory matching options motivates the introduction of an efficient simplification step.

The proposed simplification step is shown in Algorithm 1. It takes as input the array graph \(\underline{G}\), and operates as follows. First, it constructs a set L from every node in the graph with only one arc, regardless of it being an equation or a variable.

Then, for every node \(\underline{n}_1\) in the set, it uses the solveLocalMatchingProblem procedure to attempt to match the only arc adjacent to that (array) node.

In the case of multiple matching options, the simplification step skips the considered node, because at this stage any matching choice would be arbitrary, and may affect the feasibility of the array-aware matching problem. Conversely, if only one option is found that that fully matches \(\underline{n}_1\), that option is included in the solution. Attention then shifts to \(\underline{n}_2\), the only node reached by \(\underline{n}_1\).

Since we have fully matched one node and matched some variables from \(\underline{n}_2\), the simplify procedure now checks whether also \(\underline{n}_2\) is fully matched. If it is, \(\underline{n}_2\) is removed from the set in case it was there, and since having fully matched \(\underline{n}_2\) may result in neighboring nodes with only one unmatched arc, all such nodes are added to the set. Even if \(\underline{n}_2\) is not fully matched, it too may end up having only one unmatched arc (if it previously had two), and in this case \(\underline{n}_2\) is added to the set. The simplify algorithm thus recursively eliminates all nodes with a constrained match, leaving only irreducible connected components, as well as nodes where arcs have multiple matching options. After the simplify step all matches found are removed from the graph, and only the remaining part of the graph is passed to the subsequent matching algorithm.

7.2.1 Optimal Matching in Polynomial Time.

It is evident from the formulation of the simplification step that its algorithm executes in polynomial time with respect to the number of nodes and arcs in the array graph and, with a suitable data structure to represent incidence matrices, in O(1) time with respect to the size of the arrays. We also notice that there is a class of array graphs that can be completely matched by application of the simplification algorithm alone. The model of a wire previously shown in Listing 1 is an example of that. By construction of the simplification algorithm, such graphs must have a single solution to the array-aware matching problem, and that solution is therefore optimal. It follows that such models can be optimally matched in polynomial time with respect to the number of nodes and arcs in the array graph, and in O(1) time with respect to the size of the arrays.

Finding other classes of graphs which can be optimally matched in polynomial time is an open research issue.

7.3 Matching Algorithm

Now we illustrate the complete array-aware matching algorithm. This algorithm must be able to match any valid array graph, and additionally it must attempt to approximate the optimal matching as much as possible. The algorithm we present is based on the well-known Hopcroft-Karp one for bipartite graph matching, with adaptations to support array graphs. Its main procedure is shown in Algorithm 2.

As we will see in the following, procedure augmentingPaths computes a list of non-intersecting augmenting paths in the graph to be applied later by the applyPath procedure. When no augmenting paths are found, the matching is considered complete.

Each augmenting path p, when applied, adds or removes matching from the graph. It consists of a list of tuples called steps \(p_i = (\underline{n}, \underline{d}_{i,j}, m^{\prime }_{i,j})\) where \(\underline{n} \in \underline{N}\) is the starting node of the step, \(\underline{d}_{i,j}\) is the arc being traversed, and \(m^{\prime }_{i,j}\) is the incidence matrix that specifies which elements of the local matching incidence matrix \({m}_{i,j}\) are modified by the step. In a similar way to the Hopcroft-Karp algorithm, steps starting from an equation node (\(\underline{n} \in \underline{E}\)) add non-zero matrix entries to the matching, and steps starting from a variable node (\(\underline{n} \in \underline{V}\)) remove entries from the matching. In other words, in a step starting from an equation node, the elements set to 1 in \(m^{\prime }_{i,j}\) are added to the matching matrix \({m}_{i,j}\). On the contrary, in a step starting from a variable node, the elements set to 1 in \(m^{\prime }_{i,j}\) are removed from \({m}_{i,j}\).

The computation of the augmenting paths is performed through a breadth-first-search in the residual graph. In contrast to the conventional scalar matching process performed by the Hopcroft-Karp algorithm, when computing the augmenting path we also need to keep track of the matching matrices that express the set of equivalent scalar variables and equations that are being matched or un-matched.

The breadth-first-search procedure must first be seeded with an initial frontier F by collecting the list of array equation nodes with at least one free scalar equation. This initial frontier is then passed to the bfs procedure, which computes both the initial list of augmenting path candidates L and the forest of search trees \(G_{BFS}\) traversed during the search. In \(G_{BFS}\), A is the list of nodes in the forest, and B is the list of arcs. Each node \(a\in A\) is a tuple \((n, v)\) where n is a node in the matching graph, and v is a binary vector with size and dimensionality equivalent to the one of n specifying which scalar equations are being traversed in the path. Similarly, each arc \(b \in B\) is a triple \((a, m, a^{\prime })\) where a is the parent node, \(a^{\prime }\) the child node, and m is the local incidence matrix that describes the equivalent traversed path in the homomorphic scalar graph. The augmenting path candidates \(l \in L\) are actually just leaves of the \(G_{BFS}\) forest.

In contrast to the basic Hopcroft-Karp matching algorithm, but similarly to the Ford-Fulkerson flow maximization algorithm [31], the initial “flow” in the first steps of the path can be different from the final “flow” at the last step. In this context the “flow” of a given step is simply the number of scalar equations affected. To perform this operation, each path is traversed backwards – from the end to the beginning – and it is modified such that the following invariant is respected: \(\begin{equation*} \begin{array}{ll} \mathrm{flattenRows}(m) = \mathrm{flattenRows}(m^{\prime }) & \forall \; (a, m, a^{\prime }), (a^{\prime }, m^{\prime }, a^{\prime \prime }) \in A : \text{$a^{\prime }$ equation node}\\ \mathrm{flattenColumns}(m) = \mathrm{flattenColumns}(m^{\prime }) & \forall \; (a, m, a^{\prime }), (a^{\prime }, m^{\prime }, a^{\prime \prime }) \in A : \text{$a^{\prime }$ variable node}\\ \end{array} \end{equation*}\) During this process, we build the augmenting path p from the nodes, arcs, and matching matrices traversed.

Additionally, the breadth-first-search does not immediately return a set of paths that respect the non-intersection condition already present in the Hopcroft-Karp algorithm. An augmenting path intersects another if the two paths, in at least one point, traverse the same node with intersecting matching matrices. In order to guarantee this property, each augmenting path is tested against the others. If two paths are intersecting, one of the two is discarded.

Discarding intersecting paths has the effect of eliminating multiple candidates which are in mutual exclusion between each other. The specific candidates being discarded at each step influence the solution \(\Omega\), and thus how close the solution is to the optimum. Additionally, they impact the number of steps required by the matching algorithm. In our current implementation, such selection depends on the ordering of L. Our heuristic (implemented in the heuristicSort procedure) sorts L based on the number of ones in the matching matrices of the path (paths with more ones are prioritized). Other heuristics could be devised to improve the solution \(\Omega\) and the number of steps required to complete the matching, and this could be an interesting direction for future work.

The breadth-first-search procedure is shown in Algorithm 3. It operates in the conventional way, with three additional constraints:

(1)	Each move in the search is associated with a vector of tangent elements \({f}\) to the destination node \(\underline{n}\) of the move. f is a vertical vector if \(\underline{n}\) is an equation, otherwise f is an horizontal vector.
(2)	Each move in the search must be associated with an incidence matrix s called path matrix that satisfies the same conditions imposed on matching matrices \(m_{i,j}\) (see Section 5.1, conditions (3–6)).
(3)	The incidence matrix s only affects the scalar variables or equations specified by f.

These constraints ensure that each path traversed during the search corresponds to a set of one or more equivalent paths in the scalar graph. As a result, multiple different moves can start at the same moment from the same node and through the same arc, but with different path matrices. The computation of the set of possible moves at each iteration is performed by enumerating the adjacent edges to the current node, and then by using solveLocalMatchingProblem to compute the set of distinct valid ways to traverse that edge. When a path reaches a variable node, and the path matrix contains at least one non-zero column not corresponding to any matched scalar variable, then it allows to augment the matching and the search is stopped. Since the search process is stopped for the entire frontier, all paths returned by the bfs procedure have the same length in terms of number of steps.

8 CONCLUSIONS AND FUTURE WORK

We discussed the problem of translating declarative EB models into imperative code, concentrating on the crucial step of equation/variable matching. Relating our work to the research scenario, we evidenced two issues to address. First, currently established approaches to EB-to-imperative translation handle array variables and equation looping constructs in an extremely inefficient manner. Second, array-aware proposals in the literature aim for a fast translation, but not for obtaining an efficient imperative code.

We showed that to pursue both the objectives above, one needs to translate an EB model in such a way to maximise the preservation of looping constructs in the imperative code, and we introduced a metric to measure the said preservation. This led us to define the concept of optimal array-aware matching, as the one that maximises that metric.

As our main methodological contribution, we proved that the problem of computing an optimal matching in the sense just stated is NP-complete. Motivated by this completeness, we proposed an algorithm to compute an array-aware matching in polynomial time, whose rationale is to stop at the first solution found, privileging however the preservation of looping constructs when choices need to be taken.

The ideas we presented are currently being put to work within the implementation of an experimental Modelica compiler [7, 8]. Indeed, based also on the effort that such a realisation entails, we do hope that in the next years the developments we described will be adopted in the EB modelling community at large. In fact, the advantages of an array-aware EB-to-imperative translation are essential for addressing large-scale models in an industrial context.

In the future, we also plan to extend our approach to the rest of the translation pipeline, addressing other problems such as equation scheduling and SCC resolution. Additionally, a more precise characterization of the set of array graphs that can be optimally matched in polynomial time is also of great interest, as well as better heuristics for the array-aware matching algorithm.

ACKNOWLEDGEMENTS

The authors would like to thank the anonymous reviewers for their helpful suggestions that allowed us to improve the paper presentation, in particular regarding the definition of array equation and its impact on the optimality metric, as well as the opportunity to perform optimal matching in polynomial time for certain classes of models.