A Novel Mixed Finite/Infinite Dimensional Port-Hamiltonian Model of a Mechanical Ventilator

1. Introduction

The World Health Organization (WHO) has "access to oxygen" on its model list as one of the essentials required by an individual, especially when they are in a critical condition health-wise and they are unable to breathe on their own. Access to oxygen is the only medicine listed that does not have a substitute. Access to oxygen even after the worldwide pandemic is still a major dilemma in middle and low-income countries [1]. A mechanical ventilator is a medical device that is employed in assisting patients in cases where their respiratory system is not functioning well, thus the patient has challenges in breathing or has shortness of breath. The mechanical ventilator can also be used in cases where a patient has been sedated and is undergoing surgery.

The basic operation of a mechanical ventilator is to control a high-pressure region during the inspiration stage. In the inspiration stage, the mechanical ventilator is paused. During expiration, air flows out due to the lung’s natural recoil which creates a higher pressure in the alveoli. Due to the high demand for mechanical ventilators both in the past and at present, it is imperative that efficient and affordable mechanical ventilators should be researched, modelled, designed and implemented. To ensure that robust mechanical ventilators are designed, it is important to formulate new models that can be used in the research and testing stages of the mechanical ventilators.

2. Related works

Tran et al, [2], conducted research whose main goal was to design, control, model, and simulate a mechanical ventilator that is light in weight, portable, and suitable for use at home. [2], modelled the mechanical ventilator as a voltage source. El-Hadj et al [3], applied the fluid-structure Interaction (FSI) to a couple of computational fluid dynamics used for fluid flow with finite element analysis for the solid domain. This facilitated the investigation of fluid behaviour, structural behavior, and interactions. Two designs were proposed for the mechanical ventilator. The first design was reported to be uncontrollable and the second design considers Computational Fluid Dynamics (CDF) with a moving boundary which is applied to the piston cylinder based ventilator. Tharton et al [4], designed and developed a ventilator prototype to be used by professionals in medical emergencies and in any other context where the regular ventilator is not available. The mechanical ventilator was modelled using the crank rocker mechanism in order to meet specific requirements for mechanical ventilator design. Pivk et al, [5], developed an empirical model for a low-cost mechanical ventilator by observing the response of each of the ventilator subsystems. The current progress made globally in the development of affordable and effective designs is important. However, research in mechanical ventilator development lacks model development. The Port-Hamiltonian approach applied in this research work acts as a modelling template for future energy-based mechanical ventilator modelling and design. It is important to develop designs that rest on a solid understanding of a significant aspect of the design. Naggar, [6] developed a piecewise model of a mechanical ventilator which described the artificial behaviour of a mechanical ventilator. A pressure-controlled ventilator is created and simulated. The mechanical ventilator is modelled using a periodic function with inequalities to control the beginning of inspiration and expiration periods. Shi et al, [7] developed a mathematical model of volume-controlled mechanical ventilation. The model is viewed as a pneumatic system where the ventilator is regarded as an air compressor. The exhalation valve is considered as throttle. The compressor and the container represent the ventilator. Naggar et al, [8] proposed an integrated mathematical model of the mechanical ventilator and the lung. Linear quadratic and exponential equations were used to model the system. The integrated model was used to simulate volume-controlled artificial ventilation. Giri et al, [9], proposed a simplified design of a mechanical ventilator, to reduce cost and automate the mechanical ventilation process. The proposed design was simulated on the MATLAB platform. Hannon et al, [10] presented a review on the advancement in the modelling of human anatomy, physiology and pathophysiology via mathematical modelling and computer simulation. Clinical applications in various disease states were emphasized. The research work discussed the current limitations and potential of in-silico modelling. There are currently no existing Port Hamiltonian models of mechanical ventilators intergrated with a human repiratory system. The main contribution in this research work is the development of an intergrated Port-Hamiltonian model representation of a mechanical ventilator-human respiratory system. The model consists of electromechanical and electromagnetic components modelled in the finite-dimensional representation, interconnected with fluid components in the infinite-dimensional representation. As a result, this Port-Hamiltonian model is of a mixed finite/infinite dimensional nature.

3. Materials and Method

This section presents the mathematical preliminaries of the Port Hamiltonian approach that have been followed in later sections. These definitions are standard definitions from literature and are mostly derived from the following references: [11,12] and [13]

Definition 1 (Dirac Structure)

A Dirac Structure (DS) is a pair of elements, $f\in {\mathbb{R}}^n$ and $e\in {\mathbb{R}}^n$, that satisfies the set:

$$\begin{array}{c}\hfill \mathcal{D}\left\{(f,e)\in \mathcal{D},(\widehat{f},\widehat{e})\in \mathcal{D}\left|e^{\mathsf{T}}\widehat{f}=-{\widehat{e}}^{\mathsf{T}}f\right.\right\}\end{array}$$

The DS is a subset $\mathcal{D}\subset \mathcal{F}\times \mathcal{E}$, where $\mathcal{F}$ and $\mathcal{E}$ represent the flows and efforts, respectively.

The DS in electric circuits can be represented as:

Proposition 1:

The subspace

$$\begin{array}{c}\hfill \mathcal{D}\subset {\mathbb{R}}^n\times {\mathbb{R}}^n\end{array}$$

is called a DS iff, there exist $A,B\in {\mathbb{R}}^{n\times n}$ such that

$$\begin{array}{c}\hfill A{B}^{\mathsf{T}}+A{B}^{\mathsf{T}}=0\mathrm{and}\mathrm{rank}\left(\left[AB\right]\right)=n\end{array}$$

satisfying the condition

$$\begin{array}{c}\hfill \mathcal{D}=\left\{(f,e)\in {\mathbb{R}}^n\times {\mathbb{R}}^n|Af=-Be\right\}\end{array}$$

Definition 2 (Resistive Relation)

Any relation $\mathcal{R}\subset {\mathbb{R}}^n\times {\mathbb{R}}^n$ is said to be resistive if $\forall (f_{\mathcal{R}},e_{\mathcal{R}})\in \mathcal{R}$ the

$$\begin{array}{c}\hfill e_{\mathcal{R}}^{\mathsf{T}}{f}_{\mathcal{R}}\le 0\end{array}$$

is satisfied.

Definition 3 (Port-Hamiltonian (pH) System).

The set $(\mathcal{D},\mathcal{L},\mathcal{R})$, defines a pH system where:

• $\mathcal{D}\subset ({\mathcal{F}}_{\mathcal{L}}\times {\mathcal{F}}_{\mathcal{R}}\times {\mathcal{F}}_{\mathcal{P}})\times ({\mathcal{E}}_{\mathcal{L}}\times {\mathcal{E}}_{\mathcal{R}}\times {\mathcal{E}}_{\mathcal{P}})$ is a DS
• $\mathcal{L}\subset {\mathcal{F}}_{\mathcal{L}}\times {\mathcal{E}}_{\mathcal{L}}$ is a LS and
• $\mathcal{R}\subset {\mathcal{F}}_{\mathcal{R}}\times {\mathcal{E}}_{\mathcal{R}}$ a Resistive relation.

The elements of the sets:

• ${\mathcal{F}}_{\mathcal{L}}\in {\mathbb{R}}^{n_{\mathcal{L}}}$ and ${\mathcal{E}}_{\mathcal{L}}\in {\mathbb{R}}^{n_{\mathcal{L}}}$ are known as flows and efforts
• ${\mathcal{F}}_{\mathcal{R}}\in {\mathbb{R}}^{n_{\mathcal{R}}}$ and ${\mathcal{E}}_{\mathcal{R}}\in {\mathbb{R}}^{n_{\mathcal{R}}}$ are known as resistive flows and efforts
• ${\mathcal{F}}_{\mathcal{P}}\in {\mathbb{R}}^{n_{\mathcal{P}}}$ and ${\mathcal{E}}_{\mathcal{P}}\in {\mathbb{R}}^{n_{\mathcal{P}}}$ are known as external flows and efforts, respectively.

where $n_{\mathcal{L}},n_{\mathcal{R}}$ and $n_{\mathcal{P}}\in {\mathbb{N}}_0$.

The dynamics of the pH system are given by the differential inclusion

$$\begin{array}{cc}\hfill & \hfill \left(-\frac{d}{dt}x\left(t\right),f_{\mathcal{R}}\left(t\right),f_{\mathcal{P}}\left(t\right),e_{\mathcal{L}}\left(t\right),e_{\mathcal{R}}\left(t\right),e_{\mathcal{P}}\left(t\right)\right)\in \mathcal{D}\hfill \end{array}$$

where $(x\left(t\right),e_{\mathcal{L}}\left(t\right))\in \mathcal{L}$, $(f_{\mathcal{R}}\left(t\right),e_{\mathcal{R}}\left(t\right))\in \mathcal{R}$ and $(f_{\mathcal{P}}\left(t\right),e_{\mathcal{P}}\left(t\right))\in \mathcal{P}$

Definition 4 (Interconnection of n pH Systems)

Let (${\mathcal{D}}_i,{\mathcal{L}}_i,{\mathcal{R}}_i$) denote the space of the $i^{\mathrm{th}}$ pH system in a set of n pH systems that are to be interconnected. The space of flows is divided into an external part and a part to be linked and is given by

$$\begin{array}{c}\hfill {\mathcal{F}}_i={\mathcal{F}}_{{\mathcal{L}}_i}\times {\mathcal{F}}_{\mathcal{R}i}\times {\mathcal{F}}_{\mathcal{P}i}\times {\mathcal{F}}_{\mathcal{P}\mathrm{link}}\end{array}$$

Similarly, the space of efforts is divided into an external part and a part to be linked and is expressed as

$$\begin{array}{c}\hfill {\mathcal{E}}_i={\mathcal{F}}_{\mathcal{L}i}\times {\mathcal{E}}_{\mathcal{R}i}\times {\mathcal{E}}_{\mathcal{P}i}\times {\mathcal{E}}_{\mathcal{P}\mathrm{link}}\end{array}$$

The interconnection of two pH systems $({\mathcal{D}}_1,{\mathcal{L}}_1,{\mathcal{R}}_1)$ and $({\mathcal{D}}_2,{\mathcal{L}}_2,{\mathcal{R}}_2)$ with respect to the link $(F_{\mathcal{P}\mathrm{link}},E_{\mathcal{P}\mathrm{link}})$ results in a new interconnected pH system, $(\mathcal{D},\mathcal{L},\mathcal{R})$. This interconnected system is given by the expression

$$\begin{array}{c}\hfill ({\mathcal{D}}_1,{\mathcal{L}}_1,{\mathcal{R}}_1)\circ ({\mathcal{D}}_2,{\mathcal{L}}_2,{\mathcal{R}}_2)(\mathcal{D},\mathcal{L},\mathcal{R})\end{array}$$

Definition 5 (Directed Graph).

A directed graph is a quadruple $\mathcal{G}(V,E,\mathrm{l},\mathrm{r})$ where

• V is a set of vertices
• E is a set of edges
• $\mathrm{l}:E\to V$ maps each edge, e, to an initial vertex
• $\mathrm{r}:E\to V$ maps each edge, e, to a terminal vertex

Definition 5.1 (Loop-free directed Graph).

If G is a directed graph, then G is said to be loop-free if for all $e\in E$,

$$\begin{array}{c}\hfill \mathrm{l}\left(e\right)\ne \mathrm{r}\left(e\right)\end{array}$$

Definition 5.2 (Subgraphs).

Given the graphs $\mathcal{G}(V,E,\mathrm{l},\mathrm{r})$ and ${\mathcal{G}}^{\prime }(V^{\prime },E^{\prime },\mathrm{l},\mathrm{r})$, then ${\mathcal{G}}^{\prime }$ is said to be a subgraph of ${\mathcal{G}}^{\prime }$ if $E^{\prime }\subset E$ and $V^{\prime }\subset V$. Furthermore,

• A subgraph is said to be an induced subgraph on $V^{\prime }$ if $E^{\prime }{=E|}_V^{\prime }$
• A subgraph is said to be spanning if $V^{\prime }=V$
• A subgraph is said to be a proper subgraph if $E^{\prime }\ne E$
• If both V and E are finite, then $\mathcal{G}$ is said to be finite

Definition 6 (Paths, Connectivity and Cycles).

Let $\mathcal{G}=(V,E,\mathrm{l},\mathrm{r})$ be a directed finite graph.

• An n-tuple $e=(e_1,\dots ,e_n)\in {(E\cup -E)}^n$ is called a path from $\upsilon$ to $\varpi$, if

  (a)

  $\mathrm{l}\left(e_1\right),\dots ,\mathrm{l}\left(e_n\right)$ are distinct

  (b)

  $\mathrm{r}\left(e_i\right)=\mathrm{l}\left(e_{i+1}\right)$ for all $i\in \left\{1,\dots ,n-1\right\}$

  (c)

  $\mathrm{l}\left(e_1\right)=\upsilon \wedge \mathrm{r}\left(e_n\right)=\varpi$

• A path from $\upsilon$ to $\upsilon$ is called a cycle.
• Two vertices, $\upsilon$ and $\varpi$ are said to be connected if there exists a path from $\upsilon$ to $\varpi$.
• The existence of paths from vertices gives an equivalence relation on the set of vertices.
• A subgraph is a component of the graph.
• A graph with only one component is said to be connected.

Definition 7 (Incidence Matrix).

Let $\mathcal{G}=(V,E,\mathrm{l},\mathrm{r})$ be a directed graph that is finite and loop-free such that $E=e_1,\dots ,e_m$ and $V={\upsilon }_1,\dots ,{\upsilon }_n$. Then the $j^{th}$ row and $k^{th}$ columns of the incidence matrix, $A_0\in {\mathbb{R}}^{n\times m}$ of $\mathcal{G}$ is given by

$$\begin{array}{c}\hfill a_{jk}=\left\{\begin{array}{cc}\hfill 1\hfill & \hfill \mathrm{l}\left(e_k\right)={\upsilon }_j,\hfill \\ \hfill -1\hfill & \hfill \mathrm{r}\left(e_k\right)={\upsilon }_j,\hfill \\ \hfill 0\hfill & \hfill \mathrm{otherwise}.\hfill \end{array}\right.\end{array}$$

If the $\mathrm{rank}\left(A_0\right)=n-k$, then the graph $\mathcal{G}$ has $k\in \mathbb{N}$ components, such that k rows can be removed from $A_0$ resulting in matrix with same rank.

Definition 8 (Kirchhoff-Dirac Structure, Kirchhoff-Lagrange Submanifold).

The Kirchhoff-Dirac structure of $\mathcal{G}$ can thus be defined by the set

$$\begin{array}{c}\hfill {\mathcal{D}}_K^S\left(\mathcal{G}\right)\left\{(j,i,\varphi ,u)\in {\mathbb{R}}^{n-\left|\mathsf{S}\right|}\times {\mathbb{R}}^m\times {\mathbb{R}}^{n-\left|S\right|}\times {\mathbb{R}}^m\left|\left[\begin{array}{cc}I& A\\ 0& 0\end{array}\right]\left(\begin{array}{c}j\\ i\end{array}\right)+\left[\begin{array}{cc}I& A\\ 0& 0\end{array}\right]\left(\begin{array}{c}\varphi \\ u\end{array}\right)=0\right.\right\}\end{array}$$

(1)

where i and u are the currents and voltages at the edges of the graph, respectively, whereas q and $\varphi$ are the charges and potentials at the vertices of the graph, respectively.

Assuming that $S=\left\{v_1,\dots ,v_{\left|S\right|}\right\}$, the Kirchhoff-Lagrange submanifold of $\mathcal{G}$ with respect to S is defined as

$$\begin{array}{c}\hfill {\mathcal{L}}_K^S\left(\mathcal{G}\right)0\times {\mathbb{R}}^{n-\left|S\right|}\subset {\mathbb{R}}^{n-\left|S\right|}\times {\mathbb{R}}^{n-\left|S\right|}\end{array}$$

(2)

where $\mathcal{G}=(V,E,\mathrm{l},\mathrm{r})$ is a directed graph that

• Is finite
• Is loop-free
• Has an incidence matrix $A_0\in {\mathbb{R}}^{n\times m}$

If ${\mathcal{G}}_1,\dots ,{\mathcal{G}}_k$ are the components of $\mathcal{G}$ with corresponding vertices $V_1,\dots ,V_k\subset V$ so that there exists a subset $\mathcal{S}\subset V$ such that $\mathcal{S}$ containing at most one vertex from each component, i.e. $\forall s,s^{\prime },\in S,i\le k:v,v^{\prime }\in V_i\Rightarrow v=v^{\prime }$, then $A\in {\mathcal{R}}^{(n-k)\times m}$ can be constructed by deleting the rows corresponding to the vertices $\mathcal{S}$ from $A_0\in {\mathbb{R}}^{n\times m}$.

Remark. According to Proposition 1, Equations 1 and 2 indicate that ${\mathcal{D}}_K^S\left(\mathcal{G}\right)$ is a DS and ${\mathcal{L}}_K^S\left(\mathcal{G}\right)$ is a LS in the space ${\mathbb{R}}^{n-\left|S\right|}\times {\mathbb{R}}^{n-\left|S\right|}$.

Definition 8 caters for an introduction of a pH system $({\mathcal{D}}_K^S\left(\mathcal{G}\right),{\mathcal{L}}_K^S\left(\mathcal{G}\right),\left\{0\right\})$ with dynamics

$$\begin{array}{c}\hfill (-\frac{d}{dt}q\left(t\right),i\left(t\right),\varphi \left(t\right),u\left(t\right))\in {\mathcal{D}}_K^S\left(\mathcal{G}\right)\end{array}$$

(3)

where $(q\left(t\right),\varphi \left(t\right))\in {\mathcal{L}}_K^S\left(\mathcal{G}\right)$.

Using the equivalence of $(q\left(t\right),\varphi \left(t\right))\in {\mathcal{L}}_K^S\left(\mathcal{G}\right)$ to $q\left(t\right)=0$ and $\varphi \left(t\right)\in {\mathcal{R}}^{n-\left|S\right|}$, we see that equation 3 holds, iff, the condition holds.

$$\begin{array}{c}\hfill q\left(t\right)=0\wedge Ai\left(t\right)=0\wedge A^{\mathsf{T}}\varphi 2\left(t\right)=-u\left(t\right).\end{array}$$

$$\begin{array}{c}\hfill b_{jl}=\left\{\begin{array}{cc}\hfill 1\hfill & \hfill e_l\in C_j\mathrm{with}\mathrm{orientations}\mathrm{that}\mathrm{coincide},\hfill \\ \hfill -1\hfill & \hfill e_l\in C_j\mathrm{with}\mathrm{orientations}\mathrm{that}\mathrm{do}\mathrm{not}\mathrm{coincide},\hfill \\ \hfill 0\hfill & \hfill \mathrm{otherwise}\hfill \end{array}\right.\end{array}$$

On this basis, one can define a DS

$$\begin{array}{c}\hfill {\mathcal{D}}_K^{\prime }\left(\mathcal{G}\right)\left\{(j,i,\varphi ,u)\in {\mathbb{R}}^{n-\left|\mathsf{S}\right|}\times {\mathbb{R}}^m\times {\mathbb{R}}^{n-\left|S\right|}\times {\mathbb{R}}^m\left|\left[\begin{array}{cc}I& A\\ 0& 0\end{array}\right]\left(\begin{array}{c}j\\ i\end{array}\right)+\left[\begin{array}{cc}I& A\\ 0& 0\end{array}\right]\left(\begin{array}{c}\varphi \\ u\end{array}\right)=0\right.\right\}\end{array}$$

with dynamics

$$\begin{array}{c}\hfill \left(-\frac{d}{dt}\psi \left(t\right),i\left(t\right),\iota \left(t\right),u\left(t\right)\right)\in {\mathcal{D}}_K^{\prime }\left(\mathcal{G}\right)\end{array}$$

and an LS

$$\begin{array}{c}\hfill {\mathcal{L}}_K^{\prime }\left(\mathcal{G}\right)\left\{0\right\}\times {\mathbb{R}}^{n-m+k}\end{array}$$

with dynamics

$$\begin{array}{c}\hfill \left(\psi \left(t\right),\iota \left(t\right)\right)\in {\mathcal{L}}_K^{\prime }\left(\mathcal{G}\right)\end{array}$$

Together, they form a the pH system defined by the set $(D_K^{\prime }\left(\mathcal{G}\right),{\mathcal{L}}_K^{\prime }\left(\mathcal{G}\right),\left\{0\right\})$.

3.1. Dirac Structure

A key feature of a Dirac structure is the fact that the standard composition of two Dirac structures is again a Dirac structure. The implication of this statement is that any power-conserving interconnection of a Port Hamiltonian system is also a Port Hamiltonian system itself. This constitutes the foundation feature in the Port Hamiltonian approach to modelling, simulation and control of complex physical systems. The intricate Dirac structure is the guide to the algebraic constraints of the interconnected system as well as its Casimir functions [13]. The Casimir are significant in the set point regulation of Port Hamiltonian systems. The framework for the Port Hamiltonian allows for port-based modelling. Port-based modelling means that we are interconnecting many different elements through ports. Dirac structures are the tools used to connect multiple elements. These various elements are energy-storing elements, energy-dissipating elements and external elements which could be supplying energy. A diagram to demonstrate the connection structure is given in Figure 1,13]

4. Detailed Port-Hamiltonian Model of a Mechanical Ventilator

4.1. Description of the System

The overall mechanical ventilator system diagram is provided in Figure 2. The entire system consists of various subsystems. In this research work, the main subsystem that contributes to the flow of air and consequently the air pressure are discussed. Therefore, the following main subsystems will be discussed: the DC Motor subsystem, Turbine pump/blower subsystem, Pump-shaft/impeller subsystem and Solenoid valve subsystem,

4.2. Blower Model

In this section, the dynamical system model equations for the blower model are presented in the port-Hamiltonian framework. The blower model comprises three sub-systems, namely a DC Motor which drives the blower-shaft/impeller, a blower-shaft/impeller that couples the DC motor to the fluid, using the rotational motion of the motor to accelerate the fluid and finally the fluid being driven.

In the Port-Hamiltonian perspective, the state vector ${\mathcal{X}}_b\in {\mathbb{R}}^4$, given by

$$\begin{array}{c}\hfill {\mathcal{X}}_b{\left[\begin{array}{cccc}{p}_m,& {\varphi }_m,& P_b,& Q_b\end{array}\right]}^{\top }\end{array}$$

(4)

where $p_m\in {\mathbb{R}}^1$ and ${\varphi }_m\in {\mathbb{R}}^1$ are the angular momentum and magnetic flux of the DC motor, respectively, while $P_b\in {\mathbb{R}}^1$ and $Q_b\in {\mathbb{R}}^1$ are the pressure and flow rate of the blower.

The total energy of the blower is given by the Hamiltonian

$$\begin{array}{ccc}\hfill {\mathcal{H}}_b\left[{\mathcal{X}}_b\right]& =& \hfill \frac{1}{2}\left(\frac{p_m^2}{I_m}+\frac{{\varphi }_m^2}{L_m}+C_b{P}_b^2+I_b{Q}_b^2\right)\hfill \end{array}$$

(5)

where $I_m\in {\mathbb{R}}_{+}^1$ is the inertia and $L_m\in {\mathbb{R}}_{+}^1$ is the inductance of the DC motor, while $C_b\in {\mathbb{R}}_{+}^1$ and $I_b\in {\mathbb{R}}_{+}^1$ is the hydraulic capacitance and inertance of the air in the blower, respectively.

Thus, the Port-Hamiltonian model of a blower is given by

$$\begin{array}{ccc}\hfill {\dot{\mathcal{X}}}_b& =& \hfill \left({\mathcal{J}}_b-{\mathcal{R}}_b\right)\nabla {\mathcal{H}}_b\left[{\mathcal{X}}_b\right]+{\mathcal{G}}_b{u}_b\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill {\mathcal{Y}}_b& =& \hfill {\mathcal{G}}_b^{*}\nabla {\mathcal{H}}_b\left[{\mathcal{X}}_b\right]\hfill \end{array}$$

(7)

where

$$\begin{array}{ccc}\hfill {\mathcal{J}}_b& & \hfill \left[\begin{array}{cccc}0& K_m& 0& -\left(K_o{p}_m\right)/I_m\\ -K_m& 0& 0& 0\\ 0& 0& 0& 1/\left(C_b{I}_b\right)\\ \left(K_o{p}_m\right)/I_m& 0& -1/\left(C_b{I}_b\right)& 0\end{array}\right]\in {\mathbb{R}}^{4\times 4}\hfill \\ \hfill {\mathcal{R}}_b& & \hfill \mathrm{diag}\left(\left[\begin{array}{cccc}{b}_m& R_m& 0& \left(R_b{p}_m\right)/\left(I_b{I}_m\right)\end{array}\right]\right)\in {\mathbb{R}}^{4\times 4}\hfill \\ \hfill \nabla {\mathcal{H}}_b\left[{\mathcal{X}}_b\right]& & \hfill \left[\begin{array}{cccc}{p}_m/I_m& {\varphi }_m/L_m& C_b{P}_b& I_b{Q}_b\end{array}\right]\in {\mathbb{R}}^{4\times 1}\hfill \\ \hfill {\mathcal{G}}_b& & \hfill \left[\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& -(1/C_b)& 0\\ 0& 0& (1/I_b)\end{array}\right]\in {\mathbb{R}}^{4\times 3}\hfill \\ \hfill u_b& & \hfill \left[\begin{array}{ccc}{V}_m& Q_{in}& P_{in}\end{array}\right]\in {\mathbb{R}}^{3\times 1}\hfill \end{array}$$

where $K_m\in {\mathbb{R}}_{+}^1$ is the motor torque constant, $b_m\in {\mathbb{R}}_{+}^1$ is the viscous damping, $R_m\in {\mathbb{R}}_{+}^1$ is the armature resistance and $K_o\in {\mathbb{R}}_{+}^1$ is the motor angular momentum/pressure coupling constant. The inputs to the system are the DC motor voltage, the input volumetric flow rate and pressure given by $V\in {\mathbb{R}}^1$, $Q_{in}\in {\mathbb{R}}^1$ and $P_{in}\in {\mathbb{R}}^1$. One can easily show that ${\mathcal{J}}_m=-{\mathcal{J}}_m^{\top }$ and ${\mathcal{R}}_m={\mathcal{R}}_m^{\top }⪰0$.

Taking the time derivative of the Hamiltonian

$$\begin{array}{ccc}\hfill {\dot{\mathcal{H}}}_b\left[{\chi }_b\right]& =& \hfill u_b·{\mathcal{Y}}_b-\nabla {\mathcal{H}}_b\left[{\mathcal{X}}_b\right]·({\mathcal{R}}_b\nabla {\mathcal{H}}_b\left[{\mathcal{X}}_b\right])\hfill \end{array}$$

(8)

This system has power and resistive ports.

5. Solenoid Valve Subsystem

Figure 3 shows a solenoid valve in the open and closed position. It is assumed that the air gaps are sufficiently small such that the effect of fringing of the magnetic flux is negligible. Consider a solenoid in which the permeability of the core and the length of the part of the magnetic circuit inside the core are denoted by ${\mu }_c\in {\mathbb{R}}_{+}^1$ and $l_c\in {\mathbb{R}}_{+}^1$, respectively. The equivalent length of the solenoid’s magnetic circuit, $l_{eq}\left[·\right]:\mathbb{R}\to \mathbb{R}$ is dependent on the displacement of the spool, $q_s\in \mathbb{R}$, and can be written as

$$\begin{array}{ccc}\hfill l_{eq}\left[q_s\right]& =& \hfill l_c+\frac{{\mu }_c}{{\mu }_0}(q_{s_{\mathrm{tot}}}-q_s)\hfill \end{array}$$

(9)

where ${\mu }_0\in {\mathbb{R}}_{+}^1$ are the permeability of air and $q_{{\mathrm{s}}_{\mathrm{tot}}}\in {\mathbb{R}}_{+}^1$ is the total air-gap. The solenoid coordinate systems is represented in Figure 4.

Figure 4. Solenoid coordinate system

Figure 4. Solenoid coordinate system

Figure 3. Diagram of a Solenoid valve in the i) open position which allows fluid flow and ii) closed position which stops fluid flow

Figure 3. Diagram of a Solenoid valve in the i) open position which allows fluid flow and ii) closed position which stops fluid flow

Thus, the inductance of the solenoid varies with displacement of the spool and hence can be expressed by the function $L_s\left[·\right]:\mathbb{R}\to \mathbb{R}$ given by

$$\begin{array}{c}\hfill L_s\left[q_s\right]=\frac{N^2{A}_e{\mu }_c}{l_{eq}\left[q_s\right]}\end{array}$$

(10)

where $N\in {\mathbb{R}}_{+}^1$ is the number of turns in the coil of the solenoid and $A_e\in {\mathbb{R}}_{+}^1$ is the effective cross-sectional area of the path of the magnetic flux. The magnetic and mechanical subsystems in the solenoid valve are therefore coupled magnetically due to the dependence of the inductance on the displacement of the spool.

The total energy of the solenoid is given by the Hamiltonian, ${\mathcal{H}}_s\left[·\right]:{\mathbb{R}}^3\to \mathbb{R}$, which is a function of the state vector ${\mathcal{X}}_s={\left[\begin{array}{ccc}{\varphi }_s,& p_s,& q_s\end{array}\right]}^{\top }\in {\mathbb{R}}^3$, expressed as the sum of the magnetic, kinetic and potential energies denoted by ${\mathcal{H}}_{\mathrm{magnetic}}\left[{\varphi }_s,q_s\right]:{\mathbb{R}}^2\to \mathbb{R}$, ${\mathcal{H}}_{\mathrm{kinetic}}\left[p_s\right]:\mathbb{R}\to \mathbb{R}$ and ${\mathcal{H}}_{\mathrm{potential}}\left[q_s\right]:\mathbb{R}\to \mathbb{R}$, respectively. Thus

$$\begin{array}{c}\hfill {\mathcal{H}}_s\left[{\mathcal{X}}_s\right]={\mathcal{H}}_{\mathrm{magnetic}}\left[{\varphi }_s,q_s\right]+{\mathcal{H}}_{\mathrm{kinetic}}\left[p_s\right]+{\mathcal{H}}_{\mathrm{potential}}\left[q_s\right]\end{array}$$

(11)

given a magnetic flux ${\varphi }_s\in \mathbb{R}$.

Assuming that the pretension of the spring is set to $q_0\in \mathbb{R}$, the Hamiltonian is

$$\begin{array}{c}\hfill {\mathcal{H}}_s\left[{\mathcal{X}}_s\right]=\frac{1}{2}\left(\frac{{\varphi }_s^2}{L\left[q_s\right]}+\frac{p_s^2}{m_s}+k_s{(q_s+q_{s_0})}^2\right)+m_s{q}_{s}g\end{array}$$

(12)

where $p_s\in \mathbb{R}$ is the momentum of the spool, $m_s\in {\mathbb{R}}_{+}^1$ is the mass of the spool, $k_s\in {\mathbb{R}}_{+}^1$ is the spring stiffness and $g\in {\mathbb{R}}_{+}^1$ is the acceleration due to gravity.

Hence, the solenoid’s state and output dynamics are expressed in Port-Hamiltonian form in equations 13 and respectively

$$\begin{array}{ccc}\hfill {\dot{\mathcal{X}}}_s& =& \hfill \left({\mathcal{J}}_s-{\mathcal{R}}_s\right)\nabla {\mathcal{H}}_s\left[{\mathcal{X}}_s\right]+{\mathcal{G}}_s{u}_s\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill {\mathcal{Y}}_s& =& \hfill {\mathcal{G}}_s^{*}{u}_s\hfill \end{array}$$

(14)

where

$$\begin{array}{ccc}\hfill {\mathcal{J}}_s& =& \hfill \left[\begin{array}{ccc}0& 0& 0\\ 0& 0& -1\\ 0& 1& 0\end{array}\right],{\mathcal{R}}_s=\left[\begin{array}{ccc}{R}_s& 0& 0\\ 0& b_s& 0\\ 0& 0& 0\end{array}\right],{\mathcal{G}}_s=\left[\begin{array}{cc}1& 0\\ 0& 0& 0\\ 0& (A_{s_1}-A_{s_2})& (A_{s_3}-A_{s_4})\end{array}\right],\hfill \\ \hfill {\mathcal{G}}_s^{*}& =& \hfill {\mathcal{G}}_s^{\top },u_s=\left[\begin{array}{c}{V}_s\\ p_s\end{array}\right]\mathrm{and}\nabla {\mathcal{H}}_s\left[{\mathcal{X}}_s\right]=\left[\begin{array}{c}\frac{{\varphi }_s}{A_e{N}^2{\mu }_c}\left(l_c-\frac{{\mu }_c(q_s-q_{s_{\mathrm{tot}}})}{{\mu }_0}\right)\\ \frac{p_s}{m_s}\\ \frac{A_e{N}^2{\mu }_c^2{\mu }_0}{{\left({\mu }_0{l}_c-{\mu }_c(q_s-q_{s_{\mathrm{tot}}})\right)}^2}+k_s(q_s+q_{s_0})+m_{s}g\end{array}\right]\hfill \end{array}$$

where $R_s\in {\mathbb{R}}_{+}^1$ is the resistance of the coil, $b_s\in {\mathbb{R}}_{+}^1$ is the viscous damping acting on the spool and ∇ is the gradient operator. $A_{s_1},A_{s_2}\in {\mathbb{R}}_{+}^1$ are the various cross-sectional areas of the lands of the spool and the input vector $u_s\in {\mathbb{R}}^2$ consists of the input voltage $V_s\in {\mathbb{R}}^1$ and the supply pressure $p_s\in {\mathbb{R}}^1$. It can be seen that ${\mathcal{J}}_s=-{\mathcal{J}}_s^{\top }\in {\mathbb{R}}^{3\times 3}$ possesses skew-symmetry, while ${\mathcal{R}}_s={\mathcal{R}}_s^{\top }\in {\mathbb{R}}^{3\times 3}$ is positive semi-definite.

Taking the time derivative of the Hamiltonian

$$\begin{array}{ccc}\hfill {\dot{\mathcal{H}}}_s\left[{\chi }_s\right]& =& \hfill u_s·{\mathcal{Y}}_s-\nabla {\mathcal{H}}_s\left[{\mathcal{X}}_s\right]·({\mathcal{R}}_s\nabla {\mathcal{H}}_s\left[{\mathcal{X}}_s\right])\hfill \end{array}$$

(15)

This system has power and resistive ports.

5.1. Pipe Model

In this section, a port-Hamiltonian model of a single pipe segment is developed. The basis of these developments is the Navier Stokes equations for one-dimensional non-stationary flow of gas in a pipe. The following assumptions are taken for the sake of model simplification [14,15]:

• The pipe is taken as rigid (it does not expand in cross-section as a result of fluid flow).
• Frictional and gravitational effects are neglected (this will be relaxed in future works in this research area),
• The model parameters of the gas remain constant along the pipe cross-section but vary in time along the pipe length. Thus they can be averaged about the cross-section and thus the gas flow is one-dimensional.
• The temperatures of the pipe walls are assumed to be constant and equal to the ambient room temperature. Hence temperature effects are ignored.

Taking into account these assumptions, the coordinate system attached to a segment of pipe is illustrated in Figure 5.

Thus, for a given time interval $t_s\le t<t_f$ with start time $t_s$ and finish time $t_f$, the length normalized one-dimensional Euler equations for gas of density $\rho (t,z):(t_s,t_f]\times (0,\ell )\to {\mathbb{R}}^1$ flowing at a velocity $\nu (t,z):(t_s,t_f]\times (0,\ell )\to {\mathbb{R}}^1$ through a pipe of length $\ell \in {\mathbb{R}}_{+}^1$ and cross-sectional area $A_p\in {\mathbb{R}}_{+}^1$ are given by:

$$\begin{array}{ccc}& & \hfill {\partial }_t\left(\rho A_p\right)+{\partial }_z\left(\rho A_p\nu \right)=0\hfill \end{array}$$

(16a)

$$\begin{array}{ccc}& & \hfill {\partial }_t\left(\rho A_p\nu \right)+{\partial }_z(\rho A_p{\nu }^2+p{A}_p)=0\hfill \end{array}$$

(16b)

where $p(t,z):(t_s,t_f]\times (0,\ell )\to {\mathbb{R}}^1$ is the pressure and ${\partial }_i\partial /{\partial }_i$ is the partial derivative with respect to the temporal and spatial variables given by the subscripts $i\in \left\{t,z\right\}$.

5.1.1. Port-Hamiltonian Formulation of Pipe-Flow Model

The fluid dynamics can be written in terms of mass per unit length, i.e. $\varrho \rho A_p$, as well as the fluid momentum, $m\varrho \nu$. Thus equations 16a and can be expressed as

$$\begin{array}{ccc}& & \hfill {\partial }_t\varrho =-{\partial }_{z}m\hfill \end{array}$$

(17a)

$$\begin{array}{ccc}& & \hfill {\partial }_{t}m=-{\partial }_z\left(\frac{m^2}{\varrho }+p{A}_p\right)\hfill \end{array}$$

(17b)

Defining the state vector of the gas flow through a pipe segment as ${\chi }_p\left[\begin{array}{c}\varrho \\ m\end{array}\right]$, the energy of the gas can be expressed in form of a Hamiltonian ${\mathcal{H}}_p\left[{\chi }_p\right]:{\mathbb{R}}^2\to {\mathbb{R}}^1$ given by

$$\begin{array}{c}\hfill {\mathcal{H}}_p\left[{\chi }_p\right]={\int }_0^{\ell }{H}_p\left[{\chi }_p\right]\mathrm{d}x={\int }_0^{\ell }\left(\frac{m^2}{2\varrho }+\varrho U\left[\varrho /A_p\right]\right)\mathrm{d}x\end{array}$$

(18)

where $H_p\left[·\right]:{\mathbb{R}}^2\to {\mathbb{R}}^1$ is the Hamiltonian density and $U\left[·\right]:{\mathbb{R}}^1\to {\mathbb{R}}^1$ is the internal energy of the gas which in the case of an isentropic fluid, can be expressed as a function of density.

The port-Hamiltonian dynamics take the following form

$$\begin{array}{c}\hfill {\partial }_t{\chi }_p=\mathcal{J}\left[{\chi }_p\right]{\delta }_{{\chi }_p}{\mathcal{H}}_p\left[{\chi }_p\right]\end{array}$$

(19)

where $\mathcal{J}\left[{\chi }_p\right]$, the formally skew-symmetric operator and ${\delta }_{{\chi }_p}\mathcal{H}\left[{\chi }_p\right]$, the variational derivative of the Hamiltonian density are expressed as

$$\begin{array}{c}\hfill \mathcal{J}\left[{\chi }_p\right]=-\left[\begin{array}{cc}0& {\partial }_z\\ {\partial }_z& 0\end{array}\right]\mathrm{and}{\delta }_{{\chi }_p}{\mathcal{H}}_p\left[{\chi }_p\right]=\left[\begin{array}{c}h-\frac{m^2}{2{\varrho }^2}\\ \frac{m}{\varrho }\end{array}\right]\end{array}$$

(20)

where h is the enthalpy, and $\mathcal{J}\left[{\chi }_p\right]$ is a formally skew-symmetric operator.

The rate of change of the Hamiltonian can be found as

$$\begin{array}{c}\hfill {\dot{\mathcal{H}}}_p={\mathbf{u}}_p^{\top }{\mathbf{y}}_p\end{array}$$

(21)

where

$$\begin{array}{c}\hfill {\mathbf{u}}_p=W_B{R}_{\mathrm{ext}}\left[\begin{array}{c}{\left.{\delta }_{{\chi }_p}{\mathcal{H}}_p\right|}_0\\ {\left.{\delta }_{{\chi }_p}{\mathcal{H}}_p\right|}_{\ell }\end{array}\right]\mathrm{and}{\mathbf{y}}_p=W_C{R}_{\mathrm{ext}}\left[\begin{array}{c}{\left.{\delta }_{{\chi }_p}{\mathcal{H}}_p\right|}_0\\ {\left.{\delta }_{{\chi }_p}{\mathcal{H}}_p\right|}_{\ell }\end{array}\right]\end{array}$$

(22)

with components given by

$$\begin{array}{c}\hfill W_B=\frac{1}{\sqrt{2}}\left[\begin{array}{cccc}1& 0& 0& 1\\ 0& -1& 1& 0\end{array}\right],W_C=\frac{1}{\sqrt{2}}\left[\begin{array}{cccc}0& 1& 1& 0\\ 1& 0& 0& -1\end{array}\right],R_{\mathrm{ext}}=\frac{1}{\sqrt{2}}\left[\begin{array}{cccc}0& -1& 0& 1\\ -1& 0& 1& 0\\ 1& 0& 1& 0\\ 0& 1& 0& 1\end{array}\right]\end{array}$$

$$\begin{array}{c}\hfill \mathrm{and}\left[\begin{array}{c}{\left.{\delta }_{{\chi }_p}{\mathcal{H}}_p\right|}_0\\ {\left.{\delta }_{{\chi }_p}{\mathcal{H}}_p\right|}_{\ell }\end{array}\right]=\left[\begin{array}{c}{\left.\left(h-\frac{m^2}{2{\varrho }^2}\right)\right|}_0\\ {\left.\left(\frac{m}{\varrho }\right)\right|}_0\\ {\left.\left(h-\frac{m^2}{2{\varrho }^2}\right)\right|}_{\ell }\\ {\left.\left(\frac{m}{\varrho }\right)\right|}_{\ell }\end{array}\right]\end{array}$$

Expanding ${\mathbf{u}}_p$ and ${\mathbf{y}}_p$

$$\begin{array}{c}\hfill {\mathbf{u}}_p=\left[\begin{array}{c}{\left.\left(\frac{m}{\varrho }\right)\right|}_{\ell }\\ {\left.\left(h-\frac{m^2}{2{\varrho }^2}\right)\right|}_0\end{array}\right]\mathrm{and}{\mathbf{y}}_p=\left[\begin{array}{c}{\left.\left(h-\frac{m^2}{2{\varrho }^2}\right)\right|}_{\ell }\\ -{\left.\left(\frac{m}{\varrho }\right)\right|}_0\end{array}\right]\end{array}$$

(23)

Thus, the rate of change of the Hamiltonian is

$$\begin{array}{c}\hfill {\dot{\mathcal{H}}}_p={\left.\left(h-\frac{m^2}{2{\varrho }^2}\right)\right|}_{\ell }{\left.\left(\frac{m}{\varrho }\right)\right|}_{\ell }-{\left.\left(\frac{m}{\varrho }\right)\right|}_0{\left.\left(h-\frac{m^2}{2{\varrho }^2}\right)\right|}_0\end{array}$$

(24)

5.2. Electric Circuit Model of the Lung

The port Hamiltonian formulation for nonlinear electric circuits is presented in this section. Since the main focus of this article is on the mechanical ventilator model, the lung model is simplified by considering an electric circuit analogy. The model under consideration is that of a fully sedated patient who relies completely on the mechanical ventilator to breathe. The circuit model is shown in Figure 6.

The circuit can be represented in the form of a network graph using graph theory. The graph of the circuit given in Figure 6 is given in Figure 7. The model has $n=\left|\mathcal{V}\right|=7$ vertices.

The complete graph of the circuit, $A_0$, is

The ventilator input voltage is represented as a source $\mathcal{S}=({\mathcal{D}}_{\mathcal{S}},{\mathcal{L}}_{\mathcal{S}},{\mathcal{R}}_{\mathcal{S}})$. The set of sources is $\mathcal{S}=\left\{v_7\right\}$ and the dimension $\left|\mathcal{S}\right|=1$ source. The circuit also consists of $m_{\mathcal{C}}=5$ capacitors (storage elements) that are represented as ${\mathcal{C}}_i=({\mathcal{D}}_{{\mathcal{C}}_i},{\mathcal{L}}_{{\mathcal{C}}_i},{\mathcal{R}}_{{\mathcal{C}}_i})$ where the $i^{\mathrm{th}}$ index is used to distinguish between the components is $i=\left\{\begin{array}{ccccc}{C}_{\mathrm{L}},& C_{\mathrm{T}},& C_{\mathrm{B}},& C_{\mathrm{A}},& C_{\mathrm{CW}}\end{array}\right\}$. Furthermore there are $m_{\mathcal{R}}=4$ resistors (dissipative elements) that are represented as ${\mathcal{R}}_i=({\mathcal{D}}_{{\mathcal{R}}_i},{\mathcal{L}}_{{\mathcal{R}}_i},{\mathcal{R}}_{{\mathcal{R}}_i})$ where the $i^{\mathrm{th}}$ index is given by $i=\left\{\begin{array}{cccc}{R}_{\mathrm{ML}},& R_{\mathrm{IT}},& R_{\mathrm{TB}},& R_{\mathrm{BA}}\end{array}\right\}$.

Selecting node 7 as the ground, the reduced incidence matrix, $\mathcal{A}$, is given by

where the dimensions of the components of $\mathcal{A}$ are $dim\left(A_{\mathcal{R}}\right)=(n-|S\left|\right)\times m_{\mathcal{R}}=6\times 4$, $dim\left(A_{\mathcal{C}}\right)=(n-|S\left|\right)\times m_{\mathcal{C}}=6\times 5$ and $dim\left(A_{\mathcal{S}}\right)=(n-|S\left|\right)\times m_{\mathcal{S}}=6\times 1$.

The total energy of the circuit is given by the Hamiltonian

$$\begin{array}{c}\hfill {\mathcal{H}}_C=\sum _{\forall i}{\mathcal{H}}_{C_i}\left(q_{C_i}\right)=\frac{1}{2}\sum _{\forall i}\frac{q_{C_i}}{C_i}\end{array}$$

(27)

where $q_{C_i}\in {\mathbb{R}}^1$ is the charge of the $i^{\mathrm{th}}$ capacitor for $i=\left\{\begin{array}{ccccc}{C}_{\mathrm{L}},& C_{\mathrm{T}},& C_{\mathrm{B}},& C_{\mathrm{A}},& C_{\mathrm{CW}}\end{array}\right\}$

$$\begin{array}{c}\hfill \left[\begin{array}{cccc}\mathbb{I}& A_{\mathcal{C}}& A_{\mathcal{R}}& A_{\mathcal{S}}\\ \mathbb{O}& \mathbb{O}& \mathbb{O}& \mathbb{O}\\ \mathbb{O}& \mathbb{O}& \mathbb{O}& \mathbb{O}\\ \mathbb{O}& \mathbb{O}& \mathbb{O}& \mathbb{O}\end{array}\right]\left(\begin{array}{c}-\frac{d}{dt}q\\ -\frac{d}{dt}{q}_{\mathcal{C}}\\ -i_{\mathcal{R}}\\ -i_{\mathcal{S}}\end{array}\right)+\left[\begin{array}{cccc}\mathbb{O}& \mathbb{O}& \mathbb{O}& \mathbb{O}\\ -A_{\mathcal{R}}^{\top }& \mathbb{O}& \mathbb{I}& \mathbb{O}\\ -A_{\mathcal{C}}^{\top }& \mathbb{I}& \mathbb{O}& \mathbb{O}\\ -A_{\mathcal{S}}^{\top }& \mathbb{O}& \mathbb{O}& \mathbb{I}\end{array}\right]\left(\begin{array}{c}\varphi \\ u_{\mathcal{C}}\\ u_{\mathcal{R}}\\ u_{\mathcal{S}}\end{array}\right)=\mathbb{O}\end{array}$$

(28)

where $\mathbb{I}$ and $\mathbb{O}$ are appropriately sized identity and zero matrices, respectively and $q=0$ and $\varphi$ is the node potential. Expanding Equation 41

$$\begin{array}{c}\hfill A_{\mathcal{C}}\frac{d}{dt}{q}_{\mathcal{C}}+A_{\mathcal{R}}\left(A_{\mathcal{R}}^{\top }\varphi \right)+A_{\mathcal{S}}{i}_{\mathcal{S}}=0\end{array}$$

(29)

$$\begin{array}{c}\hfill -A_{\mathcal{C}}^{\top }\varphi +u_{\mathcal{C}}=0\end{array}$$

(30)

$$\begin{array}{c}\hfill -A_{\mathcal{S}}^{\top }\varphi +u_{\mathcal{S}}=0\end{array}$$

(31)

The Dirac structure is given by the set where $\left(\mathcal{D},\mathcal{L},\mathcal{R}\right)$

$$\begin{array}{c}\hfill \left(\begin{array}{cccccccc}-\frac{d}{dt}q,& -\frac{d}{dt}{q}_{\mathcal{C}},& -i_{\mathcal{R}},& -i_{\mathcal{S}},& \varphi ,& u_{\mathcal{C}},& u_{\mathcal{R}},& u_{\mathcal{S}}\end{array}\right)\in \mathcal{D}\end{array}$$

(32a)

$$\begin{array}{c}\hfill \left(\begin{array}{cccc}q,& q_{\mathcal{C}},& \varphi ,& u_{\mathcal{C}}\end{array}\right)\in \mathcal{L}\end{array}$$

(32b)

$$\begin{array}{c}\hfill \left(\begin{array}{cc}-i_{\mathcal{R}},& e_{\mathcal{R}}\end{array}\right)\in \mathcal{R}\end{array}$$

(32c)

6. Model Network Topology

Definition 1.(Directed graph) [17] A directed graph, $\mathcal{G}$ denotes the pair $(\mathcal{V},\mathcal{A})$, where $\mathcal{V}\left(\mathcal{D}\right)$ and $\mathcal{A}\left(\mathcal{D}\right)$ denote the set of vertices and arcs respectively. An arc is a distinct ordered pair of vertices.

The directed graph of a mechanical ventilator model is composed of the disjoint union of the vertices associated with the patient, as well as the inspiratory and expiratory limbs of the ventilator given by $v_O,{\mathcal{V}}_I$ and ${\mathcal{V}}_E$ respectively. Similarly, the patient, inspiratory and expiratory arcs given by $a_O,{\mathcal{A}}_I$ and ${\mathcal{A}}_E$, respectively. Each arc and vertex can further be subdivided into those associated with blowers, valves, pipes and electric circuits present in the mechanical ventilator. Thus the totality of vertices and arcs are given by

$$\begin{array}{c}\hfill \mathcal{V}\left(\mathcal{D}\right)=\left\{v_O\right\}\cup {\mathcal{V}}_I\cup {\mathcal{V}}_E=\left\{v_O\right\}\cup ({\mathcal{V}}_{I_b}\cup {\mathcal{V}}_{I_v}\cup {\mathcal{V}}_{I_p}\cup {\mathcal{V}}_{I_e})\cup ({\mathcal{V}}_{E_b}\cup {\mathcal{V}}_{E_v}\cup {\mathcal{V}}_{E_p}\cup {\mathcal{V}}_{E_e})\end{array}$$

and

$$\begin{array}{c}\hfill \mathcal{A}\left(\mathcal{D}\right)=\left\{a_O\right\}\cup {\mathcal{A}}_I\cup {\mathcal{A}}_E=\left\{a_O\right\}\cup ({\mathcal{A}}_{I_b}\cup {\mathcal{A}}_{I_v}\cup {\mathcal{A}}_{I_p}\cup {\mathcal{A}}_{I_e})\cup ({\mathcal{A}}_{E_b}\cup {\mathcal{A}}_{E_v}\cup {\mathcal{A}}_{E_p}\cup {\mathcal{A}}_{E_e})\end{array}$$

respectively, where the subscripts $i\in \left\{O,p,b,v,p,e\right\}$ are used to indicate the patient, blowers, valves, pipes and electric circuits components, respectively.

Figure 9 shows a detailed graph of the mechanical ventilator given in Figure 2 indicating the patient vertex ${\mathcal{V}}_O$, the arcs and vertices along the inspiratory path given by

$$\begin{array}{c}\hfill {\mathcal{A}}_I\left\{{\mathcal{A}}_{I_{V12}},{\mathcal{A}}_{I_{V23}},{\mathcal{A}}_{I_{V3B}},{\mathcal{A}}_{I_{V23}},{\mathcal{A}}_{I_{BV4}},{\mathcal{A}}_{I_{V4S1}},{\mathcal{A}}_{I_{BV5}},{\mathcal{A}}_{I_{V5S2}}\right\}\end{array}$$

and

$$\begin{array}{c}\hfill {\mathcal{V}}_I\left\{{\mathcal{V}}_{I_{V1}},{\mathcal{V}}_{I_{V2}},{\mathcal{V}}_{I_{V3}},{\mathcal{V}}_{I_B},{\mathcal{V}}_{I_{V4}},{\mathcal{V}}_{I_{S1}},{\mathcal{V}}_{I_{V5}},{\mathcal{V}}_{I_{S2}}\right\}\end{array}$$

respectively, as well as the arcs and vertices along the expiratory path given by

$$\begin{array}{c}\hfill {\mathcal{A}}_E\left\{{\mathcal{A}}_{O{E}_{V1}},{\mathcal{A}}_{E_{V12}},{\mathcal{A}}_{E_{V23}}\right\}\mathrm{and}{\mathcal{V}}_E\left\{{\mathcal{V}}_{E_{V1}},{\mathcal{V}}_{E_{V2}},{\mathcal{V}}_{E_{V3}}\right\}\end{array}$$

respectively.

Figure 8. Simplified graph of the mechanical ventilator given in Figure 2 indicating the elements belonging to the patient, given by the arcs ${\mathcal{A}}_{O_{21}}\in {\mathcal{A}}_O$ and vertices $\left\{{\mathcal{V}}_{O_1},{\mathcal{V}}_{O_2}\right\}\in {\mathcal{V}}_O$ as well as elements of the inspiratory and expiratory arcs and vertices given by and ${\mathcal{V}}_I,{\mathcal{V}}_E$ and ${\mathcal{A}}_{IO},{\mathcal{A}}_{OE}$, respectively.

Figure 8. Simplified graph of the mechanical ventilator given in Figure 2 indicating the elements belonging to the patient, given by the arcs ${\mathcal{A}}_{O_{21}}\in {\mathcal{A}}_O$ and vertices $\left\{{\mathcal{V}}_{O_1},{\mathcal{V}}_{O_2}\right\}\in {\mathcal{V}}_O$ as well as elements of the inspiratory and expiratory arcs and vertices given by and ${\mathcal{V}}_I,{\mathcal{V}}_E$ and ${\mathcal{A}}_{IO},{\mathcal{A}}_{OE}$, respectively.

Figure 9. Detailed graph of the Mechanical Ventilator given in Figure 8 indicating the patient vertex ${\mathcal{V}}_O$, as well as the arcs and vertices along the inspiratory and expiratory paths.

Figure 9. Detailed graph of the Mechanical Ventilator given in Figure 8 indicating the patient vertex ${\mathcal{V}}_O$, as well as the arcs and vertices along the inspiratory and expiratory paths.

7. Model Interconnection/Coupling Conditions

In this section, the coupling conditions of the port-Hamiltonian network model of the mechanical ventilator are given. The types of interconnections occurring in this model are given below:

• Pump-to-pipe interconnection The pressure and flow rate of the fluid exiting the pump, $P_b$ and $Q_b$ respectively, are equal to the pressure and flow rate at the inlet of the pipe given by $p\left(0\right)$ and $m\left(0\right)$, respectively. Thus

  $$\begin{array}{c}\hfill P_b=p_0\mathrm{and}{Q}_b=m\left(0\right)\end{array}$$
• Pipe to valve interconnection The pressure and flow rate of the fluid entering/exiting a valve, $P_v$ and $Q_v$ respectively, are equal to the pressure and flow rate at the inlet/outlet of the pipe given by $p\left(0\right)$ and $m\left(0\right)$ for the inlet and $p(\ell )$ and $m(\ell )$. Thus

  $$\begin{array}{ccc}& & \hfill P_v=p\left(0\right)\mathrm{and}{Q}_v=m\left(0\right)\mathrm{at}\mathrm{the}\mathrm{inlet}\hfill \\ & & \hfill P_v=p(\ell )\mathrm{and}{Q}_v=m(\ell )\mathrm{at}\mathrm{the}\mathrm{outlet}\hfill \end{array}$$
• Pipe to circuit interconnection The pressure and fluid flow-rate at the outlet of a pipe can act as inputs to a circuit model, thus

  $$\begin{array}{c}\hfill p(\ell )=v_{c_{\mathrm{in}}}\mathrm{and}m(\ell )=i_{c_{\mathrm{in}}}\end{array}$$
  On the other hand, the output voltage and current of a circuit can be interconnected to a fluid pipe at the inlet of the pipe. In this case, the output voltage and or current of the circuit should be equal to the inlet pressure and inlet flow rate respectively. This relation can expressed mathematically as:

  $$\begin{array}{c}\hfill v_{c_{out}}=p\left(0\right)\mathrm{and}{i}_{c_{out}}=m\left(0\right)\mathrm{at}\mathrm{the}\mathrm{inlet}\end{array}$$
  The Hamiltonian of the complete system is given by the sum of the Hamiltonian’s of the individual systems

  $$\begin{array}{c}\hfill \mathcal{H}={\mathcal{H}}_b+{\mathcal{H}}_s+{\mathcal{H}}_p+{\mathcal{H}}_{\mathcal{C}}\end{array}$$

  (33)
  The rate of change of energy of the complete system is

  $$\begin{array}{ccc}\hfill \dot{\mathcal{H}}=& & \hfill {\dot{P}}_b({\partial }_{z}p\left(0\right)+F_p)+{\dot{Q}}_b({\partial }_{z}m\left(0\right)+M_p)\hfill \\ & & \hfill +{\dot{P}}_{v_{in}}({\partial }_{z}p\left(0\right)+F_{v_{in}})+{\dot{Q}}_{v_{in}}({\partial }_{z}m\left(0\right)+M_{v_{in}})\hfill \\ & & \hfill +{\dot{P}}_{v_{out}}({\partial }_{z}p(\ell )+F_{v_{out}})+{\dot{Q}}_{v_{out}}({\partial }_{z}m(\ell )+M_{v_{out}})\hfill \\ & & \hfill +{\dot{P}}_{v_{out}}({\partial }_{z}p(\ell )+F_{v_{out}})+{\dot{Q}}_{v_{out}}({\partial }_{z}m(\ell )+M_{v_{out}})+\hfill \\ & & \hfill {\dot{\mathcal{H}}}_b+{\dot{\mathcal{H}}}_s+{\dot{\mathcal{H}}}_c\hfill \end{array}$$

  (34)
  The terms $F_p,M_p,F_{v_{in}},M_{v_{in}},F_{v_{out}}$ and $M_{v_{out}}$ are the external pressure and flow rates acting on the system. They should be equal to zero to complete the interconnection.

8. Structure Preserving Discretization

The port Hamiltonian model of the pipe is a partial differential equation, continuous in space. As such it is difficult to simulate the dynamics of a pipe section. In order to do this, it is necessary to approximate the model with a discrete model, in this case, a finite difference model which is an approximation of the original system. Within the context of port Hamiltonian systems, an additional requirement is the need to ensure that the discrete approximation maintains the structural properties of the original system e.g. skew symmetry etc.

Each system state can be replaced with a discrete approximation consisting of a total on n elements as can be seen in Figure 10. As such, $\varrho \approx {\varrho }_d=\left[\begin{array}{cccc}{\varrho }_1,& {\varrho }_2,& \cdots ,& {\varrho }_n\end{array}\right]\in {\mathbb{R}}^{n\times 1}$ and $m\approx m_d=\left[\begin{array}{cccc}{m}_1,& m_2,& \cdots ,& m_n\end{array}\right]\in {\mathbb{R}}^{n\times 1}$. As such, state vector ${\chi }_p$ can be replaced by a discrete approximation ${\chi }_d=\left[\begin{array}{cccccc}{\varrho }_1,& \cdots ,& {\varrho }_n,& m_1,& \cdots ,& m_n\end{array}\right]\in {\mathbb{R}}^{2n\times 1}$. The $i^{\mathrm{th}}$ element of ${\varrho }_i,m_i\in {\mathbb{R}}^{1\times 1}$ is located at $z=\left\{\Delta (i-1),\Delta (i-1/2)\right\}$, where $\Delta$ is the fixed discrete step size between points and $i=1,2,\cdots ,n$. In addition, the efforts at the boundaries are given by ${\delta }_{{\chi }_0}{\mathcal{H}}_0$ and ${\delta }_{{\chi }_n}{\mathcal{H}}_n$. Thus the Hamiltonian given in Equation 18 can be approximated by a discrete approximation such that ${\mathcal{H}}_p\left[{\chi }_p\right]\approx \Delta {\mathcal{H}}_d\left[{\chi }_d\right]$ so that the discrete system effort is now ${\delta }_{{\chi }_d}\left({\mathcal{H}}_d\right)$. A finite difference approximation of the spatial derivatives at the $i^{\mathrm{th}}$ point is

$$\begin{array}{c}\hfill {\left.\frac{\partial }{\partial z}\varrho (t,z)\right|}_i\approx \frac{1}{\Delta }(\varrho (t,z_{i+0.5})-\varrho (t,z_{i-0.5}))\mathrm{and}{\left.\frac{\partial }{\partial z}m(t,z)\right|}_i\approx \frac{1}{\Delta }(m(t,z_{i+1})-m(t,z_i))\end{array}$$

(35)

The central difference approximation at the $i^{\mathrm{th}}$ point is

$$\begin{array}{c}\hfill {\partial }_t\left[\begin{array}{c}{\varrho }_i\\ m_i\end{array}\right]=\frac{1}{\Delta }\left(\left[\begin{array}{c}{\delta }_m\mathcal{H}\left[{\chi }_i\right]\\ {\delta }_{\varrho }\mathcal{H}\left[{\chi }_{i+1}\right]\end{array}\right]-\left[\begin{array}{c}{\delta }_m\mathcal{H}\left[{\chi }_{i-1}\right]\\ {\delta }_{\varrho }\mathcal{H}\left[{\chi }_i\right]\end{array}\right]\right)\end{array}$$

(36)

In matrix form this is

$$\begin{array}{ccc}\hfill {\partial }_t{\varrho }_d& =& \hfill \frac{1}{\Delta }\left[\begin{array}{cccc}-1& & & \\ 1& -1& & \\ & \ddots & \ddots & \\ & & 1& -1\end{array}\right]{\delta }_m\mathcal{H}\left[{\chi }_d\right]+\frac{1}{\Delta }\left[\begin{array}{c}1\\ 0\\ ⋮\\ 0\end{array}\right]{\delta }_m\mathcal{H}\left[{\chi }_d\right]\hfill \end{array}$$

(37)

$$\begin{array}{ccc}\hfill {\partial }_t{m}_d& =& \hfill \frac{1}{\Delta }\left[\begin{array}{cccc}1& -1& & \\ & \ddots & \ddots & \\ & & 1& -1\\ & & & 1\end{array}\right]{\delta }_{\varrho }\mathcal{H}\left[{\chi }_d\right]+\frac{1}{\Delta }\left[\begin{array}{c}0\\ ⋮\\ 0\\ -1\end{array}\right]{\delta }_{\varrho }\mathcal{H}\left[{\chi }_0\right]\hfill \end{array}$$

(38)

which can be re-written as

$$\begin{array}{ccc}\hfill {\partial }_t{\varrho }_d& =& \hfill D{\delta }_m\mathcal{H}\left[{\chi }_d\right]+{\mathcal{G}}_{\varrho }{\delta }_m\mathcal{H}\left[{\chi }_0\right]\hfill \end{array}$$

(39a)

$$\begin{array}{ccc}\hfill {\partial }_t{m}_d& =& \hfill -D^{\top }{\delta }_{\varrho }\mathcal{H}\left[{\chi }_d\right]+{\mathcal{G}}_m{\delta }_{\varrho }\mathcal{H}\left[{\chi }_n\right]\hfill \end{array}$$

(39b)

where

$$\begin{array}{c}\hfill D=\frac{1}{\Delta }\left[\begin{array}{cccc}-1& & & \\ 1& -1& & \\ & \ddots & \ddots & \\ & & 1& -1\end{array}\right],{\mathcal{G}}_{\varrho }=\frac{1}{\Delta }\left[\begin{array}{c}1\\ 0\\ ⋮\\ 0\end{array}\right]\mathrm{and}{\mathcal{G}}_m=\frac{1}{\Delta }\left[\begin{array}{c}0\\ ⋮\\ 0\\ -1\end{array}\right]\end{array}$$

The skew-symmetric operator is clear from Equation Section 8.

$$\begin{array}{ccc}\hfill \left[\begin{array}{c}{\partial }_t{\varrho }_d\\ {\partial }_t{m}_d\end{array}\right]& =& \hfill \left[\begin{array}{cc}& D\\ -D^{\top }\end{array}\right]\left[\begin{array}{c}{\delta }_{\varrho }\mathcal{H}\left[{\chi }_d\right]\\ {\delta }_m\mathcal{H}\left[{\chi }_d\right]\end{array}\right]+\left[\begin{array}{cc}& {\mathcal{G}}_m\\ {\mathcal{G}}_{\varrho }\end{array}\right]\left[\begin{array}{c}{\delta }_{\varrho }\mathcal{H}\left[{\chi }_0\right]\\ {\delta }_m\mathcal{H}\left[{\chi }_n\right]\end{array}\right]\hfill \end{array}$$

(40)

Figure 10. Staggered grid discretization of the one dimensional port Hamiltonian pipe dynamic model

Figure 10. Staggered grid discretization of the one dimensional port Hamiltonian pipe dynamic model

9. Results and Discussion

The pipe and solenoid valve model parameters used in this work are given in Table 1 and Table 2, respectively.

9.1. Model Validation

The model was validated using parameters obtained from the literature. The simulated results provided in the results and discussion section were found to be comparable to those in the existing literature.

9.2. Simulation Environment

Simulations were conducted using MATLAB. Three conditions were used to simulate the mechanical ventilator behaviour with the following lung conditions: compliance 50 mL/cmH2O and resistance 5cmH2O-s/L, compliance 20 mL/cmH2O and resistance 20cmH2O-s/L and compliance 10 mL/cmH2O and resistance 50cmH2O-s/L. For each condition, all the sizes were tested for 2 min per size and the ventilation curves air volume, The conditions were chosen in order to have a fair comparison with results that exist in literature that used similar lung compliance. pressure and flow over time were obtained. The simulation results are shown in Figure 6.10 and 6.11, where A, B, C, D, and E represent various cams sizes of the mechanical ventilator from extra small (A), small(B), Medium(C), Large (D) and extra-large(E)

Figure 11. Air pressure, volume and flow versus time graph

Figure 11. Air pressure, volume and flow versus time graph

Figure 12. Mechanical Ventilator Simulation results

Figure 12. Mechanical Ventilator Simulation results

The figure show ventilation curves for a lung under two conditions, for a partially damaged lung whose compliance and resistance of 20 mL/cmH2O and 20cmH2O-s/L respectively. The lung with complete damage is represented by lung compliance and resistance of 10 mL/cmH2O and resistance 50cmH2O-s/L. As the cam size increases the ventilation curves maintain the same sinusoidal behavior. For a damaged lung case study, minimum ventilation values of flow over time and air volume are achieved. In a damaged lung, the compliance and resistance of the lung increases the pressure.

10. Conclusions and Recommendations

In these research work, the formulation of a detailed mechanical ventilator in the port Hamiltonian framework. This is followed by a Port Hamiltonian model of the respiratory system and thereafter, these two systems integrated. The work conducted demonstrates the the Port Hamiltonian approach is a valid method in the modelling of integrated Mechanical ventilator human respiratory system.