# Port-Hamiltonian formulation and stabilizing controller for a liquid propelled rocket engine

Jules Gibart, Hélène Piet-Lahanier, Francois Farago

► To cite this version:

Jules Gibart, Hélène Piet-Lahanier, Francois Farago. Port-Hamiltonian formulation and stabilizing controller for a liquid propelled rocket engine. Control Engineering Practice, 2025, 163, pp.106389. 10.1016/j.conengprac.2025.106389. hal-05095585

HAL Id: hal-05095585

https://hal.science/hal-05095585v1

Submitted on 3 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Port-Hamiltonian formulation and stabilizing controller for a liquid propelled rocket engine

Jules Gibartª, Hélène Piet-LahanierªID,*, François Faragoᵇ

## ARTICLE INFO

**Keywords:**

Liquid propellant rocket engine

Stability analysis

Port-Hamiltonian framework

Non-linear control design

## ABSTRACT

Reusable technology in the field of space launchers requires complex maneuvers to land a launcher first stage, requiring variable thrust from the rocket engine. The developments in electrical actuators allowed the introduction of closed-loop controllers for liquid propelled rocket engines (LPRE). While closed-loop controllers have been suggested in the literature with robustness guarantees, few stability proofs have been given. The LPRE is a complex non-linear system, rendering a direct approach to determine a Lyapunov function complex. In this paper, a reformulation of the state-space equations into a model more adapted to stability analysis is proposed, and a passivity approach is derived to prove the stability. In addition, a closed-loop controller that enforces the passivity of the system is designed, with a new equilibrium assignment. Simulated results illustrate the performances of the closed-loop controlled engine.

## 1. Introduction

In the wake of reusable technology, launchers are now requested to operate more complex trajectories with enhanced maneuverability. Initial controllers for rocket engines for non-reusable launchers, were first designed in open-loop and sought for a single operating point to maintain a constant thrust for the duration of the ascend. The development of reusable engines, the space shuttle engine being one of the earliest engine of this type, has made it necessary to develop controls for an extended throttling range. Liquid propelled rocket engines (LPRE) prove to be a suitable solution to perform the variations needed in the power provided and the development of electrical actuators enables thrust control of the engine. These requirements on the functioning points make it necessary to introduce closed loop control and extend the functioning domain of the engine. The design of control laws for such complex non-linear system induces guarantees, around the different operating points and during the state transitions between such operating points (Casiano, Hulka, & Yang, 2010). In Wu, Zhang, Hui, and Zhou (2018), $H_\infty$ techniques have been suggested while in Dresia, Waxenegger-Wilfing, Dos Santos Hahn, Deeken, and Oschwald (2021) the use of reinforcement learning was introduced to provide an efficient control accounting for the variability of the environment. Model Predictive control design was suggested in Pérez-Roca, Marzat, Flayac, Piet-Lahanier, Langlois, Farago, Galeotta, and Gonidec (2019). A fuzzy PID control is defined in Cha and de Oliveira (2021) to follow trajectories for the terminal landing phase. In Gibart,

Piet-Lahanier, Farago, and Galeotta (2023), a contractive controller is proposed for a simplified rocket engine model. These approaches face representativity issues with the high number of state equations and lack emphasizing the problematic interactions that occur between the different elements. In this case proof of the system's stability is no simple task, e.g. by providing a Lyapunov function which has proven difficult (Perez Roca, 2020). However, to the best of our knowledge, there exists a very limited number of fully nonlinear or hybrid frameworks allowing control for a wider throttling range. Most of the control designed are thus limited to low correction margins. This type of control improvements, with more sophisticated monitoring could be a real asset for the lifetime extension of an LPRE, which is a key aspect for their reusability. Port-Hamiltonian systems theory (Duindam, Macchelli, Stramigioli, & Bruyninckx, 2009; Lohmiller, 1999) offers an alternative to classical state-space modeling, which is more adapted to represent complex non-linear systems involving several physical domains. Port-Hamiltonian systems additionally make use of physical invariants, in particular the total energy of the system to derive a storage function that can be used to prove the stability of a system. Such function is in practice easier to obtain than a classical Lyapunov function. The particular geometry of the Port-Hamiltonian framework is particularly relevant for passivity and stability study of a system. A combustion system being easily represented by energy exchanges, the passivity analysis is a well-fitted approach to the LPRE stability problematic.

<footnote>ª DTIS, ONERA, University Paris-Saclay, F-91123 Palaiseau, France</footnote>

<footnote>ᵇ CNES, Paris, 75012, France</footnote>

<footnote>* Corresponding author.
E-mail addresses: jules.gibart@orange.fr (J. Gibart), helene.piet-lahanier@onera.fr (H. Piet-Lahanier), francois.farago@cnes.fr (F. Farago).</footnote>
A rocket engine is composed of several elements each representative of different physical domains (chemical, mechanical, hydrodynamical, ...). Each subsystem composing the engine can be modeled independently under the form of a Port-Hamiltonian structure using its physical features as presented in Rashad, Califano, van der Schaft, and Stramigoli (2020). In Eberard and Maschke (2004), Eberard, Maschke, and Van Der Schaft (2007), Port-Hamiltonian modeling of irreversible thermodynamics has been discussed. Indeed in Esquivel-Sancho, Muñoz-Arias, Phillips-Brenes, and Pereira-Arroyo (2022), a reversible hydropump model is developed. Finally, a model of chemical reactions in a reactor has been developed in García-Sandoval, González-Álvarez, and Calderón (2015).

In this paper, a new formulation of an LPRE's dynamics in Port-Hamiltonian framework is proposed. This new expression allows for a better comprehension of the different interactions that are of major importance for stability analysis. This modeling is used to derive a controller to perform trajectory tracking. It also provides an expression for a Lyapunov function of the system that can be used to prove the stability of the controlled system. A step-by-step design procedure is described for the controller. It allows the LPRE to follow the desired trajectory, while conserving stability properties established during the design. An introductory example is presented on a simplified model, and results are presented on a complex model.

The paper is organized as follows. The description of an LPRE is given in Section 2. A brief presentation of Port Hamiltonian systems is provided in Section 3, and the development of the Port Hamiltonian formulation of the LPRE system is depicted in Section 4. Section 5 explicits the design of a passivity based control law. Numerical results illustrating the performances of the controlled system are presented in Section 6. Concluding remarks close the paper.

## 2. Presentation of a liquid propelled rocket engine

An LPRE system is a powerful combustion system that provides a thrust to the rocket body it is attached to. The thrust intensity is related to both the chamber pressure and the velocity of the mass blown by the combustion out of the engine. The main components are represented on Fig. 1. A main combustion chamber (CC) provides the closed environment for the species to be consumed. The species chosen for this model are liquid hydrogen LH and liquid oxygen LOX. As the injection requires high pressured fluids to further elevate the chamber pressure, pumps are placed upstream for both species. The pumps are powered via a secondary chamber called a gas generator (GG), that redirects exhaust gases into two independent turbines, thus driving the pumps. The secondary chamber is fueled by redirection of a portion of the species into it. Finally, several valves are placed along the pipes to regulate the different mass flows:

* Two valves $S_{VCH}, S_{VCO}$ are located on the main combustion chamber lines. Both of them are considered opened at all times during the permanent regime of the LPRE. However, in the case of a transient behavior (see Perez Roca 2020), the valves have to be controlled accordingly. This paper will only consider the permanent regime of the LPRE, when the valves present fixed surfaces $A_{VCH}$ and $A_{VCO}$.

* Two valves $S_{VGH}, S_{VGO}$ are located on the gas generator lines, and will be used as actuators in the control section of the paper. Both valves regulate the amount of propellants in the gas generator and thus define the pump's power, directly affecting the pressure in the main chamber.

* A last valve $S_{HG}$ is placed to distribute the exhaust gases of the gas generator, regulating the fraction of powered allowed for both propellant lines. During the permanent regime $S_{HG}$ is set to a constant value. This valve affects the mixture ratios of species in both chambers, that has a crucial impact on the combustion's efficiency.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Representation of a gas generator cycle of an LPRE.</figcaption></figure>

Contrary to the transient behavior (Pérez-Roca, Langlois, Marzat, Piet-Lahanier, Galeotta, Farago, & Le Gonidec, 2018), the permanent regime involves a homogeneous time-scale for the variables' evolution. In this paper, hypothesis are made on the behavior of the LPRE, to simplify the expressions of the different state evolutions.

* The temperature of the chambers is fixed to a constant value. The temperature of the combustion chambers mainly depend on the mixture ratio of the species. If well-regulated, the mixture ratio hardly deviates from its optimal value, therefore, the changes in temperature can be considered negligible in front of the other variables.

* The density $\rho$ of the species is constant in the feeding lines. Apart from the hot gases, the species in the fueling lines can be considered liquid thus displaying a constant density.

* The thermic phenomena are neglected in this approach. The main thermic exchange in the LPRE takes place between the combustion chamber interface and a cooling circuit that prevents the chamber's edges from overheating. However the phenomenon has little impact on the chamber's pressure evolution, and during the permanent regime the effect on the mass flow in the cooling circuit can be represented using a static temperature increase of the fluid.

* Finally the turbo-pumps' dynamics, i.e. the dynamic for the resistive torque and the motor torque of the turbo-pump are approximated by second-order polynomial functions using manufacturer data.

**Remark 1.** The assumption of constant density corresponds to a simplified version of the regenerative circuit as it is represented with an increase in the temperature of the fuel LH2, and a resistive term in the mass flow equation. In Liu, Zhang, Wei, and Haidn (2023) a RANS based numerical simulation of a combustion chamber is proposed, to model the heat flux circulation in the regenerative circuit of an LPRE. The main result shows that the heat flow has a strong impact near the cooling walls, but does not have a noticeable impact on the main flow.

The different components of the LPRE can be categorized following the physical domain to which they belong. Three main domains are derived (as the thermic domain is considered static):

* The mechanical domain, containing the turbopump's evolution.

* The hydrodynamical domain, containing the mass flow and pressures evolution.

* The chemical domain, containing the combustion's evolution.

The main variables used in the description of the engine are described in Table 1.
<table><caption>**Table 1**

Definition of main variables.</caption><thead><tr><th>Name</th><th>Description</th></tr></thead><tbody><tr><td>γ</td><td>specific heat ratio [–]</td></tr><tr><td>ω</td><td>rotational speed [rad s<sup>-1</sup>]</td></tr><tr><td>ρ</td><td>density [kg m<sup>-3</sup>]</td></tr><tr><td>Θ</td><td>Temperature [K]</td></tr><tr><td>A<sub>v</sub></td><td>Surface [m<sup>2</sup>]</td></tr><tr><td>h</td><td>enthalpy [J kq<sup>-1</sup>]</td></tr><tr><td>I</td><td>hydraulic inertia [m<sup>-1</sup>]</td></tr><tr><td>J<sub>i</sub></td><td>angular inertia [kg m<sup>2</sup>]</td></tr><tr><td>ρ</td><td>pressure [Pa]</td></tr><tr><td>ṁ</td><td>mass flow [kg s<sup>-1</sup>]</td></tr><tr><td>R<sub>g</sub></td><td>gas constant [J kg<sup>-1</sup> K<sup>-1</sup>]</td></tr><tr><td>T</td><td>torque [N m]</td></tr><tr><td>W</td><td>work [J]</td></tr><tr><td>Z<sub>res</sub></td><td>equivalent resistive coefficient [kg<sup>-1</sup>m<sup>-1</sup>]</td></tr></tbody></table>

## 2.1. System equations

This subsection presents the main equations of the LPRE model. Here the symmetry of the hydrogen and oxygen lines is used, and the hydrogen line is presented. The subscripts are omitted for simplification.

### 2.1.1. Mechanical domain

The mechanical domain is governed by the conservation of the momentum equation,

$$J_i \dot{\omega} = \sum_i T_i \qquad (1)$$

where $J_i$ is the total inertia of the turbo-pump rotational elements, $T_i$ are the exterior torque applied to the turbo-pump and $\omega$ is the rotational speed. Two exterior torques are applied to the system. A resistive torque $T_r$ results of the mass flow action on the pump,

$$T_r = - \left| \frac{a_c}{\rho} \dot{m}_{tot}^2 + b_c \omega \dot{m}_{tot} + c_c \rho \omega^2 \right|, \qquad (2)$$

where $a_c, b_c, c_c$ are manufacturer parameters, $\rho$ is the density of the fluid, and $\dot{m}_{tot}$ is the total mass flow in the pump,

$$\dot{m}_{tot} = \dot{m}_{cc} + \dot{m}_{gg}, \qquad (3)$$

where $\dot{m}_{cc}$ is the main chamber line mass flow and $\dot{m}_{gg}$ is the gas generator line mass flow. The motor torque is described by the equations

$$\begin{align*} N_R &= \frac{\omega d_{tu}}{\sqrt{\gamma R_g \Theta}}, \\ \dot{m}_{out} &= \frac{p_{tu} A_{tu}}{\gamma} \sqrt{\frac{R_g \Theta}{\gamma} \left(\frac{2}{\gamma+1}\right)^{-\frac{\gamma+1}{\gamma-1}}}, \\ W &= \dot{m}_{out} \frac{\sqrt{\gamma R_g \Theta}}{\gamma-1} \left(1 - \pi_T^{\gamma}\right) d_{tu}, \\ ST &= (a_1 + N_R(a_2 + N_R a_3) + \pi_T(a_4 + \pi_T a_5), \\ & \quad + a_6 \pi_T N_R + a_7 \ln \pi_T + a_8 \ln N_R). \text{Corr}, \\ T_m &= ST \cdot W. \end{align*} \qquad (4)$$

where $N_R$ is equivalent to a reduced rotational speed, $W$ is the mass flow work on the turbine and $ST$ denotes the specific torque of the turbine. Additionally, $d_{tu}$, $A_{tu}$, $\pi_{tu}$ designate the diameter of the turbine, the intake area and the area ratio of the turbine. Coefficients $a_1, ..., a_8$ are manufacturer parameters and $Corr$ is a correlation factor computed in Perez Roca (2020). $p_{tu}$ is the pressure of the hot gases in the turbine intake. $\gamma$ is the heat capacity ratio, $R$ is the molar gas constant and $\Theta$ is the temperature of the fluid. In the following, this expression for the motor torque is denoted by

$$T_m = p_{tu}ST(\omega)w_{tu}, \qquad (5)$$

$$\text{where } w_{tu} = \frac{A_{tu}}{\gamma} \sqrt{\frac{R_g \Theta}{\gamma} \left(\frac{2}{\gamma+1}\right)^{-\frac{\gamma+1}{\gamma-1}}}. $$

### 2.1.2. Hydrodynamical domain

The hydrodynamical domain contains the pressure's evolution due to hydrodynamics phenomena and the mass flow's evolution. The mass flow $\dot{m}$'s evolution is governed by the momentum conservation equation,

$$\dot{m} = \frac{1}{I} (p_{in} - p_{out} - Z_{res pipe} \dot{m}^2) \qquad (6)$$

where $p_{in}, p_{out}$ are respectively the up and downstream pressures, $Z_{res pipe}$ is the equivalent resistance of the pipe, defined by the sum of the resistive inertia in a line, and $I$ is the inertia of the pipe. In particular, for a line containing a valve, $Z_{res pipe} = Z_{res line} + Z_{res valve}$. For an LPRE controlled by valve apertures, the term $Z_{res valve} = \frac{1}{2\rho A^2}$ is considered as the control variable. In the case of the pump outlet pressure, it is modeled using manufacturer data as

$$p_{sp} = - \left| \frac{a_p}{\rho} \dot{m}_{tot}^2 + b_p \omega \dot{m}_{tot} + c_p \rho \omega^2 \right|. \qquad (7)$$

The dynamics of a pressure $p$ follows the energy conservation equation,

$$\frac{\dot{p}V}{\gamma - 1} = (h_{in} \dot{m}_{in} - h_{out} \dot{m}_{out}), \qquad (8)$$

where $V$ is the volume of the cavity (or chamber). The hydrodynamical enthalpy consists in the kinetic energy of the entering and exiting flows. In the case of the main combustion chamber, as the exiting mass flow is not considered, it is possible to approximate it using a choked nozzle equation (see Pérez-Roca, Marzat, Piet-Lahanier, Langlois, Galeotta, Farago, and Le Gonidec 2020). In this work, the generic equation for the pressure is finally obtained using first order Taylor expansion (Perez Roca, 2020) of Eq. (8),

$$\dot{p}_{\bullet} = (k_{4,\bullet} + k_{5,\bullet} + (k_{2,\bullet} - k_{1,\bullet})\Theta_{\bullet})\dot{m}_{in} - k_{3,\bullet}\sqrt{\Theta_{\bullet}}p_{\bullet}, \qquad (9)$$

where the coefficient $k_{2,\bullet}, k_{3,\bullet}$ represent the pressure variations due to the volume of gas entering or exiting the chamber, $k_{4,\bullet}, k_{5,\bullet}$ represent the consumption and creation of species resulting from the combustion and $k_{1,\bullet}$ represents the enthalpy loss due to the heating of the species in the chamber. Note that for the turbine pressure, all combustion induced terms are set to 0, since no combustion takes place in the turbine intake, therefore the coefficients $k_{1,\bullet}, k_{4,\bullet}, k_{5,\bullet}$ are equal to 0,

$$\dot{p}_{tu} = (k_{2tu}\Theta)\dot{m}_{tu} - k_{3tu}\sqrt{\Theta_{gg}}p_{tu}. \qquad (10)$$

### 2.1.3. Chemical domain

The chemical domain contains the evolution of the pressure resulting from the combustion. Here the evolution of the chamber pressures is considered, which follows Eq. (8). Two phenomena are to be considered in the combustion process, the enthalpy difference between the reacting species and the products in the chemical reaction (we suppose that the reaction is perfect, meaning that all species are consumed and the reaction is instantly starting) and the loss of enthalpy due to the heating of the reactive species.

## 2.2. Stability analysis

The objective of this paper is to describe how the reformulation of the LPRE using Port-Hamiltonian representation allows to analyze the stability of the LPRE system described in Section 2.1, and provide a design method for a stabilizing controller. In this subsection the stability properties aimed for in the study are briefly recalled.

Stability analysis of non-linear systems often relies on the definition of a Lyapunov function for the system. If a Lyapunov function can be found for the system's equilibrium, then this equilibrium is globally stable. In the case of an LPRE, determining a Lyapunov function from the state equations previously described in Section 2.1 proves to be a difficult task (Perez Roca, 2020). An alternative to requiring a Lyapunov function is to use passivity properties of a system. The
passivity analysis of a system aims to characterize the energy increase of a system compared to the energy flow delivered to the system. From here, an energetic approach is proposed, originating from the simple observation that an LPRE is a physically dissipative (Khalil, 2002) system.

A system

$$ 
\begin{aligned}
\dot{x} &= f(x, u) \\
y &= h(x, u),
\end{aligned}
\quad (11) $$

is dissipative with supply rate $w(u, y)$ if there exists a storage function $S(x)$, $S(0) = 0$, such that for all $x$,

$$ S(x) \ge 0 \text{ and } S(x(T)) - S(x(0)) \le \int_0^T w(u(t), y(t)). \quad (12) $$

A system is passive if it is dissipative with supply rate $w(u, y) = u^\top y$. If the storage function $S$ is differentiable, then Eq. (12) for a passive system can be written,

$$ \dot{S}(x(t)) \le u^\top y. \quad (13) $$

Moreover, for $u = 0$ in (13), $S(x)$ is a candidate Lyapunov function, and if $x = 0$ is an equilibrium of the system, then $x = 0$ is globally stable. A storage function is in practice difficult to obtain for an LPRE, therefore a reformulation of the state-space into a Port-Hamiltonian system is performed to obtain it.

# 3. Port-Hamiltonian systems

## 3.1. Introduction to Port-Hamiltonian systems

Let a system represented by the state-space equation,

$$ 
\begin{aligned}
\dot{x} &= f(x) + g(x)u \\
y &= h(x, u).
\end{aligned}
\quad (14) $$

Port-Hamiltonian modeling aims at transforming the system (14) under the form

$$ 
\begin{aligned}
\dot{x} &= (J(x) - R(x)) \frac{\partial H}{\partial x}(x) + G(x)u \\
y &= G^\top(x)
\end{aligned}
\quad (15) $$

where the different terms composing the equation are explained in this section.

Port-Hamiltonian systems framework, described in Maschke and van der Schaft (1993) and Duindam et al. (2009), derived from bond graphs theory (Breedveld, 1984), makes use of physical invariants of a system, such as the total energy, to highlight the energy transfers between the different elements composing a system. Originally established for mechanical systems, Port-Hamiltonian framework has been extended to various systems such as hydrodynamics in Pasumarthy and van der Schaft (2006) or chemical processes in Eberard and Maschke (2004), Eberard et al. (2007), Ramirez, Maschke, and Sbarbaro (2013). Port-Hamiltonian theory uses the elementary interactions between components to reconstruct the equations that describe a system's behavior. Interactions between elements follow Dirac structures (Duindam et al., 2009) which respect a power balance property (see (16)). The different nodes of the structures contain two connecting elements: the effort $e$ and the flow $f$ that represent coupled physical quantities (e.g. pressure and flow). Elements can be classified into three distinct functions:

* Energy-storing, denoted by $e_s$, $f_s$, are elements that store energy (kinetic, potential, ...).

* Energy-dissipating, denoted by $e_r$, $f_r$ dissipate the energy, using friction for example.

* Energy-routing, denoted by $e_t$, $f_t$ transfers energy, with no storage nor loss. These elements are by construction perfect transformers, therefore are seldom used in real systems.

Additionally, interactions with the exterior environment are denoted by $e_p$, $f_p$. The product $e^\top f$ is the power through a node defined by the flow/effort couple $(e, f)$. To achieve a Dirac structure $D$, the sum of the powers involved in $D$ (e.g. $\forall(e, f) \in D$) is equal to 0. This requirement implies that the energy balance of a Port Hamiltonian system is

$$ e_s^\top f_s + e_r^\top f_r + e_p^\top f_p = 0. \quad (16) $$

where $f_s$ is taken equal to the rate of change of the states with a minus sign $-\dot{x}$ per convention. By definition, the dissipative elements yield $e_r^\top f_r \le 0$, therefore

$$ -e_s^\top f_s \le e_p^\top f_p. \quad (17) $$

Eq. (17) is a classical equation for the passivity of a system where the system energetic input $u$ is denoted by $e_p$ and the output $y$ is $f_p$ (e.g. $-e_s^\top f_s \le u^\top y$). In a Port-Hamiltonian system where the state is denoted by $x$, let the function $H(x)$ be defined by

$$ \frac{dH}{dt} = -e_s^\top f_s = \frac{\partial H}{\partial x} \dot{x}. \quad (18) $$

The function $H$ is called the Hamiltonian of the system. $H$ represents the total energy of the system by definition as its derivatives correspond to the stored power in the system. Using (17) and for $H(x) > 0$ for all $x$, $H$ is a candidate storage function for the passivity of a system (see Duindam et al. 2009,2.6).

From these preliminary definitions, a Port Hamiltonian system is defined under the structure of the form (15), where the matrix $J(x)$ expresses the energy storage terms, with $J(x) = -J(x)^\top$ and $R(x)$ expresses the energy dissipation with $R(x) \ge 0$.

**Remark 2.** In practice, the expression of a system under Port-Hamiltonian form differs from the classical physical equations from which the state equations are usually derived. To obtain the fixed geometric structure of Port Hamiltonian systems a reformulation of the system is often required. However, such reformulation yields the Hamiltonian function, which is a powerful tool for the analysis of a system, justifies the need of a new modeling tool.

## 3.2. Interconnection of Port-Hamiltonian systems for multi-physical systems

In multi-physical systems, the different subsystems are expressed independently, however it is necessary to interconnect the individual subsystems to build the complete model in Port-Hamiltonian framework (Cervera, van der Schaft, & Baños, 2007; Duindam et al., 2009). Let two Port-Hamiltonian subsystems

$$ 
\begin{aligned}
\dot{x}_1 &= (J_1 - R_1) \frac{\partial H_1}{\partial x_1} + G_1(x_1)x_2 \\
\dot{x}_2 &= (J_2 - R_2) \frac{\partial H_2}{\partial x_2} + G_2(x_2)x_1,
\end{aligned}
\quad (19) $$

where $J_1, R_1, H_1, J_2, R_2, H_2$ are respectively the storage, dissipation matrices and Hamiltonian function of the considered systems and $G_1, G_2$ are the exterior interaction matrices. The following interconnection of the subsystems in (19) is considered,

$$ \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} = ([J'] - [R']) \begin{bmatrix} \frac{\partial H_1}{\partial x_1} \\ \frac{\partial H_2}{\partial x_2} \end{bmatrix}. \quad (20) $$

From (20), it is necessary to find $J'$ and $R'$ in (20) to conserve the equations in (19). The method presented in Cervera et al. (2007) proposes to use shared space of variables. To interconnect the turbo-pump and the hydrodynamical system, this method is depicted in Section 4.2.

The author in Cervera et al. (2007) furthermore proves that the composition of two Dirac structures remains a Dirac structure, and that the interconnection remains a Port Hamiltonian system, with new Hamiltonian,

$$ H'(x) = H_1(x_1) + H_2(x_2). \quad (21) $$
The properties of interconnections for Port Hamiltonian systems enable the user to build subsystems separately and interconnect them in the following, if a shared space of variables is found. In the practical, the construction of the shared space is realized using equivalent quantities, take for example the reformulation of a torque into an equivalent pressure in Esquivel-Sancho et al. (2022), where the rotational speed is formulated into an equivalent mass flow.

In the next section, the Port Hamiltonian representation of the LPRE model is depicted, first by describing the subsystems that compose each physical domain, and then by interconnecting each subsystem.

# 4. Port-Hamiltonian modeling of an LPRE

## 4.1. Presentation of the subsystems

The different physical domains can be reformulated into several Port-Hamiltonian subsystems which can be deduced from their state-space formulation. In the literature, the couples effort/flow have been defined for several physical domains (Rashad et al., 2020). In this subsection, Port Hamiltonian models are proposed for the three physical domains of the engine.

Mechanical equations make use of the coupled quantities rotational speed / torque (see for example Aoues, Cardoso-Ribeiro, Matignon, & Alazard, 2017). In the case of the LPRE, (1) translates into

$$J_{\dot{\omega}} = -T_r + T_m \qquad (22)$$

where both torques are given using (2) and (4) and are considered exterior interactions with other subsystems.

Hydrodynamical equations involve the coupled quantities volume flow / pressure (Esquivel-Sancho et al., 2022). For the sake of conserving the hydrodynamical Eq. (6), the coupled quantity in this paper will be taken as mass flow/pressure. Additionally, the pressure evolution can be fully written under hydrodynamical form and (6) translates into

$$\begin{bmatrix} I \dot{m} \\ k_p \dot{p} \end{bmatrix} = \begin{bmatrix} -R_m & -1 \\ 1 & -R_p \end{bmatrix} \begin{bmatrix} \dot{m} \\ p \end{bmatrix} + \begin{bmatrix} p_{ep} \\ 0 \end{bmatrix} \qquad (23)$$

where the different terms can be expressed following (6), (8), and yield $R_m = Z_{res}J_{\dot{m}}$, $k_p = \frac{V}{(\gamma-1)h_{tu}}$. The resistive term $R_p$ differs depending on the pressure considered.

## 4.2. Interconnection of the subsystems and Hamiltonian

The interconnection of the hydrodynamical part of the system can be driven directly, by corresponding both flow and efforts of the individual subsystems, as these coupled quantities are already part of a shared space of variables (e.g. the mass flows and pressures). Furthermore, the chemical equation uses the same space of variables as the hydrodynamical Eqs. (6), (8). Therefore the chemical domain can be interconnected to the hydrodynamical domain directly.

The mechanical system requires a distinct space of variables (e.g. the rotational speed and torque) and cannot be connected directly to the hydrodynamical domain (which uses pressures and mass flows). A proposition driven in Esquivel-Sancho et al. (2022) is to formulate the torque on the mechanical end into an equivalent pressure. This results in equating the output pressure of the pump with the resistive torque equivalent pressure. However, the main hypothesis of the reversible hydro-pump system in Esquivel-Sancho et al. (2022) is a linear relationship between the mass flow and the rotational speed $\dot{m} = k_{q\omega}\omega$, and that both the pressure increase $p_{sp}$ and the resistive torque $T_r$ solely depend on respectively $\dot{m}$ and $\omega$, e.g. $p_{sp} = k_{psp}q^2$, $T_r = k_{T_r}\omega^2 = k_{T_r}k_{m\omega}\dot{m}^2$. As a consequence, the resistive torque and the pressure increase can be expressed using the mass flow,

$$T_r = \frac{k_{m\omega}k_{T_r}}{k_{psp}} p_{sp} \qquad (24)$$

In this work, the turbo pump studied presents no such simplification, as the resistive torque and the pressure increase follow distinct polynomial functions:

$$\begin{aligned} p_{sp} &= - \left| \frac{a_p}{\rho} \dot{m}_{tot}^2 + b_p \omega \dot{m}_{tot} + c_p \rho \omega^2 \right|, \\ T_r &= - \left| \frac{a_c}{\rho} \dot{m}_{tot}^2 + b_c \omega \dot{m}_{tot} + c_c \rho \omega^2 \right|. \end{aligned} \qquad (25)$$

To face this issue, it is possible to identify the driving terms of the polynomial functions to approximate an equivalent pressure. Consider a lower value for the resistive torque, expressed as $\frac{a_c}{\rho}\dot{m}^2$ and an upper value for the pressure increase, expressed as $b_p\dot{m}\omega$. Define

$$\alpha_q = \frac{b_p \rho}{a_c}, \qquad (26)$$

where the coefficients $b_p$, $a_c$ are taken as in (2) and (7). With the notations $\bar{p}_{sp} = b_p\dot{m}\omega$ and the resistive torque value $\underline{T}_r = \frac{1}{a_m}b_p\dot{m}^2$, the torque and pressure increase are rewritten

$$\begin{aligned} p_{sp} &= \bar{p}_{sp} - \left(\frac{a_p}{\rho}\dot{m} - c_p\rho\frac{\omega^2}{\dot{m}}\right)\dot{m} + p_{ep} \\ T_r &= \underline{T}_r - \alpha_q(b_c\dot{m} + c_c\omega)\omega. \end{aligned}$$

Then the equations of $\omega$ and a mass flow $\dot{m}$ through the pump read:

$$\begin{aligned} \alpha_m J \dot{\omega} &= -b_p \dot{m}^2 - \alpha_m (b_c \dot{m} + c_c \omega) \omega + \alpha_m T_m \\ I \dot{m} &= b_p \dot{m} \omega - (Z_{eq} \dot{m} - \frac{a_p}{\rho} \dot{m} - c_p \rho \frac{\omega^2}{\dot{m}}) \dot{m} + p_{ep}, \end{aligned}$$

which gives the following interconnection:

$$\begin{bmatrix} \alpha_m J \dot{\omega} \\ I \dot{m} \end{bmatrix} = \begin{bmatrix} -R_{\omega} & -b_p \dot{m} \\ b_p \dot{m} & -R_m \dot{m} \end{bmatrix} \begin{bmatrix} \omega \\ \dot{m} \end{bmatrix} + \begin{bmatrix} \alpha_m T_m \\ p_{ep} \end{bmatrix}, \qquad (27)$$

where $R_q = Z_{eq}q - \frac{a_p}{\rho}\dot{m} - c_p\rho\frac{\omega^2}{\dot{m}}$ and $R_{\omega} = \alpha_m(b_c\dot{m} + c_c\omega)$ have to be positive to satisfy the upper and lower value hypothesis. Here two conditions are derived for the validity of the Port-Hamiltonian formulation of the system (namely $R_m > 0$, $R_{\omega} > 0$). These conditions can be used as validation test for the stability of the system.

For the turbine part, a new shared space of variables is proposed. Let the effort $\alpha_m T_m = \alpha_m ST(\omega)w_{tu}p_{tu}$ (27). Suppose that a minimal value for the pressure loss in the turbine is $\alpha_m ST(\omega)w_{tu}\omega$. A proposition to close the system in Port-Hamiltonian framework is to reformulate the term $\sqrt{\Theta_{gg}k_{3tu}p_{tu}}$ into:

$$\sqrt{\Theta_{gg} k_{3tu} p_{tu}} = \alpha_t \omega + \beta_t p_{tu}, \qquad (28)$$

where $\alpha_t = \alpha_m ST(\omega)w_{tu}$ and $\beta_t$ is computed from $\sqrt{\Theta_{gg}k_{3tu}p_{tu}} - \alpha_t\omega$. Hence the following interconnection:

$$\begin{bmatrix} \alpha_m J_t \dot{\omega} \\ k_{ptu} \dot{p}_{tu} \end{bmatrix} = \begin{bmatrix} 0 & \alpha_t \\ -\alpha_t & -R_{ptu} \end{bmatrix} \begin{bmatrix} \omega \\ p_{tu} \end{bmatrix}. \qquad (29)$$

Here $\alpha_t\omega$ is assumed to be a lower value of the term $k_{3tu}p_{tu}$. The term $\beta_t p_{tu}$ has to be positive, with $R_{ptu} = \beta_t \ge 0$. This leads to a third condition on the system for the Port-Hamiltonian formulation.

## 4.3. Expression of the Port-Hamiltonian system

Using the subsystems developed from (22) to (29), the resulting model of the LPRE as a Port-Hamiltonian system expresses as

$$\begin{bmatrix} \alpha_m J_t \dot{\omega} \\ I_{cc} \dot{m}_{cc} \\ I_{gg} \dot{m}_{gg} \\ I_{tu} \kappa_{tu} \dot{m}_{tu} \\ k_{pcc} \dot{p}_{cc} \\ k_{pgg} \dot{p}_{gg} \\ k_{ptu} \dot{p}_{tu} \end{bmatrix} = [J(X) - R(X)] \begin{bmatrix} \omega \\ \dot{m}_{cc} \\ \dot{m}_{gg} \\ \dot{m}_{tu} \\ p_{cc} \\ p_{gg} \\ p_{tu} \end{bmatrix} + \begin{bmatrix} 0 \\ 0 \\ -\dot{m}_{gg}^2 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} Z_{vgg}, \qquad (30)$$
where the matrices J and R express as

$$
J(X) = \begin{bmatrix}
0 & -\lambda & -\lambda & 0 & 0 & 0 & \alpha_t \\
\lambda & 0 & 0 & 0 & -1 & 0 & 0 \\
\lambda & 0 & 0 & 0 & 0 & -1 & 0 \\
0 & 0 & 0 & 0 & 0 & \kappa & -\kappa \\
0 & 1 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & -\kappa & 0 & 0 & 0 \\
-\alpha_t & 0 & 0 & \kappa & 0 & 0 & 0
\end{bmatrix}
$$

where $\lambda = b_p(\dot{m}_{cc} + \dot{m}_{gg})$, and

$$
\mathrm{diag}([R_{\omega} \quad R_{mcc} \quad R_{mgg} \quad R_{mtu} \quad R_{pcc} \quad R_{pgg} \quad R_{ptu}]). \quad R(X) =
$$

where the resistive terms are computed with

$$
R_{\omega} = \alpha_{\dot{m}}(b_{c}(\dot{m}_{cc} + \dot{m}_{gg}) - c_{c}\rho\omega)
$$

$$
R_{mcc} = Z_{eqcc} \dot{m}_{cc} - \frac{a_p}{\rho} (\dot{m}_{cc} + \dot{m}_{gg}) - c_p \rho \frac{\omega^2}{\dot{m}_{cc}}
$$

$$
R_{mgg} = Z_{eqgg} \dot{m}_{gg} - \frac{a_p}{\rho} (\dot{m}_{cc} + \dot{m}_{gg}) - c_p \rho \frac{\omega^2}{\dot{m}_{gg}}
$$

$$
R_{mtu} = \kappa_{tu} Z_{hg} \dot{m}_{tu}
$$

$$
R_{pcc} = \sqrt{\Theta_{cc}} k_{3cc} k_{pcc}
$$

$$
R_{pgg} = 0
$$

$$
R_{ptu} = \beta_t
$$

Finally, the constant scaling terms of the system can be approximated
with the following

$$
\kappa_{tu} = \sqrt{\Theta_{gg}} k_{3gg} k_{pgg}
$$

$$
k_{pcc} = \frac{1}{k_{4cc} + k_{5cc} + (k_{2cc} - k_{1cc})\Theta_{cc}}
$$

$$
k_{pgg} = \frac{1}{k_{4gg} + k_{5gg} + (k_{2gg} - k_{1gg})\Theta_{gg}}
$$

$$
k_{ptu} = \frac{\kappa_{tu}}{k_{2tu}\Theta_{gg}}
$$

Expression of the interconnection terms:

$$
\alpha_{\dot{m}} = \frac{b_p \rho}{a_c}, \quad \alpha_t = \alpha_{\dot{m}} S T(\omega) w_{tu}
$$

4.4. Expression of the Hamiltonian function

Using the Port-Hamiltonian model equations, a Hamiltonian func-
tion is derived directly using (18).

$$
H(x) = \frac{1}{2} J_t \omega^2 \alpha_m + \frac{1}{2} I_{cc} \dot{m}_{cc}^2 + \frac{1}{2} I_{gg} \dot{m}_{gg}^2 + \frac{1}{2} I_{tu} \dot{m}_{tu}^2 \\
+ \frac{1}{2} k_{pcc} p_{cc}^2 + \frac{1}{2} k_{pgg} p_{gg}^2 + \frac{1}{2} k_{ptu} p_{tu}^2. \tag{31}
$$

Remark 3. Given the formulation of the Hamiltonian function (31), the condition $H(x) > 0$ is verified. From Eqs. (17), (18), the Hamiltonian function verifies

$$
H(x) \le u^T y \tag{32}
$$

which proves that the Hamiltonian *H* is a candidate function for
the passivity, and the LPRE system is stable at the minimum of the
Hamiltonian *H*.

The Hamiltonian function is minimal at the origin where all the
states are at 0 and therefore the stability of the system is only proven at
this equilibrium. In order to obtain stability properties for the system at
a desired equilibrium, a controller is designed in the following section.
The controller ensures a passive behavior around a desired functioning
point and provides a new Hamiltonian function, proving stability of the
engine during the permanent regime at a required thrust.

5. Passivity-based controller

A passivity-based controller (PBC) makes use of the properties of passive systems to stabilize a system. In particular, such controller makes the closed-loop system passive regarding the energetic input of the system. A passivity-based controller is proposed in the following to regulate the LPRE system in (30), around a designed functioning point as the equilibrium for the closed-loop system.

5.1. Passivity-based control theory

PBC (Ortega, Van Der Schaft, Maschke, & Escobar, 2002) is a
classical control design method in Port-Hamiltonian systems theory. It
is briefly recalled here.

Let a Port-Hamiltonian system of the form

$$
\dot{x} = (J - R) \frac{\partial H}{\partial x} + g.u(x, t), \quad (33)
$$

with the associated desired system behavior given by

$$
\dot{x}_d = (J_d - R_d) \frac{\partial H_d}{\partial x}, \quad (34)
$$

where the matrices $J_d = -J_d^\top$, $R_d \ge 0$ represent the desired storage and damping matrices, and $H_d$ is the resulting Hamiltonian function.

The principle of this method is to define a new stable equilibrium for the system, by assigning a new total energy represented by the function $H_d(x)$. This new energy function should be minimal at the desired equilibrium and consequently be a candidate function for the passivity of the new closed-loop system, on the one hand with $H_d(x) \ge 0$ and on the other hand by conserving the Port-Hamiltonian structure $J_d(x) = -J_d(x)^\top$, $R_d(x) \ge 0$. While $J_d$ and $R_d$ are defined by the user, the Hamiltonian function $H_d$ is not known a priori.

The developed method consists in searching for a controller *u* such that the dynamics of the closed-loop system (33) converges to the dynamics of (34). This translates into

$$
(J - R) \frac{\partial H}{\partial x} + gu = (J_d - R_d) \frac{\partial H_d}{\partial x}, \quad (35)
$$

which is the generation equation for the controller. First, it is possible
to compute the new Hamiltonian function using (35), where the control
action is nullified,

$$
g^{\perp}(J_d - R_d) \frac{\partial H_d}{\partial x} = g^{\perp}(J - R) \frac{\partial H}{\partial x} \quad (36)
$$

where $g^{\perp}$ is a full rank left annihilator of $g$, i.e. $g^{\perp}g = 0$.

Assuming we can find $u(x)$, $R_a(x)$, $J_a(x)$ and $\frac{\partial H_a}{\partial x}$, with properties specified in the following, the closed loop system (33) is a Port Hamiltonian system with dissipation $R_d(x)$ and a stable equilibrium at a specified $x = x_0$ (see Ortega, Van Der Schaft, Castanos, and Astolfi 2008, Prop. 1.)

First define

$$
\begin{align}
R_a(x) &= R_d(x) - R(x), \quad R_a(x) \ge 0 \nonumber \\
J_a(x) &= J_d(x) - J(x), \quad J_a(x) = -J_a(x)^\top. \tag{37}
\end{align}
$$

$H_a(x)$ is defined by $H_a(x) = H_d(x) - H(x)$. While the knowledge of $J_a$, $R_a$ is granted from $J_d$, $R_d$, $H_a(x)$ is to be determined (equivalent to finding $H_d(x)$). $H_a(x)$ is the Hamiltonian function added to the system that must respect the following existence conditions,

* $\frac{\partial H_a}{\partial x}$ is the gradient of a scalar function, and should be integrable, therefore $H_a(x)$ must respect an integrability condition

$$
\frac{\partial^2 H_a}{\partial x^2} = \left[ \frac{\partial^2 H_a}{\partial x^2} \right]^\top . \qquad (38)
$$

* The equilibrium assignment at $x = x_0$ defines the minimum of $H_d(x)$. In other terms, $\min_x(H_d(x)) = H_d(x_0)$. Paired with the passivity conditions ($H_d(x) \ge 0$, see Section 3.1), the minimal
value of $H_d$ is set to $H_d(x_0) = 0$. The function $H_a(x)$ is then defined at $x_0$ by

$$ 
\begin{aligned}
H_a(x_0) &= H(x_0) \\
\frac{\partial^2 H_a}{\partial x^2}(x_0) &> 0.
\end{aligned}
\quad (39) $$

**Remark 4.** $H_a(x)$ is used as an intermediary to compute the Hamiltonian function $H_d(x)$. $H_a(x)$ introduces an integrability condition (39), in addition to the partial derivative equation (PDE) on $H_d(x)$ (36). The integrability condition introduces enough equations to solve the PDE.

Following the definitions of $H_a$, $R_a$, $J_a$, (36) can be written

$$ g^{\perp}(J + J_a - R - R_a) \frac{\partial H_a}{\partial x} = -g^{\perp}(J_a - R - a) \frac{\partial H}{\partial x} \quad (40) $$

The controller $u(x)$ is defined directly using

$$ u = (g^{\top} g)^{-1} g^{\top} \left( (J_d - R_d) \frac{\partial H_d}{\partial x} - (J - R) \frac{\partial H}{\partial x} \right) \quad (41) $$

where $(g^{\top} g)^{-1} g^{\top}$ is a left inverse of the control matrix $g$.

## 5.2. Application to a simplified LPRE model

To illustrate the design of the passivity-based controller, we present it using a simplified model of an LPRE, which has been depicted in Gibart et al. (2023) and using motopumps instead of turbo-pumps. Let the system

$$ \begin{bmatrix} J_i \alpha_m \dot{\omega} \\ I \dot{m} \\ k_p \dot{p} \end{bmatrix} = (J(x) - R(x)) \frac{\partial H}{\partial x} + g(x) T_m + \begin{bmatrix} 0 \\ 0 \\ p_{ep} \end{bmatrix}, \quad (42) $$

where

$$ J(x) = \begin{bmatrix} 0 & -\lambda & 0 \\ \lambda & 0 & -1 \\ 0 & 1 & 0 \end{bmatrix}, \quad R(x) = \begin{bmatrix} -R_{\omega} & 0 & 0 \\ 0 & -R_{\dot{m}} & 0 \\ 0 & 0 & -R_p \end{bmatrix} \quad (43) $$

$$ \text{and } g(x) = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}. \quad (44) $$

### 5.2.1. Preliminaries

The application of PBC to the simplified LPRE model aims to stabilize the system to an equilibrium that differs from the origin of the system ($x_0 \neq 0$) and retrieve a new Hamiltonian function $H_d$ for the system. A first step for the design of the PBC controller is to establish the desired effect of the controller on the system. The equilibrium is defined by a known chamber pressure $p_0$, for which the corresponding $\omega_0$, $\dot{m}_0$ can be computed (therefore are supposed known). As the initial Hamiltonian function of an LPRE is affine in the states $x$ in (31), the Hamiltonian $H_d$ is chosen under an affine form,

$$ H_d(x) = x^{\top} Q x \quad (45) $$

where $Q$ is a positive symmetric square matrix that represents the equivalent inertia of the system. Additionally, let the state vector be of the form,

$$ x_d = \begin{bmatrix} J_i \alpha_m \omega - J_i \alpha_m \omega_0 \\ I \dot{m} - I \dot{m}_0 \\ k_p p - k_p p_0 \end{bmatrix}, \quad (46) $$

for an equilibrium defined by $x_0 = [J_i \alpha_m \omega_0 \ I q_0 \ k_p p_0]^\top$. The state vector at the equilibrium is then defined by $x_d(x_0) = [0 \ 0 \ 0]^\top$. In this first approach, define $J_d = J$ and $R_d = R$, (leading to $J_a = R_a = 0$) as there are no other requirements than equilibrium assignment.

### 5.2.2. Developed approach for the design

The first step to design the controller is the determination of the Hamiltonian function $H_a(x)$ or equivalently the matrix $Q$ in Eq. (45). This step makes use of Eqs. (38), (40) and (45) to express the Hamiltonian function.

Let a full rank left annihilator of $g$ be

$$ g^{\perp} = \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}. \quad (47) $$

Eq. (40) yields

$$ \begin{bmatrix} \lambda & -R_{\dot{m}} & -1 \\ 0 & 1 & -R_p \end{bmatrix} \frac{\partial H_a}{\partial x} = 0, \quad (48) $$

which translates into two equations for the partial derivatives of $H_a$,

$$ \begin{aligned} \frac{\partial H_a}{\partial I \dot{m}} &= R_p \frac{\partial H_a}{\partial k_p p} \\ \frac{\partial H_a}{\partial J_i \alpha_m \omega} &= \frac{(R_{\dot{m}} R_p + 1)}{\lambda} \frac{\partial H_a}{\partial k_p p}. \end{aligned} \quad (49) $$

The new Hamiltonian $H_d$ is proposed with an affine form in (45). Additionally, (49) implies cross terms within the different components. Then $H_d$'s partial derivatives can be written,

$$ \begin{aligned} \frac{\partial H_d}{\partial J_i \alpha_m \omega} &= Q_{\omega}(\omega - \omega_0) + k_{\omega \dot{m}}(\dot{m} - \dot{m}_0) + k_{\omega p}(p - p_0), \\ \frac{\partial H_d}{\partial I \dot{m}} &= Q_{\dot{m}}(I \dot{m} - I \dot{m}_0) + k_{\dot{m} \omega}(\omega - \omega_0) + k_{\dot{m} p}(p - p_0), \\ \frac{\partial H_d}{\partial k_p p} &= Q_p(p - p_0) + k_{p \omega}(\omega - \omega_0) + k_{p \dot{m}}(\dot{m} - \dot{m}_0), \end{aligned} \quad (50) $$

where $Q$, is the equivalent inertia coefficient for a state and $k_\cdot$ is a cross term coefficient in the inertia matrix $Q$. Here the states are expressed as in (46) to ensure the equilibrium assignment. The matrix $Q$ expresses as:

$$ Q = \begin{bmatrix} Q_{\omega} & k_{\omega \dot{m}} & k_{\omega p} \\ k_{\dot{m} \omega} & Q_{\dot{m}} & k_{\dot{m} p} \\ k_{p \omega} & k_{p \dot{m}} & Q_p \end{bmatrix} \quad (51) $$

The relation $H_a(x) = H_d(x) - H(x)$ is used to translate (50) into

$$ \begin{aligned} \frac{\partial H_a}{\partial J_i \alpha_m \omega} &= (Q_{\omega} - 1)(\omega - \omega_0) + k_{\omega \dot{m}}(\dot{m} - \dot{m}_0) + k_{\omega p}(p - p_0) + \omega_0, \\ \frac{\partial H_a}{\partial I \dot{m}} &= (Q_m - 1)(\dot{m} - \dot{m}_0) + k_{\dot{m} \omega}(\omega - \omega_0) + k_{\dot{m} p}(p - p_0) + \dot{m}_0, \\ \frac{\partial H_a}{\partial k_p p} &= (Q_p - 1)(p - p_0) + k_{p \omega}(\omega - \omega_0) + k_{p \dot{m}}(\dot{m} - \dot{m}_0) + p_0, \end{aligned} \quad (52) $$

The integrability condition in (38) gives

$$ \begin{aligned} k_{\omega \dot{m}} &= k_{\dot{m} \omega}, \\ k_{\omega p} &= k_{p \omega}, \\ k_{\dot{m} p} &= k_{p \dot{m}} \end{aligned} \quad (53) $$

From (49), (52), (53) it is possible to express every unknown term in terms of a single one. Let $Q_{\dot{m}}$ endorse such role. The other terms of the equation expressed as a function of $Q_{\dot{m}}$ read as,

$$ \begin{aligned} Q_{\omega} &= \frac{(R_{\dot{m}} R_p + 1)^2}{\lambda^2 R_p^2} (Q_{\dot{m}} - 1) + 1 \\ k_{\omega \dot{m}} &= \frac{R_{\dot{m}} R_p + 1}{\lambda R_p} (Q_{\dot{m}} - 1) \\ k_{\omega p} &= \frac{R_{\dot{m}} R_p + 1}{\lambda R_p^2} (Q_{\dot{m}} - 1) \\ Q_p &= \frac{1}{R_p^2} (Q_{\dot{m}} - 1) + 1 \\ k_{qp} &= \frac{1}{R_p} (Q_{\dot{m}} - 1) \end{aligned} \quad (54) $$

The equilibrium assignment (39) can be easily verified from (50).
The controller can finally be expressed using (41),

$$u = -R_{\omega} \frac{\partial H_a}{\partial J_i \alpha_m \omega} - \lambda \frac{\partial H_a}{\partial I \dot{m}}, \quad (55)$$

or in a developed form,

$$u(x) = -R_{\omega} \left[ \frac{(R_p R_m + 1)}{\lambda R_p} (Q_m - 1) \left( \frac{(R_p R_m + 1)}{\lambda R_p} (\omega - \omega_0) + (\dot{m} - \dot{m}_0) + \frac{1}{R_p} (p - p_0) \right) + \omega_0 \right] - \lambda \left[ (Q_m - 1) \left( (\dot{m} - \dot{m}_0) + \frac{(R_p R_m + 1)}{\lambda R_p} (\omega - \omega_0) + \frac{1}{R_p} (p - p_0) \right) + \dot{m}_0 \right] \quad (56)$$

### 5.2.3. Controller tuning with the parameter $Q_m$

The tuning parameter of the controller $Q_m$ is derived from conditions on $H_a$ summarized by (39). To ensure a minimum of the Hamiltonian function at the desired equilibrium $x_0$, the parameter $Q_m$ is set to satisfy

$$\begin{bmatrix} Q_{\omega} & k_{\omega \dot{m}} & k_{\omega p} \\ k_{\omega \dot{m}} & Q_{m} & k_{m p} \\ k_{\omega p} & k_{m p} & Q_{p} \end{bmatrix} > 0 \quad (57)$$

A simple choice considering the eigenvalues of (57) is to use $Q_m = 1 + v(x)$ where $v(x)$ is a positive function. The eigenvalues of the matrix $Q$ are

$$\begin{bmatrix} Q_m + \frac{1}{R_p^2 \lambda^2} (Q_m - 1) (R_m^2 + 2R_m R_p + \lambda^2 + 1) \\ 1 \\ 1 \end{bmatrix} \quad (58)$$

Additionally, to obtain the desired behavior the function $v(x)$ can be tuned to scale down the highest eigenvalue of the system (to the other terms of the rotational speed equation). In this system, the term $\frac{1}{R_p^2}$ represents the highest gain of the system ($0 < R_p < 1$). Therefore the equivalent inertia of the Hamiltonian is defined as

$$Q_m = 1 + k_u R_p^2 \quad (59)$$

where the parameter $k_u > 0$ can vary and influences the speed and precision of the controller. Fig. 2 illustrates the behavior of the pressure for different tuning parameters $k_u$ in the range [0.1, 10] in response to a step function. The time-response of the controller diminishes with the value of $k_u$. The requirement for an integrator is illustrated by the static errors that remain for the different parameters.

In this design, the controller obtained yields the closed-loop (34). The introduction of an integral action in the definition of the equilibrium $x_0$ in (46) suppresses the resulting static error while conserving the Port-Hamiltonian structure.

# 6. Application to the complete system and results

## 6.1. Validation of the Port-Hamiltonian model

The first analysis consists in comparing the evolution of the main variables of the LPRE obtained either using the state-space modeling and the Port-Hamiltonian representation. Both models were integrated in Matlab Simulink simulations. The state-space modeling corresponds to the one developed in Perez Roca (2020), where the detailed equations are presented. Figs. 3(a) and 3(b) feature the evolution of the states in response to a step in the actuators value. At time $t = 1$ s, the valves' surfaces spike from $A_{v0}$ to $1.1A_{v0}$. The step is illustrated with the evolution of $A_{v,gg,H}$, $A_{v,gg,O}$ on Fig. 3(a). In this comparison, the actuators dynamics are not modeled, and therefore the step on the actuators is represented by a step on the equivalent resistive inertia $Z_{resgg}$.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Evolution of the pressure for different tuning parameters $Q_m$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Open-loop response to a step function: State space model.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Open-loop response to a step function: Port-Hamiltonian model.</figcaption></figure>

Figs. 4(a) and 4(b) display the evolution of the pressures and the mass flows responding to the same step function with similar initial conditions, obtained using the Port Hamiltonian modeling. One can observe that both models provide the same behavior on the different state components, therefore validating the use of the Port-Hamiltonian model. Other comparisons have been performed using different initial conditions and input values for complementary validation.

## 6.2. Addition of a damping matrix $R_a(x)$

In this section, the results for the application of the controller designed in Section 5 for a complete LPRE system (30) are presented. For the sake of simplicity, the computations to design the controller are not written here.

Similarly to the simplified system, the desired Hamiltonian function $H_d$ for the complete LPRE is computed using the variables separation method as in Zhou, Wu, Hu, Li, and Wang (2021). While for the simplified system, the equation closure is granted for (49), in the case of the complex system, Eq. (36) has no solution for $R_a = J_a = 0$, with the method conducted in Section 5.2. To find a solution for the PDE of
the complex model, it is then necessary to find a new suitable matrix $R_u$. Let a diagonal matrix $R_u \neq 0$,

$$R_u(x) = \text{diag}(R_{a\omega}, R_{aqcc}, R_{a\text{gg}}, R_{a\text{qtu}}, R_{apcc}, R_{apgg}, R_{aptu}). \quad (60)$$

The addition of 7 parameters to Eq. (36) allows the user to solve the equations. In particular, 4 conditions on the added resistive matrix $R_u(x)$ are derived from solving the equation,

$$\begin{aligned} R_{aqcc} &= R_{mcc} \\ R_{apcc} &= 0 \\ R_{a\omega} &= R_{\omega} - R_{pcc} \\ R_{apgg} &= R_{m\omega} - R_{am\omega} \end{aligned} \quad (61)$$

Using the newly defined matrix $R_u(x)$ with conditions stated in (61), the number of independent equations derived from the PDE (36) is now 1 less than the number of parameters, which constitutes a sufficient condition for solving the PDE system. A more powerful method is presented in Ortega et al. (2002), to solve partial derivative equations using the method of characteristics. However, in the case of a separation method, a sufficient condition on the closure of the equations (illustrated by the addition of parameters) is enough to solve the system.

The equations are solved by expressing every parameter as a function of $Q_{m\text{gg}}$, which is the equivalent impedance for the state $\dot{m}_{gg}$ in the desired Hamiltonian. As the evolution of $\dot{m}_{gg}$ is directly affected by the controller in (30), $Q_{gg}$ is chosen as the tuning parameter. In this case the equilibrium assignment condition (39) yields the following parameter $Q_{m\text{gg}}$

$$Q_{m\text{gg}} = 1 + \nu(x, t) \quad (62)$$

where $\nu(x, t)$ remains a positive function.

The following simulation illustrate the step response of the controlled system. Fig. 5, shows the controlled system in closed-loop, responding to a step function. Similarly to the previous example using the simplified model, a proportionate-integral (PI) controller with $k_I = 20, k_P = 1$ has been added to the designed controller to eliminate the static error. The addition of the integrator does not affect the stability of the closed-loop, as the equation of the controlled system expresses as,

$$\dot{x} = (J_d - R_d)\frac{\partial H_d}{\partial x} + G(x)k_I \int_0^t (p_{ccref} - p_{cc})dt. \quad (63)$$

where the Hamiltonian $H_d$ is by construction a storage function for the passivity of the system.

The step function is represented by two equilibrium assignments, at $x_0$ from time $t = 0$ to 1 s and then at $1.1x_0$ at time $t = 1$ s. The actuators dynamics are not modeled here, hence the non-physical behavior of the state $\dot{m}_{gg}$. The tuning of the controller is taken as $k_u = 0.1$, as a higher $k_u$ introduces a spike in the actuator, at time $t = 1$ s. Additionally, in order to introduce actuators dynamics in the model, which are modeled using a second order system, it is required to have a longer time-response to not outpace the actuators.

The actuators dynamics, are now introduced under the expression

$$\dot{Z}_v = -2\xi w_0 \dot{Z}_v - w_0^2 Z_v, \quad (64)$$

for $\xi = 0.7$ and the pulsation $w_0 = 25$ rad s$^{-1}$. The use of a controller with a slow $k_u = 0.1$ to maintain stability is justified by the time-response of the actuators. On Fig. 6, the actuators dynamics are added after the PBC design. The second-order behavior introduces an overshoot and oscillations that attenuate rapidly. Note that for a higher $k_u = 1$ for example, the controller outpaces the actuators and the oscillations increase in amplitude.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Step response of the controlled system for $k_u = 0.1$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Step response of the controlled system with a second order actuator for $k_u = 0.1$.</figcaption></figure>

### 6.3. Comparisons with a PI controller

A comparison is performed with a classical PI controller, to illustrate the extension of the stable zone due to the new dynamics proposed by $R_d$. On Fig. 7, the results of a simulation for a quick adjustment in the pressure is provided. The scenario used in this simulation corresponds to short-time changes in the operating point, to enhance the response performance of the controller. Additionally, the scenario overlaps the unstable zone of the polynomial functions, to illustrate the stabilizing behavior of the PBC controller. On the bottom half, the chamber pressures $p_{cc}$ and $p_{gg}$ are provided along the reference (ref) and the gas generator valve's surface $A_{vg,H}$. On the top half, the evolution of the feeding mass flows $\dot{m}_{cc,H}, \dot{m}_{gg,H}$ next to the rotational speed $\omega_H$ are figured. Fig. 7 displays the resulting trajectories of the states for both controllers. The PBC controller is able to remain stable for lower values of the states than the PI controller, e.g. $p_{ccmin,PBC} < p_{ccmin,PID}$. Note that although the pressure is lower, $p_{cc,PBC} = 0.8712$, the mass flow $\dot{m}_{cc}$ does not decrease lower than the threshold value $\dot{m}_{cc,H} = 0.86$. This is due to the changes applied to the resistive matrices, that influence the dynamics of the system.

This is explained by the fact that the PBC controller maintains a higher gas generator mass flow $\dot{m}_{gg}$ than the PI controller (see the orange and purple evolutions on Fig. 7). This allows to maintain a higher
<figure><img src="image_8.png" /><figcaption>Fig. 7. Comparison between PI and PBC responses to an unstable reference.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Real-life scenario using the PBC controller and a second-order actuator.</figcaption></figure>

rotational speed $\omega$ and higher main chamber pressure $p_{cc}$ with less variations in $\dot{m}_{cc}$. The same principle applies for low values, and result in the comparison, where the pressure $p_{cc,PBC} = 0.8712$ attains lower values than the pressure $p_{cc,PID} = 0.8999$.

## 6.4. Realistic scenario

On a realistic powered descent arc, the vehicle guidance system would react to certain kinds of variation in the navigation readings or in the mission profile by updating the commanded thrust level. In this paragraph we present a simplified reference profile alternating between a high and low thrust level with ramps connecting both levels. The scenario proposed here corresponds to the type of trajectories encountered when landing a rocket first stage. During the landing maneuver, the first stage detaches from the rocket body, and after performing a U-turn to orient the rocket engine thrust opposite to the velocity, a first burn at maximum thrust is performed to slow down the rocket engine for atmospheric re-entry. In a second time, the engine's thrust is minimized to save fuel, and allow for aerobraking. Once the velocity of the first stage is low enough thanks to friction in the atmosphere, a last burn is performed at maximum thrust to land the engine. The detailed guidance trajectory is presented in Leparoux, Hérissé, and Jean (2022).

On Fig. 8, such scenario has been applied to the PBC controlled LPRE. The main states of the engine have been represented, namely the gas generator and main combustion chamber mass flows, and the main chamber pressure. The controller is set to regulate the main combustion chamber pressure around a reference represented along $p_{cc}$ (see (63)). Fig. 8 shows that the controller is able to follow the trajectory with no overshoot for the following parameters: $k_I = 20$, $k_P = 1$, $Q_{mgg} = 1 + \frac{0.1}{R_{pc}^2}$. For a slower reference, the actuators oscillations vanish, and the trajectory tracking proves to be efficient, as the error between the reference and the real trajectory is limited.

## 7. Conclusion

Using a direct approach, proof of the stability or design of a stable controller for an LPRE system proves difficult. In this paper, its representation in a Port-Hamiltonian framework is presented that makes it possible to obtain a straightforward proof of stability of the LPRE. Moreover, this framework can be used to design a non-linear controller and a passivity storage function associated with closed-loop stability guarantees for the LPRE system. The controller allows efficient tracking of a given desired trajectory between operating points. Additionally, the tuning of the controller remains simple, using a single tuning parameter $k_u$. This controller is applied first to a simplified and then to a complex model of an LPRE, showing that the developed methodology can be extended to a large variety of rocket engines.

This stability analysis and control design can be applied to a diverse scope of complex systems for which classical Lyapunov function design proves difficult. There are particularly efficient when dealing with multi-physical systems.

As the resulting controller consists in a state-feedback function, it requires availability of measurements of the states or estimation of the unmeasured states. In the case of LPRE, the different mass flows and pressures are commonly not measured and thus require to be estimated. Future work includes the design of observers in the Port-Hamiltonian framework and further analysis of the stability properties of the controller-observer loop.

### CRediT authorship contribution statement

**Jules Gibart:** Writing – review & editing, Writing – original draft, Methodology, Investigation, Formal analysis, Conceptualization. **Hélène Piet-Lahanier:** Writing – review & editing, Methodology, Investigation. **François Farago:** Writing – review & editing, Methodology, Investigation, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

The author is an Associate Editor for Control Engineering Practice and was not involved in the editorial review or the decision to publish this article.

### Acknowledgments

This study was partly conducted under co-funding by CNES, France and ONERA, France, to whom the authors are thankful.

### References

* Aoues, S., Cardoso-Ribeiro, F. L., Matignon, D., & Alazard, D. (2017). Modeling and control of a rotating flexible spacecraft: A port-Hamiltonian approach. *IEEE Transactions on Control Systems Technology*, 27(1), 355-362.

* Breedveld, P. C. (1984). *Physical systems theory in terms of bond graphs*. (Ph.D. thesis), Technische Hogeschool Twente.

* Casiano, M. J., Hulka, J. R., & Yang, V. (2010). Liquid-propellant rocket engine throttling: A comprehensive review. *Journal of Propulsion and Power*, 26(5), 897-923.

* Cervera, J., van der Schaft, A. J., & Baños, A. (2007). Interconnection of port-Hamiltonian systems and composition of Dirac structures. *Automatica*, 43(2), 212-225.

* Cha, J., & de Oliveira, É. J. (2021). Robust guidance and control of liquid-propellant rocket engines for landing of reusable stages using fuzzy PID control. In *72th IAC IAF*.

* Dresia, K., Waxenegger-Wilfing, G., Dos Santos Hahn, R. H., Deeken, J. C., & Oschwald, M. (2021). Nonlinear control of an expander-bleed rocket engine using reinforcement learning.

* Duindam, V., Macchelli, A., Stramigioli, S., & Bruyninckx, H. (2009). *Modeling and control of complex physical systems: the port-Hamiltonian approach*. Springer Science & Business Media.
Eberard, D., & Maschke, B. (2004). Port Hamiltonian systems extended to irreversible systems: the example of the heat conduction. *IFAC Proceedings Volumes*, 37(13), 243-248.

Eberard, D., Maschke, B., & Van Der Schaft, A. (2007). An extension of Hamiltonian systems to the thermodynamic phase space: Towards a geometry of nonreversible processes. *Reports on Mathematical Physics*, 60(2), 175-198.

Esquivel-Sancho, L. M., Muñoz-Arias, M., Phillips-Brenes, H., & Pereira-Arroyo, R. (2022). A reversible hydropump-turbine system. *Applied Sciences*, 12(18), 9086.

García-Sandoval, J. P., González-Álvarez, V., & Calderón, C. (2015). Stability analysis and passivity properties for a class of chemical reactors: Internal entropy production approach. *Computers & Chemical Engineering*, 75, 184-195.

Gibart, J., Piet-Lahanier, H., Farago, F., & Galeotta, M. (2023). Regulation of a liquid propelled rocket engine using contraction theory. *IFAC-PapersOnLine*, 56(1), 307-312.

Khalil, H. (2002). Nonlinear systems. Pearson education, Prentice Hall.

Leparoux, C., Hérissé, B., & Jean, F. (2022). Optimal planetary landing with pointing and glide-slope constraints. In *2022 IEEE 61st conference on decision and control* (pp. 4357-4362). IEEE.

Liu, J., Zhang, S., Wei, J., & Haidn, O. J. (2023). RANS based numerical simulation of a GCH4/GO2 rocket engine combustion chamber with film cooling and improvement of wall heat flux prediction. *Applied Thermal Engineering*, 219, Article 119544.

Lohmiller, W. (1999). *Contraction analysis of nonlinear systems* (Ph.D. thesis), MIT.

Maschke, B. M., & van der Schaft, A. J. (1993). Port-controlled Hamiltonian systems: modelling origins and systemtheoretic properties. In *Nonlinear control systems design 1992* (pp. 359-365). Elsevier.

Ortega, R., Van Der Schaft, A., Castanos, F., & Astolfi, A. (2008). Control by interconnection and standard passivity-based control of port-Hamiltonian systems. *IEEE Transactions on Automatic Control*, 53(11), 2527-2542.

Ortega, R., Van Der Schaft, A., Maschke, B., & Escobar, G. (2002). Interconnection and damping assignment passivity-based control of port-controlled Hamiltonian systems. *Automatica*, 38(4), 585-596.

Pasumarthy, R., & van der Schaft, A. (2006). A port-Hamiltonian approach to modeling and interconnections of canal systems. In *16th international symposium on mathematical: vol. 2*, (pp. 1436-1443).

Perez Roca, S. (2020). *Model-based robust transient control of reusable liquid-propellant rocket engines* (Ph.D. thesis), Université Paris-Saclay.

Pérez-Roca, S., Langlois, N., Marzat, J., Piet-Lahanier, H., Galeotta, M., Farago, F., et al. (2018). Derivation and analysis of a state-space model for transient control of liquid-propellant rocket engines. In *9th ICMA* (pp. 58-67).

Pérez-Roca, S., Marzat, J., Flayac, É., Piet-Lahanier, H., Langlois, N., Farago, F., et al. (2019). An MPC approach to transient control of liquid-propellant rocket engines. *IFAC-PapersOnLine*, 52(12), 268-273, 21st IFAC Symposium on Automatic Control in Aerospace ACA 2019. URL https://www.sciencedirect.com/science/article/pii/S2405896319312054.

Pérez-Roca, S., Marzat, J., Piet-Lahanier, H., Langlois, N., Galeotta, M., Farago, F., et al. (2020). Model-based robust transient control of reusable liquid-propellant rocket engines. *IEEE Transactions on Aerospace and Electronic Systems*, 57(1), 129-144.

Ramirez, H., Maschke, B., & Sbarbaro, D. (2013). Irreversible port-Hamiltonian systems: A general formulation of irreversible processes with application to the CSTR. *Chemical Engineering Science*, 89, 223-234.

Rashad, R., Califano, F., van der Schaft, A. J., & Stramigioli, S. (2020). Twenty years of distributed port-Hamiltonian systems: a literature review. *IMA Journal of Mathematical Control and Information*, 37(4), 1400-1422.

Wu, Y., Zhang, H., Hui, J., & Zhou, X. (2018). Non-fragile $H_\infty$ controller for combustion process in rocket motors. *International Journal of Automation and Control*, 12(3), 381-398.

Zhou, W., Wu, Y., Hu, H., Li, Y., & Wang, Y. (2021). Port-Hamiltonian modeling and IDA-pbc control of an IPMC-actuated flexible beam. In *Actuators: vol. 10*, (no. 9), (p. 236). MDPI.