Recursive modular modelling methodology for lumped-parameter dynamic systems

1. Introduction

Mathematical modelling has become a fundamental tool for scientific and technological projects in several areas of knowledge. Moreover, it has been acquiring increasing importance in multidisciplinary studies. Although the models of lumped-parameter¹ dynamic systems can always be expressed in terms of a finite number of differential–algebraic equations (DAEs), the modelling methodologies for obtaining them still depend on the nature of the system being modelled and on the simplifying hypotheses adopted. Thus, there is an extensive framework of methodologies for the mathematical modelling of lumped-parameter dynamic systems from various areas of knowledge. On the one hand, this can be understood as an advantage, once it is possible to use modelling procedures that are suitable for each type of system. On the other hand, this might represent an obstacle when someone tries to conciliate the use of different approaches in the modelling of a complex system.

An aspect that might be explored is the modularity of systems, which, in the context of this text, can be understood as the ability of conceiving a complex system from simpler systems by the enforcement of constraints. A modular methodology is capable of transforming the mathematical models according to the constraints that are imposed, regardless of how the original models have been obtained. Consider, for instance, a multibody mechanical system, which can be conceived as a set of constrained subsystems. Such a concept not only explores the topological symmetries that often exist in these systems (figure 1) but also allows the use of already known mathematical models of subsystems (which might have been obtained by any other formalism) as the starting point for the modelling procedure. Therefore, a modular methodology might not be understood simply as another alternative for the already extensive framework of existing ones, but rather as a methodology that aims to conciliate the use of others from which the models of the modules may have been obtained. Some examples of modular modelling methodologies originally developed for multibody systems can be found on works by Orsino [1], Orsino & Hess-Coelho [2] and Udwadia et al. [3–8]. The latter works are typically referred to in the literature as the Udwadia–Kalaba equation.

• Download figure
• Open in new tab
• Download PowerPoint

The objective of this paper is to propose a general modular modelling methodology that must be applicable to any lumped-parameter dynamic system (i.e. to any dynamic system whose model can be expressed by a finite set of DAEs), which is also recursive. This last property will allow exploration of the modularity of a system in many levels, leading to a hierarchical representation of it in which each level has its own mathematical model, but only the model at the top one is supposed to adequately describe the dynamics of the system. Such a new approach not only proves the equivalence between the previously cited modular methodologies but also generalizes them to a recursive form. Constraint enforcement algorithms are introduced for obtaining the dynamic model at a given level of a hierarchy, given the model at the inferior level and a mathematical description of the constraints to be enforced.

The generality of the methodology is ensured by the fact that the enforcement algorithms are based on no other hypotheses than the models being expressed by a finite number of DAEs and the virtual work of the generalized constraint forces associated to admissible variations in the state of the system being either zero or a known function given by a constitutive relation. There are no restrictions regarding:

• (i) the formalisms used in the derivation of the models at the inferior level of a hierarchy (the models may either follow from the application of physical principles using some analytical formalism or can be phenomenological);

• (ii) the order of the differential equations involved, which can be higher than two;²

• (iii) the criteria for selecting variables, which allows exploration of the advantages of formulations based on the use of redundant variables and on reparameterizations of time rates [1,2];

• (iv) the nature of the constraints and the techniques adopted for their description; and

• (v) how the user chooses to conceive the system modularly and how many levels there are in a given hierarchical representation.

This lack of restrictions allows the user to specialize the methodology for particular applications according to conveniences (e.g. when already available mathematical models are used) or personal preferences, leading to a great variety of possibilities for computational implementations based on the approach presented in this text. On the other hand, it makes the performance of these implementations (in terms of computational complexity and numerical stability) highly dependent on the techniques adopted for these specializations. Such a discussion is out of the scope of this paper, and the reader is invited to check some textbooks and a paper on the theme [11–14]. Basically, this paper aims to address the following topics:

• (i) prove the equivalence among projection methods for elimination of generalized constraint forces from equations of motion in the literature and the strategy used in the modular approach originally introduced by Orsino & Hess-Coelho³ [2];

• (ii) present a recursive form for the modular approach, allowing to apply the methodology directly to complex systems represented by multiple-level hierarchies; and

• (iii) prove that a recursive algorithm based on the Udwadia–Kalaba equation is identical to the algorithm of a Kalman filter for the estimation of the state of a static process subjected to a sequence of measurements modelled by noiseless observation equations.

Recursive algorithms are among the most popular ones for the modelling of complex dynamic systems, once they do not have to consider the whole complexity of the system at once. Particularly, in the area of multibody systems, recursive algorithms based on Newton–Euler equations stand out [11,14–17]. In these algorithms, the systems are represented by graphs, each node corresponding to a rigid body and each edge to a joint. In the case of inverse simulations, the motion is specified as a function of time and force and torque inputs are the desired unknown variables, which must be computed by the algorithms along with the constraint forces and torques. In the case of forward simulations, the initial value problem associated to the equations of motion of a given system must be solved, provided that all the force and torque inputs are specified as functions of time; the corresponding algorithms also require the values of the constraint forces and torques to be found at each time instant. Rodriguez [17] has explored the similarities between the formulations based on Newton–Euler equations of forward and inverse problems for open-loop kinematic chains and the algorithms of the Kalman filter and the Bryson–Frazier smoother. An analogy with the Kalman filter is also discussed in this text on the formulation of a recursive algorithm based on the Udwadia–Kalaba equation. Also, some of these recursive methods propose the use of extra algorithms for computing operators that eliminate the constraint forces from the equations of motion, before applying them to find the desired unknowns. The natural orthogonal complement (NOC) and its decoupled variant (DeNOC) [14–16] are some examples of operators obtained by such algorithms. This text discusses some methods, based on the same fundamental theorem (theorem 3.2), for computing such operators using orthogonal complements, generalized inverses or orthonormalization procedures.

This paper is divided into five sections, the first being this introduction and the last reserved for concluding remarks. Section 2 introduces a novel descriptive approach for lumped-parameter dynamic systems based on the use of hierarchies of mathematical models. Section 3 presents the fundamental theoretical results of this paper concerning the general structure of a recursive modular modelling algorithm. Section 4 develops a recursive algorithm based on the Udwadia–Kalaba equation and proves that it is identical to the algorithm of a Kalman filter for the estimation of the state of a static process subjected to a sequence of noiseless measurements representing the constraint equations. See electronic supplementary material, appendix S1, where the steps required for modelling a planar 3RRR parallel mechanism applying the proposed recursive modular methodology are presented in detail. This elucidating example highlights the main advantages of applying the proposed methodology for modelling lumped-parameter dynamic systems discussed in this text.

The derivation of the results presented in this paper is based on classical definitions and theorems from Linear Algebra and on the theory of generalized inverses. The reader is invited to check a textbook on generalized inverses by Campbell & Meyer [18], as well as their applications to Analytical Dynamics in the textbook by Udwadia & Kalaba [3]. An n×m matrix A^g is said to be a generalized inverse of an m×n matrix A if it satisfies at least the first of the following properties:

If A^g satisfies property P₁, it is said to be a {1}-inverse of A; if both P₁ and P₄ are satisfied, it is a {1,4}-inverse of A; if P₁, P₂ and P₄ are satisfied, A^g is a {1,2,4}-inverse of A and so on. There is a unique generalized inverse matrix that satisfies the four properties simultaneously: it is called the Moore–Penrose pseudoinverse of A and is specially denoted by A⁺ [3,18].

Finally, it is worth commenting that the notation adopted in this text is not only compatible with the desired generality for the methodology proposed, but also tries to resemble as much as possible the notations from previous works from the author [1,2] and from texts on the Udwadia–Kalaba equation and on the Kalman filter. Once all the variables in the presented equations are matrices, no distinctive boldface notation is adopted. Lowercase letters are reserved for column matrices which, in this text, are assumed to be equivalent to tuples. The exceptions are the letter t used to denote time (scalar) and the lowercase letters used as subscript or superscript indexes. Uppercase notation is reserved for matrices that may not be representable by a single row or column of complex entries. The identity matrices are denoted by I and the zero matrices and zero column matrices by 0. Script letters ($\mathcal{A},\mathcal{B},\dots$) are used to denote dynamic systems. Subscript indexes refer to the level of the hierarchical representation to which the given matrix or column matrix refer. Matrices and column matrices which are the same for all levels do not receive subscripts. Underline notation is used to denote functions. Superscripts (⋅)*, (⋅)^g and (⋅)⁺ denote, respectively, the conjugate transpose, a generalized inverse and the Moore–Penrose pseudoinverse of a given matrix. The index superscripts (⋅)^〈k〉 and (⋅)^{〈〈k〉〉} are associated with the order of a given set of generalized variables or of functions of them and are properly defined in §2. Table 1 presents the list of symbols adopted in this text.

2. Hierarchical description of lumped-parameter dynamic systems

This section proposes a special form of conceiving lumped-parameter dynamic systems. Assume that a system is described by a hierarchy of mathematical models of increasing complexity, such that the hierarchy corresponds to the equations of motion that better model the dynamic behaviour of the system. Consider also that all the models within this hierarchy are expressed in terms of the same set of variables, and that the increase in complexity from one level to its superior is associated with the enforcement of constraints.

A multibody system, for instance, is a lumped-parameter system which can be conceived as a set of assembled mechatronic modules, each of them being either a single body or a simpler multibody system. In order to exemplify how to describe a multibody system as a hierarchy of mathematical models, consider the 3RRR parallel mechanism (system $\mathcal{M}$) illustrated in figure 2. An intuitive form of describing this mechanism is to conceive it as a set of four assembled modules: ${\mathcal{H}}_{\mathtt{A}}$, ${\mathcal{H}}_{\mathtt{B}}$ and ${\mathcal{H}}_{\mathtt{C}}$ representing active RR kinematic chains and $\mathcal{L}$ representing its platform (end-effector). If the mathematical models of these four subsystems are already known, a model for the 3RRR parallel mechanism can be directly obtained. If not, as illustrated in figure 3, each of the three active RR kinematic chains, generically denoted by ${\mathcal{H}}_{\mathtt{K}}$, can be further partitioned, for example, in three modules: ${\mathcal{A}}_{\mathtt{K}}$ representing an actuator, ${\mathcal{U}}_{\mathtt{K}}$ representing a pivoted rigid bar and ${\mathcal{B}}_{\mathtt{K}}$ representing a generic two-dimensional bar element. In this latter case, the 3RRR parallel mechanism can be conceived as a hierarchy of three levels of mathematical models as shown in figure 4. At level 1, the model is constituted by a system of DAEs describing the decoupled dynamics of the modules ${\mathcal{A}}_{\mathtt{K}}$, ${\mathcal{U}}_{\mathtt{K}}$, ${\mathcal{B}}_{\mathtt{K}}$ (K=A,B,C) and $\mathcal{L}$. The modelling at this level consists of simply gathering the systems of equations of motion associated to each of these modules, neglecting all the constraints among them, which means that these systems of equations are decoupled from each other. Also, the set of variables in terms of which the models at superior levels are described is given by the union of the sets of variables adopted for the modelling of each of these modules. At level 2, the model corresponds to the system of DAEs describing the decoupled dynamics of the modules ${\mathcal{H}}_{\mathtt{A}}$, ${\mathcal{H}}_{\mathtt{B}}$, ${\mathcal{H}}_{\mathtt{C}}$ and $\mathcal{L}$, which can be obtained from the model at level 1 by enforcing the constraints associated to the assembly of each ${\mathcal{H}}_{\mathtt{K}}$ (i.e. the constraints among the modules ${\mathcal{A}}_{\mathtt{K}}$, ${\mathcal{U}}_{\mathtt{K}}$ and ${\mathcal{B}}_{\mathtt{K}}$). Note that such an approach allows to explore the topological symmetries among the three active RR kinematic chains in the application of the constraint enforcement algorithm, such that the systems of equations of motion associated to the modules ${\mathcal{H}}_{\mathtt{A}}$, ${\mathcal{H}}_{\mathtt{B}}$, ${\mathcal{H}}_{\mathtt{C}}$ become similar and remain decoupled from each other, as well as from the DAEs associated to the module $\mathcal{L}$. Finally, the enforcement of the constraints among these four modules leads to the mathematical model of the 3RRR parallel mechanism, which corresponds to the level 3 (the hierarch) of this hierarchy.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

Denote by $\mathcal{S}$ a generic lumped-parameter system and assume it is described as a hierarchy of mathematical models satisfying the conditions stated above. These models are constituted by parameters, generalized variables and DAEs.

The parameters represent a finite set of constant or known time-variant quantities that describe, according to the adopted simplifying assumptions, the properties of each component of a system. Constant or known time-variant quantities within the expressions of constitutive equations are also parameters of the system. In the case of a multibody system, parameters describe physical properties such as geometry, inertia, stiffness, damping, etc.

The generalized variables can be understood as a finite indexed family of variables which, along with the parameters, is able to provide a description of any quantity that can possibly be associated to a system. If x stands for any quantity related to the description of configurations⁴ of $\mathcal{S}$, there must be a tuple of variables q^〈0〉, called generalized coordinates, such that x can be expressed as a function of q^〈0〉 and time (along with the parameters), i.e. $x=\underset{\_}{x}(t,q^{⟨0⟩})$. If v stands for any quantity associated to some rate of variation of configuration,⁵ some further variables, typically called quasi-velocities, are required for the parametric description of this quantity. Denoting the tuple of quasi-velocities of the system $\mathcal{S}$ by q^〈1〉 and by q^{〈〈1〉〉}=(q^〈0〉,q^〈1〉), it can be stated that $v=\underset{\_}{v}(t,q^{⟨0⟩},q^{⟨1⟩})=\underset{\_}{v}(t,q^{⟨⟨1⟩⟩})$. A trivial choice of quasi-velocities is to take $q^{⟨1⟩}={\dot{q}}^{⟨0⟩}$. Higher-order generalized variables can be defined, recursively, in a similar way. Let q^{〈〈k〉〉}=(q^〈0〉,q^〈1〉,…,q^〈k〉). Formally, define the (k+1)-th order generalized variables of a given system, q^〈k+1〉, as a finite indexed family of variables such that: (i) q^〈k+1〉 can be expressed as a function of $(t,q^{⟨⟨k⟩⟩},{\dot{q}}^{⟨k⟩})$; (ii) ${\dot{q}}^{⟨k⟩}$ can be expressed as a function of (t,q^{〈〈k〉〉},q^〈k+1〉)=(t,q^{〈〈k+1〉〉}). We say that q^〈k+1〉 is trivially defined if $q^{⟨k+1⟩}={\dot{q}}^{⟨k⟩}$.

If a generic quantity is associated to a module of $\mathcal{S}$, which, at some level of the hierarchy, has its dynamics decoupled from other modules, then, at that level, the description of this quantity requires only the subset of generalized variables associated to the particular module. Also, a given set of generalized variables is said to be redundant at a given level of the hierarchy if it has to satisfy some conditions in order not to violate any of the constraints imposed to the system at that level. These conditions, when given by equalities, are called constraint equations. At a given level r of the hierarchy they might be expressed as follows:

The superscript k of ${\underset{\_}{h}}_r^{⟨k⟩}$ denotes that the corresponding constraint equations can be expressed as a function of time and of the generalized variables of $\mathcal{S}$ up to k-th order. Again, if some of these constraints correspond to a set of modules which, up to level r, have their dynamics decoupled from some others, the generalized variables originally associated to the latter modules shall not be present in these constraint equations (2.1). If the expressions in ${\underset{\_}{h}}_r^{⟨k⟩}$ are integrable forms, the constraint equations can alternatively be expressed as a function of time and generalized variables up to (k−1)-th order, i.e. ${\underset{\_}{h}}_r^{⟨k-1⟩}(t,q^{⟨⟨k-1⟩⟩})=0$. The time derivative of ${\underset{\_}{h}}_r^{⟨k⟩}$, on the other hand, can be expressed as a function of time and generalized variables up to (k+1)-th order, i.e. ${\underset{\_}{h}}_r^{⟨k+1⟩}(t,q^{⟨⟨k+1⟩⟩})=0$, which is affine with respect to the higher-order variables.

Considering this, choose an integer k such that the following conditions are satisfied.

• (i) All the generalized variables of $\mathcal{S}$ above k-th order are trivially defined, i.e. $q^{⟨l+1⟩}={\dot{q}}^{⟨l⟩}$ for all l≥k.

• (ii) All the constraint equations at any level of the hierarchy can be expressed in the following affine form:

  $${\underset{\_}{A}}_r(t,q^{⟨⟨k-1⟩⟩})q^{⟨k⟩}={\underset{\_}{b}}_r^{⟨k-1⟩}(t,q^{⟨⟨k-1⟩⟩}).$$2.2

• (iii) The dynamic equations of motion associated to any level of the hierarchy can be expressed in the following form:

  $$\underset{\_}{M}(t,q^{⟨⟨k-1⟩⟩})q^{⟨k⟩}={\underset{\_}{f}}^{⟨k-1⟩}(t,q_{\mathcal{N}}^{⟨⟨k-1⟩⟩})+{\gamma }_r,$$2.3
  with γ_r representing the generalized constraint forces acting in $\mathcal{S}$ at the level r of the hierarchy (i.e. γ_r stands for generalized forces due to all the constraints enforced at that level).

In order to obtain a complete system of DAEs of motion associated to a given level r of the hierarchy representing system $\mathcal{S}$, it is necessary to establish some assumptions for the computation or elimination of the generalized constraint forces from equation (2.3). Assume that there is a matrix $S_r={\underset{\_}{S}}_r(t,q^{⟨⟨k-1⟩⟩})$ such that ${\underset{\_}{S}}_r^{\ast }(t,q^{⟨⟨k-1⟩⟩}){\gamma }_r$ is equal to a known function of time and of the generalized variables of $\mathcal{S}$ up to (k−1)-th order, which will be denoted by ${\underset{\_}{S}}_r^{\ast }(t,q^{⟨⟨k-1⟩⟩})\underset{\_}{\phi }(t,q^{⟨⟨k-1⟩⟩})$, following Udwadia & Kalaba [5,6]. If $\mathcal{S}$ respects the extended version of D’Alembert’s principle⁶ [7,19], $\underset{\_}{\phi }(t,q^{⟨⟨k-1⟩⟩})=0$ and $S_r={\underset{\_}{S}}_r(t,q^{⟨⟨k-1⟩⟩})$ must be a matrix whose columns span the kernel of $A_r={\underset{\_}{A}}_r(t,q^{⟨⟨k-1⟩⟩})$, i.e. $\mathrm{im}(S_r)=\ker (A_r)$. For the cases in which the total virtual work of the generalized constraint forces is not zero, i.e. in which a constraint is associated to an energy source or sink, it is assumed that this effect can be expressed by a constitutive relation defined by the function $\underset{\_}{\phi }(t,q^{⟨⟨k-1⟩⟩})$ [6].

A complete system of DAEs of motion for $\mathcal{S}$, associated to the level r of the hierarchy, can be expressed as follows⁷:

Therefore, given a state (t,q^{〈〈l+1〉〉}), the values obtained for q^〈k〉 and, consequently, for the time derivatives of the state variables will depend on the constraints enforced at the particular level of the hierarchy to which equations (2.4) correspond. Thus, renaming q^〈k〉 as α_r in these equations, the algebraic equations of this system can be rewritten in the following matrix form: A necessary condition for ξ_r to be uniquely determined is that L_r must be a full-rank matrix. Provided that the ordinary differential equations (ODEs) ${\dot{q}}^{⟨l⟩}={\underset{\_}{\dot{q}}}^{⟨l⟩}(t,q^{⟨⟨l+1⟩⟩})$, l=0,1,…,k−1, satisfy the hypotheses of Picard’s existence and uniqueness theorem, the uniqueness of ξ_r implies the uniqueness of solutions of the equations of motion associated to the level r of the hierarchy.

It is worth noting that once the generalized variables are the same for all the levels of the hierarchy, so must be the matrix M and the column matrix f. Indeed, the entries of these matrices are associated to the chosen set of variables and can be obtained already at the first (lowest) level of the hierarchy, which means that their expressions are subject to the same modular decoupling present at this first level. Back to the example of the 3RRR parallel mechanism, considering the hierarchy presented at figure 4,⁸ it can be stated that M is a block-diagonal matrix whose blocks are the generalized inertia matrices associated to the models of the modules $\mathcal{L}$, ${\mathcal{A}}_{\mathtt{K}}$, ${\mathcal{U}}_{\mathtt{K}}$ and ${\mathcal{B}}_{\mathtt{K}}$ (K=A,B,C), and f is a column matrix in which the rows associated to the generalized variables of a given module (among these 10) provide the expressions of the generalized active and gyroscopic forces corresponding to the respective module in terms of its own related variables.

Moreover, such a modular structure might also be noticed in the expressions of the matrices A_r. If a block of rows of A_r describes a constraint among a given set of modules of a system, only the columns corresponding to the variables associated to these modules have non-zero entries. In the example of the 3RRR parallel mechanism, a given row of A₁ (which is related to the constraints at level 1) cannot have non-zero entries in two columns associated to different modules among the 10 at this level; similarly, a given row of A₂ (related to level 2) cannot have non-zero entries in two columns associated to different modules among $\mathcal{L}$, ${\mathcal{H}}_{\mathtt{A}}$, ${\mathcal{H}}_{\mathtt{B}}$ and ${\mathcal{H}}_{\mathtt{C}}$. An analogous reasoning applies to the entries of the corresponding column matrices b_r.

Such a block-diagonal structure of the matrices related to the adoption of a modular description for multibody systems has already been explored in previous methodologies proposed by Orsino [1] and Orsino & Hess-Coelho [2]. As previously mentioned, the developments shown in this paper can be understood as a recursive generalization of these methodologies for arbitrary lumped-parameter dynamic systems.

3. General recursive algorithm for obtaining explicit equations of motion of constrained lumped-parameter dynamic systems

As previously discussed, it is possible to represent any lumped-parameter dynamic system as a hierarchy of mathematical models organized in order of increasing complexity, such an increase being associated to the enforcement of some constraints from one level of the hierarchy to its superior level. In such a representation, the hierarch is assumed to be the most precise mathematical model for the system. In the previous section, it was shown that a mathematical model at any level of such a hierarchy can be expressed by a system of DAEs constituted by a subsystem of ODEs and a subsystem of linear algebraic equations (in the corresponding variables α_r and γ_r), whose coefficients are functions of the state of the system.

This section concentrates on presenting the fundamentals for the development of a general recursive algorithm for obtaining explicit equations of motion for each level of a hierarchy that represents a lumped-parameter dynamic system. This algorithm traverses these hierarchies from the lowest level (which is indexed as level 1) to the top one (level s). Let the subsystem of linear algebraic equations associated to a given level r of a hierarchy be represented by the following matrix equations:

Pre-multiplying equation (3.1) by $S_r^{\ast }$ and using the identity in (3.3), it can be stated that α can be directly obtained as a solution of the following equation: This procedure can be interpreted as a generalization of the idea of projecting the components of the system of forces (assuming that the inertial effects are conceived as inertial forces) along directions in which the components of the constraint forces are known. In this case, $S_r^{\ast }$ represents an operator that performs this projection.

Consider the following identity valid for any matrices A and B:

According to this property, it is also possible to rewrite equation (3.4) as follows: Therefore, it is always possible to adopt a Hermitian operator that performs these projections.

Assuming that the columns of S_r span the kernel of A_r, as required by the extended version of D’Alembert’s principle, the following proposition presents some alternatives for computing S_r.

Proposition 3.1.

Let A be an arbitrary m×n matrix with complex entries, such that$\ker (A)\ne \{0\}$. It can be stated that the following statements are true.

• (i) There is a full-rank matrix Z defining a linear transformation onto$\ker (A)$. In this case, Z is shortly referred to as an orthogonal complement⁹of A.

• (ii) If A^gdenotes a {1}-inverse of A, Y =I−A^gA is a projector onto$\ker (A)$.

• (iii) If A^gdenotes a {1,4}-inverse of A, X=I−A^gA is the orthogonal projector onto$\ker (A)$.¹⁰

• (iv) Let$\hat{A}$be a matrix obtained by an orthonormalization procedure applied to the rows of A. The orthogonal projector onto$\ker (A)$can be alternatively computed by the expression$X=I-{\hat{A}}^{\ast }\hat{A}$.

Proof.

Once $\ker (A)\ne \{0\}$, there are non-zero $z\in {\mathbb{C}}^n$ such that Az=0. Also, as $\ker (A)$ is a linear subspace of ${\mathbb{C}}^n$, if $z\in \ker (A)$, $\kappa z\in \ker (A)$ for any $\kappa \in \mathbb{C}$.

• (i) Take any basis for $\ker (A)$, i.e. take a finite set of linearly independent $z_i\in \ker (A)$ such that any $z\in \ker (A)$ can be expressed as a linear combination of them. By definition, a matrix Z whose columns are given by these z_i is a full-rank matrix that defines a linear transformation onto $\ker (A)$, i.e., Z is an orthogonal complement of A.

• (ii) Denote by A^g a {1}-inverse of A and let Y =I−A^gA. For any $w\in {\mathbb{C}}^n$, $y=Yw\in \ker (A)$. Indeed, AY =A−AA^gA=0. Also, YY =I−A^gA−A^gA+A^gAA^gA=I−A^gA=Y . If $y=Yw\in \ker (A)$, $Yy=YYw=Yw=y\in \ker (A)$ and if Yy=y, Ay=AYy=0, which means that $y\in \ker (A)$. Therefore, Y is a projector onto $\ker (A)$.

• (iii) Denote by A^g a {1,4}-inverse of A and let X=I−A^gA. From property P₄, equation (1.1), it follows that X is Hermitian. Also, from the uniqueness of the orthogonal projector onto a given linear space,¹¹ it follows that X is the orthogonal projector onto $\ker (A)$.

• (iv) Let $\hat{A}$ be a matrix obtained by an orthonormalization procedure applied to the rows of A and define $X=I-{\hat{A}}^{\ast }\hat{A}$. Note that $\hat{A}$ is a full-rank matrix and its rows constitute a basis for the linear subspace spanned by the rows of A, i.e. $\mathrm{im}({\hat{A}}^{\ast })=\mathrm{im}(A^{\ast })$. Thus, there exists a matrix U such that $A=U\hat{A}$. Also, $AX=U\hat{A}-U\hat{A}{\hat{A}}^{\ast }\hat{A}$, and once $\hat{A}$ was obtained by an orthonormalization procedure, $\hat{A}{\hat{A}}^{\ast }=I$, which means that AX=0. In this case, X is not only a Hermitian matrix, but also a projector, once $XX=I-{\hat{A}}^{\ast }\hat{A}-{\hat{A}}^{\ast }\hat{A}+{\hat{A}}^{\ast }\hat{A}{\hat{A}}^{\ast }\hat{A}=I-{\hat{A}}^{\ast }\hat{A}=X$. For any $w\in {\mathbb{C}}^n$, $x=Xw\in \ker (A)$. If $x=Xw\in \ker (A)$, $Xx=XXw=Xw=x\in \ker (A)$ and if Xx=x, Ax=AXx=0, which means that $x\in \ker (A)$. Therefore, $X=I-{\hat{A}}^{\ast }\hat{A}$ is the orthogonal projector onto the kernel of A.

This completes the proof. ▪

If S_r is chosen to be an orthogonal complement of A_r, equation (3.4) becomes the one addressed in previous developments of a modular modelling metodology for multibody systems by Orsino [1] and Orsino & Hess-Coelho [2].¹² On the other hand, if $S_r=I-A_r^{+}{A}_r$, with $A_r^{+}$ denoting the Moore–Penrose pseudoinverse¹³ of A_r, equation (3.4) becomes the one addressed by Udwadia & Phohomsiri [4] in their derivation of explicit equations of motion for constrained mechanical systems with singular mass matrices (which led to applications also to the modular modelling of multibody systems).¹⁴

Consider also that the constraints specifically reinforced at level r of this hierarchy can be expressed by the equation ${\tilde{A}}_r{\alpha }_r={\tilde{b}}_r$. In this case, $A_1={\tilde{A}}_1$, $b_1={\tilde{b}}_1$ and the indexed families of matrices {A₁,A₂,…,A_r,…,A_s} and of column matrices {b₁,b₂,…,b_r,…,b_s} must satisfy the following properties, for r=1,2,…,s−1:

In order to derive recursive algorithms for computing S_r+1 given S_r, the results stated in the following theorem are needed.

Theorem 3.2.

Let {A₁,A₂,…,A_r,…,A_s} be an indexed family of matrices with complex entries satisfying the following condition:

Let S_rdenote a linear operator onto the kernel of A_r. It can be stated that if C_r+1defines a linear operator onto the kernel of$B_{r+1}={\tilde{A}}_{r+1}{S}_r$, then S_r+1=S_rC_r+1is a linear operator onto the kernel of A_r+1.

Proof.

Owing to condition (3.8), it can be stated that: $\ker (A_{r+1})\subset \ker (A_r)\subset {\mathbb{C}}^n$. Denote by S_r a linear operator onto $\ker (A_r)$ and assume that S_r+1=S_rC_r+1 is a linear operator onto $\ker (A_{r+1})$. It can be stated that

Therefore, if C_r+1 defines a linear transformation onto $\ker ({\tilde{A}}_{r+1}{S}_r)$, then S_r+1=S_rC_r+1 defines a linear transformation onto $\ker (A_{r+1})$. ▪

This theorem proves that, given an operator S_r onto $\ker (A_r)$ for a level r of a hierarchy representing a lumped-parameter dynamic system, the operator associated to the superior level, i.e. an operator S_r+1 onto $\ker (A_{r+1})\subset \ker (A_r)$, can be obtained as a composition of S_r with an operator C_r+1 onto $\ker ({\tilde{A}}_{r+1}{S}_r)$. It is worth noting that $B_{r+1}={\tilde{A}}_{r+1}{S}_r$ represents the operator corresponding to the new constraints to be enforced at level r+1, ${\tilde{A}}_{r+1}$, restricted to the kernel of A_r, in which the variations of any solution of the equations of motion associated to level r must lie.

Particularly, when the operator S_r onto $\ker (A_r)$ is chosen to be the orthogonal projector onto this subspace, the corollary below follows from theorem 3.2.

Corollary 3.3.

Let {A₁,A₂,…,A_r,…,A_s} be an indexed family of matrices satisfying the condition (3.8).

Let$S_r=I-A_r^{\mathrm{g}}{A}_r$, for some r, be the orthogonal projector onto the kernel of A_r, with (⋅)^gdenoting a {1,4}-inverse. It can be stated that

with$G_{r+1}={({\tilde{A}}_{r+1}{S}_r)}^{\mathrm{g}}$, is the orthogonal projector onto the kernel of A_r+1and satisfies$S_{r+1}=I-A_{r+1}^{\mathrm{g}}{A}_{r+1}$.

Proof.

From proposition 3.1, once $A_r^{\mathrm{g}}$ stands for a {1,4}-inverse of A_r, $S_r=I-A_r^{\mathrm{g}}{A}_r$ is the orthogonal projector onto $\ker (A_r)$. Moreover, let $C_{r+1}=I-{({\tilde{A}}_{r+1}{S}_r)}^{\mathrm{g}}{\tilde{A}}_{r+1}{S}_r$ be the orthogonal projector onto $\ker ({\tilde{A}}_{r+1}{S}_r)$, which satisfies the statement of theorem 3.2. Taking S_r+1=S_rC_r+1, it can be stated that

Once $S_{r+1}=S_{r+1}{S}_{r+1}^{\ast }$, S_r+1 is a Hermitian matrix. Also, replacing $S_{r+1}^{\ast }$ by S_r+1, it follows that S_r+1=S_r+1S_r+1, which proves that S_r+1 is also idempotent. Thus, S_r+1 is the orthogonal projector onto $\ker (A_{r+1})$. Once $S_{r+1}=S_{r+1}^{\ast }=C_{r+1}^{\ast }{S}_r^{\ast }=C_{r+1}{S}_r$, it can be stated that The second equality in equation (3.10) follows from $S_{r+1}=S_{r+1}{S}_{r+1}^{\ast }$, replacing S_r+1 in the right-hand side of this identity by the expression obtained in equation (3.12). Finally, the fact that $S_{r+1}=I-A_{r+1}^{\mathrm{g}}{A}_{r+1}$ follows from the uniqueness of the orthogonal projector onto a given linear space. ▪

A second corollary that provides a recursive algorithm for obtaining {1,4}-inverses for an indexed family of matrices satisfying the condition (3.8), which will be useful in the derivations of the following section, is presented below.

Corollary 3.4.

Let {A₁,A₂,…,A_r,…,A_s} be an indexed family of matrices satisfying the condition (3.8) and denote by (⋅)^ga {1,4}-inverse. It can be stated that, given a {1,4}-inverse of A_r, a {1,4}-inverse of A_r+1can be computed as follows:

with$G_{r+1}={({\tilde{A}}_{r+1}{S}_r)}^{\mathrm{g}}$.

Proof.

Observing the statements of P₁ and P₄ in equation (1.1), one could reinterpret them as properties of the matrix A^gA, without any loss in generality. Thus, in order to prove that $[(I-G_{r+1}{\tilde{A}}_{r+1})A_r^{\mathrm{g}}\hspace{0.17em}{G}_{r+1}]$ is a {1,4}-inverse of A_r+1, it is enough to show that $[(I-G_{r+1}{\tilde{A}}_{r+1})A_r^{\mathrm{g}}\hspace{0.17em}{G}_{r+1}]A_{r+1}=A_{r+1}^{\mathrm{g}}{A}_{r+1}$ for any {1,4}-inverse of A_r+1, $A_{r+1}^{\mathrm{g}}$. Let $A_r^{\mathrm{g}}{A}_r=I-S_r$ and $A_{r+1}^{\mathrm{g}}{A}_{r+1}=I-S_{r+1}$. According to equation (3.10), it can be stated that $G_{r+1}{\tilde{A}}_{r+1}{S}_r=S_r-S_{r+1}$. Therefore,

which completes the proof. ▪

It is worth noting that even if $A_r^{\mathrm{g}}$ is the Moore–Penrose pseudoinverse of A_r, the algorithm described by equation (3.13) only ensures that $A_{r+1}^{\mathrm{g}}$ will be the {1,4}-inverse of A_r+1. Particularly, if some row of ${\tilde{A}}_{r+1}$ can be expressed as a linear combination of the rows of A_r, $A_{r+1}^{\mathrm{g}}$ computed according to corollary 3.4 will not be the Moore–Penrose pseudoinverse of A_r+1 [21]. On the one hand, it might seem useful to work with Moore–Penrose pseudoinverses only, once they are uniquely defined; on the other, there are no practical disadvantages on adopting generic {1,4}-inverse matrices in the constraint enforcement algorithms as it is shown in the next section, which emphasizes the importance of corollary 3.4.

4. Recursive algorithm based on the Udwadia–Kalaba equation

Assume that M is a positive-definite Hermitian matrix¹⁵ and note that equation (3.3) can be rewritten in the following form: $S_r^{\ast }({\gamma }_r-\phi )=0$. This last fact means that γ_r−φ lies in $\ker (S_r^{\ast })=\mathrm{im}(A_r^{\ast })$, i.e.:

for some λ_r. Thus, according to equation (3.1), it can be stated that with α₀=M⁻¹(f+φ). Moreover, replacing (4.2) into equation (3.2), the following expressions for λ_r are obtained, with (⋅)^g standing for a {1,4}-inverse and w_r for an arbitrary column matrix: and Considering that M is Hermitian, so must be M^−1/2. Thus, defining H_r=A_rM^−1/2, it can be stated that $A_r{M}^{-1}{A}_r^{\ast }=H_r{H}_r^{\ast }$ and that $M^{-1}{A}_r^{\ast }=M^{-1/2}{H}_r^{\ast }$. Therefore, using (4.4), equation (4.2) can be rewritten as follows: This last equation can be further simplified by considering the following identities, valid for (⋅)^g denoting {1}-inverses: and Indeed, according to property P₁, equation (1.1): (HH*(HH*)^g−I)HH*=0. Using the conjugate transpose of equation (3.5), it can be stated that (HH*(HH*)^g−I)HH*=0 ⇔ (HH*(HH*)^g−I)H=0 ⇔ H(H*(HH*)^g)H=H. This proves identity (4.6). Also, due to property P₁, equation (1.1): HH*(I−(HH*)^gHH*)=0; therefore, again due to equation (3.5), H*(I−(HH*)^gHH*)=0, which proves identity (4.7). Thus, equation (4.5) can be simplified to the following form: Moreover, adopting the transformation of variables α_r=M^−1/2a_r and α₀=M^−1/2a₀, equation (4.8) can be ultimately expressed in the form This is the classical Udwadia–Kalaba equation [5,8].

In order to derive a recursive modular algorithm based on the Udwadia–Kalaba equation, let ${\tilde{H}}_{r+1}={\tilde{A}}_{r+1}{M}^{-1/2}$. According to corollary 3.4, it can be stated that

with $K_{r+1}={({\tilde{H}}_{r+1}{P}_r)}^{\mathrm{g}}$ and $P_r=I-H_r^{\mathrm{g}}{H}_r$. Also, according to corollary 3.3, $P_{r+1}=I-H_{r+1}^{\mathrm{g}}{H}_{r+1}$ can be computed as follows: Note that $P_r=I-H_r^{\mathrm{g}}{H}_r$ is an idempotent Hermitian matrix. Whenever H_r+1 is a full-rank matrix, so that the rows of ${\tilde{H}}_{r+1}{P}_r$ are linearly independent, it can be stated that K_r+1 can be computed as follows: These considerations lead to the following algorithm, whose inputs are a₀, P₀, and all the ${\tilde{H}}_{r+1}$:

It is noticeable that, in this form, the algorithm for updating the values of a_r according to the information provided by the constraint equations at each level of the hierarchical structure that represents a dynamic system is identical to the algorithm of a Kalman filter for the estimation of the state of a static process subjected to a sequence of measurements modelled by the noiseless observation equation ${\tilde{b}}_{r+1}={\tilde{H}}_{r+1}a$. Matrix P_r, which represents the orthogonal projector into $\ker (H_r)$, being directly related to the description of the constraints of the system, is analogous to the estimate covariance matrix.

The analogy between this recursive algorithm based on the Udwadia–Kalaba equation and filtering algorithms can be well explored, both in terms of the use of existing computational implementations for simulations of the models of lumped-parameter dynamic systems and in terms of the development of new modelling and simulation approaches that may arise from this correspondence. For example, it might be possible to treat constraint clearances as if they were measurement noise. In this case, if R_r represents the covariance matrix associated to the noise variables, the corresponding K_r+1 and P_r+1 can be obtained by the following expressions:

and

5. Conclusion

This paper proposes a novel approach for the modelling of lumped-parameter dynamic systems, which consists of representing a given system by a hierarchy of mathematical models, expressed in terms of the same set of variables, and ordered according to the increase in complexity due to the enforcement of constraints. In such representation, the hierarch stands for the model that better represents the dynamic behaviour of the system. This approach allows to better explore the modularity of the systems if, for instance, it is possible to conceive them as a set of constrained subsystems whose mathematical models are already known. This was the context for which previous developments of similar methodologies for the modelling of multibody systems, by Orsino [1] and Orsino & Hess-Coelho [2], were proposed. They can be understood as particular cases of the general approach discussed in this text. Also, whenever a refinement of some model can be associated with the enforcement of further constraints, the methodology presented in this paper also provides a guideline for a modelling procedure. This might be the case for lumped-parameter approximations for distributed systems, for example, in which the number of levels of the hierarchy could be understood as the number of iterations required for the model to adequately represent, according to some established criteria, the original system.

The methodology presented in this text is based on the development of recursive algorithms for finding projection operators onto the kernel of the ones associated to the constraints at a particular level of a hierarchy representing a dynamic system. The modularity of the systems can also be explored in these algorithms once, particularly at the lower levels of a hierarchy, operators associated to the constraints and, consequently, the ones that make projections onto their kernels might be represented by block-sparse or even block-diagonal matrices. Exploring properly this fact, numerically efficient algorithms can be implemented. Within the recursive approach proposed, there is only need to obtain projection operators associated to new constraints that are enforced from one level to the superior one. The algorithms used for this can be based on finding orthogonal complement matrices, on computing generalized inverses or even on implementing orthonormalization procedures.

A remarkable result that came along with the development of this new methodology is the equivalence between the formulations by Orsino [1] and Orsino & Hess-Coelho [2] and the formulation by Udwadia & Phohomsiri [4]. The investigation concerning the derivation of a recursive algorithm based on the Udwadia–Kalaba equation, under the same conditions adopted for the general methodology, revealed an even more relevant result: the algorithm is equivalent to the one of a Kalman filter for the estimation of the state of a static process in which the constraint equations play the role of the observation ones, and the associated orthogonal projectors play the role of the covariance matrices. This allows not only to reuse already existing routines in which the Kalman filter algorithm is implemented for forward simulations¹⁶ of lumped-parameter dynamic systems, but also to effectively make use of the analogies to perform, for example, simulations in which constraint clearances are considered. Another aspect that should be explored in future works is how the interpretations of the Kalman filter within the theory of stochastic processes can be helpful in enhancing the framework of methodologies for the modelling of dynamic systems.

Data accessibility

This article has no additional data.

Competing interests

I declare I have no competing interests.

Funding

I acknowledge the post-doctoral grant no. 2016/09730-0, São Paulo Research Foundation (FAPESP).

Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3849613.