LMI-based passivation of LTI systems with application to marine structures

Funding information H2020 Marie Skłodowska-Curie Actions, Grant/Award Number: 101024372; Science Foundation Ireland, Grant/Award Number: 13/IA/1886 Abstract Due to the inherent relevance of passive (physically representative) models for control, state-estimation, and motion simulation in the field of marine systems, in this paper, an optimisation-based approach to passivation of linear time-invariant (LTI) systems is proposed with application to physically consistent dynamical modelling of marine structures. In particular, the presented strategy is based upon the introduction of a suitably designed perturbation, computed via minimisation of a linear objective subject to a specific set of linear matrix inequalities (LMIs). The performance of the passivation technique is showcased in terms of two case studies: An offshore platform, with a frequency-domain response computed by means of hydrodynamic codes, and a 1:20 scale wave energy converter (WEC), is characterised in terms of real experimental data.


INTRODUCTION
Availability of mathematical representations of physical processes, that is, dynamical systems, is a fundamental pillar for a large class of tasks/operations, including process performance assessment, system optimisation, process control, and state/fault estimation, among others. Such abstract representations of real processes must (to some extent, driven by the specific application involved) respect the underlying physical properties of the process under scrutiny, that is, they must be physically consistent. This includes, for instance, prescribed stability properties, zero dynamics, relative degree, controllability/observability, and passivity. The latter, particularly challenging to enforce, can be of paramount importance for a number of applications, specially those involving/necessitating an energy-based analysis of the corresponding process, for example, [1][2][3][4].
One specific discipline where passivity becomes fundamental is marine engineering. Physically representative models of marine structures are key tools for the development of training simulators, hardware-in-the-loop testing simulators, motion control systems, and model-based fault detection and diagnosis techniques (see, for instance, [2,5,6]). A common practice is to obtain these mathematical models based on frequency-domain data, carefully collected from the specific marine structure in question (either experimentally or numerically), by means of system identification techniques (see, e.g. [7]). In particular, an increase in interest in the use of linear time-invariant (LTI) models in recent years is evident, where such models are obtained via black-box system identification, using frequency-domain data provided by so-called hydrodynamic codes, such as [8]. This specific approach allows the computation of models using only limited information about the marine structure, that is, hull form and approximate mass distribution. Nonetheless, even if the physical problem is known to be passive, it is often the case that identified models do not effectively reflect this property, since enforcing passivity within the identification process is particularly challenging 1 . Not respecting this property can lead to a number of undesirable consequences, such as numerical instability in simulation [2], non-convexity of energy-based control procedures [4,11], and unpredictable coupling between different modes of motion [12], among others.
Motivated by this, a number of techniques have been proposed to enforce the passivity of identified models, that is, to achieve passivation, both within, and beyond, the scope of marine engineering (see, for instance, [13][14][15][16]). Aiming to keep this paper reasonably self-contained, we provide a brief discussion on the aforementioned techniques in the following. The passivity enforcement method presented in [13] considers the detection of specific ranges (in frequency) where passivity violations occur, and introduces a perturbation in terms of an augmented system, specifically designed to correct the real-part of the frequency-response mapping of the target nonpassive marine system, while preserving internal stability. Note that, since the order (dimension) of the augmented system is strictly linked to the 'number' of frequency intervals where passivity violations occur, the technique discussed in [13] introduces an additional computational cost to the final (passivised) model, which can be challenging for real-time applications (e.g. for control/state-estimation purposes). [14] and [16] assume a pole-residue description for the target model, and introduce a constrained optimisation procedure to enforce passivity by perturbing the associated residues, with a-priori knowledge of the specific location of passivity violations in the frequency-domain. Finally, [15] introduces a passivation method based upon perturbation of the Hamiltonian matrices associated with the specific target system. In particular, the strategy described in [15] begins by detecting passivity violations via the set of purely imaginary eigenvalues of the Hamiltonian, and introduces a corresponding perturbation. We note that the strategy is based on the assumption that any imaginary eigenvalue of the associated Hamiltonian matrices is either simple, or characterised by a complete sets of eigenvectors (i.e. the corresponding algebraic and geometric multiplicities associated with these eigenvalues are equal). To summarise, the vast majority of the available strategies are designed based on the detection and quantification of passivity violations (by, for instance, frequency-sweeping [17] or spectral analysis of associated Hamiltonian matrices [15]), and the introduction of a specified (commonly structured) perturbation to enforce such a property.
In this paper, in contrast, we propose an optimisation-based approach to the passivisation of LTI systems. In particular, we introduce a suitably designed perturbation, computed via minimisation of a linear objective, subject to a specific set of linear matrix inequalities (LMIs), which can be solved efficiently using state-of-the-art and readily available LMI solvers. The set of LMIs enforces the so-called Kalman-Yakubovich-Popov (KYP) lemma (sometimes referred to as the positive real lemma), by means of a perturbation on the associated state-space matrices. We note that, unlike the available methods described in the previous paragraph, the proposed strategy does not require the a-priori localised detection and quantification of passivity violations, hence being both straightforward to apply, and efficient to solve, given the linear nature of the optimisation objective. Furthermore, unlike the passivation technique for marine systems presented in [13], the strategy proposed in this paper does not require an augmented state-space description, hence being especially appealing for design and synthesis of real-time model-based control/state-estimation algorithms. Finally, and motivated by the importance of passivity in physically consistent modelling of marine structures, we showcase the performance of the proposed passivation technique in terms of two case studies: An offshore platform, characterised in the frequency-domain by hydrodynamic codes, and a 1:20th scale wave energy converter (WEC), characterised by real data collected during an experimental campaign at Aalborg University, Denmark [18].
The remainder of this paper is organised as follows. Section 1.1 introduces the notation utilised throughout our study. Section 2 recalls fundamental preliminaries of passivity for LTI systems, while Section 3 describes the proposed optimisationbased passivisation algorithm. Section 4 showcases the performance of our strategy for two different marine structures, while Section 5 encompasses the main conclusions of our study.

Notation
Standard notation is considered throughout this paper, with any exception detailed in this section. ℝ + (ℝ − ) denotes the set of non-negative (non-positive) real numbers. ℂ <0 denotes the set of complex numbers with negative real part. The symbol 0 stands for any zero element, dimensioned according to the context. The space of symmetric matrices of order n is denoted as n ⊂ ℝ n×n . If P ∈ n , the notation P > 0 (P ≥ 0) is used to say that P is positive (semi-positive) definite. Likewise, the notation P < 0 (P ≤ 0) is used to denote that P is negative (semi-negative) definite. The operator ‖ ⋅ ‖ 2 ∶ ℝ k×u → ℝ + , M ↦ ‖M ‖ 2 , denotes the spectral norm in ℝ k×u . The symbol n denotes the identity matrix in ℂ n×n . The spectrum of a matrix A ∈ ℝ n×n , that is, the set of its eigenvalues, is denoted by (A). If (A) ⊂ ℝ, then min {A} = min (A). The convolution between two functions f and g is denoted as f * g. Finally, the superscripts ⊺ and ⋆ denote the transposition, and Hermitian transposition operators, respectively.

PRELIMINARIES
This section briefly recalls fundamentals behind the concept of passivity for linear time-invariant systems. In particular, we consider a square, finite-dimensional, continuous-time, system, written, for t ∈ ℝ + , in terms of the following set of differential equations 2 , where (1) is minimal, that is, controllable and observable, x(0) = 0, and let (A) ⊂ ℂ <0 , that is, Σ is asymptotically stable. Let the complex-valued mapping W ∶ ℂ → ℂ p×p , s ↦ W (s), be the transfer function associated with (1), that is, Remark 1. From now on, for simplicity, we denote system (1) as We now recall key results with respect to the concept of passivity for LTI systems. The interested reader is referred to, for instance, [19,20], for a detailed and formal discussion on passivity (and positive realness) for a general class of systems. In particular, we recall that the asymptotically stable system ∀T ∈ ℝ + , if and only if its transfer function W is positive real (PR). Analogously, system (A, B, C , D) is strictly passive, that is, ∀T ∈ ℝ + , if and only if its transfer function is strictly positive real (SPR). The connection between passivity and positive realness, for the case of system (1), facilitates the passivity analysis of system Σ in terms of the well-known Kalman-Yakubovich-Popov lemma, also referred to as the positive real lemma (see, e.g. [21,Chapter 8]): The transfer function W of system (1) is PR if and only if there exists a matrix P ∈ n , with P > 0, such that holds. W is SPR if the inequality (5) holds strictly.

The case of (A, B, C, 0)
If the feedthrough matrix in (1) is such that D=0 (which is virtually always the case for real marine structure modelling [2,3]) the result arising from the KYP lemma requires a slight reformulation. In particular, suppose Σ D=0 = (A, B, C , 0), with (A, B, C ) as in (1), and let W D=0 ∶ ℂ → ℂ q×q , s ↦ W D=0 (s), be its associated transfer function. Then, W D=0 is PR if and only if [21,Chapter 8] there exists a matrix P ∈ n , with P > 0, such that holds. The transfer function W D=0 is SPR if the Lyapunov inequality in (6) holds strictly.

LMI-BASED PASSIVATION
As discussed throughout Section 1, the set of matrices (A, B, C , D), characterising the system Σ in (1), can be determined by a plethora of potential system identification procedures applied to experimentally/numerically generated datasets characterising the real (physical) process. Even if the physical problem is known to be passive, it is often the case that the identified model does not reflect this property, leading to physically inconsistent representations. This directly motivates the following LMI-based passivation technique, built upon the results recalled in Section 2. Suppose a passivity-check (such as those described in, for instance, [13,15]) has been applied to system (1), declaring Σ to be non-passive. Let a perturbed system Σ Δ be (compactly) defined as Σ Δ = (A, B, C + ΔC , D), with ΔC ∈ ℝ p×n . The objective posed in this paper is to design ΔC such that the perturbed system Σ Δ is passive, with ‖ΔC ‖ being as small as possible with regard to some suitable matrix norm. Before introducing the proposed passivation method, the following key remark is introduced.
We now propose an optimisation-based approach for the computation of ΔC , in terms of the following finitedimensional linear LMI-constrained problem.
Problem 1 (LMI-based passivation method). Find the perturbed model Σ Δ = (A, B, C + ΔC , D) such that ΔC is the solution of the minimisation problem min ∈ℝ + ,ΔC ∈ℝ p×n , subject to: Remark 3. Problem 1 is based on the idea of finding the smallest (in terms of the spectral norm) perturbation ΔC such that system Σ Δ = (A, B, C + ΔC , D) is passive. The latter property is achieved by enforcing the LMI conditions of the KYP lemma, that is, by forcing the transfer function of Σ Δ to be PR (see Section 2). Note that strict passivity can be directly considered via Problem 1, by enforcing the transfer function of Σ Δ to be SPR instead of PR.
Remark 4. The selection of the output matrix C for perturbation is not arbitrary: Since Σ is asymptotically stable by assumption, preserving A directly preserves its set of poles, that is, preserves internal stability. Furthermore, the feedthrough matrix D is preserved to respect the response of Σ at high frequencies. That said, note that one could perturb the input matrix B instead. Nonetheless, this would require a specific change of variables, within Problem 1, to avoid existence of bilinear terms in (8).
Remark 5. Problem 1 can be solved efficiently by means of state-of-the-art (and readily available) LMI solvers, such as those detailed in [22,23].
Remark 6. We stress that, in contrast to state-of-the-art passivation techniques, for example, [13][14][15][16] (see also the discussion provided in Section 1), Problem 1 does not require a-priori identification of the specific (frequency) location of any passivity violation, but rather directly computes the minimum norm perturbation required to ensure passivity without such knowledge, hence easing the application of the proposed technique. Furthermore, unlike the technique developed in [13] for marine systems, the presented strategy does not require the definition of an augmented system, being especially suitable for real-time control/state-estimation design and synthesis applications.

The case of (A, B, C, 0)
Suppose now that we consider system Σ D=0 , and let the corresponding perturbed system be defined as Σ Δ D=0 = (A, B, C + ΔC , 0). While one can straightforwardly define an analogous formulation to that posed in Problem 1, based upon the KYP lemma presented in Section 2.1, some readily available LMI toolboxes can have difficulties solving the equality constraint defined in (6) for the perturbed system, that is, PB = (C + ΔC ) ⊺ . To overcome this potential implementation issue, this equality constraint is approximated following [24,Chapter 3], that is, in terms of the minimisation of its associated spectral norm: for any matrix P ∈ n . We are now ready to introduce the analog to Problem 1 for systems with zero feedthrough.
Remark 7. Analogously to Problem 1, the methodology presented in Problem 2 is based on the idea of finding the smallest (in terms of the spectral norm) perturbation ΔC such that system Σ Δ D=0 = (A, B, C + ΔC , 0) is passive. The latter property is achieved by enforcing the conditions of the KYP lemma for systems with zero feedthrough (see Section 2.1), that is, by forcing the transfer function of Σ Δ D=0 to be PR. Note that strict passivity can be achieved by tightening the latter condition to be SPR.
Remark 8. In contrast to Problem 1, the minimisation objective is now defined as the sum of two different contributions, that is, + , where refers to the minimisation of ‖ΔC ‖ 2 , and accounts for the equality constraint as defined in (9). The 'weights', w and w , are added to Problem (10) to assist in trading-off these two objectives: w > w 'prioritises' minimisation of ‖ΔC ‖ over satisfaction of (9), and vice versa. These two weights can be determined by the user via numerical experience, depending on the specific application involved.

Algorithmic overview
We summarise the strategy, presented throughout Section 3, in the following paragraphs. In particular, an algorithmic overview of the proposed passivation method is presented in Figure 1, including a detailed description of each step involved. As can be appreciated from Figure 1, we start with either a set of numerical, or experimental, data characterising the dynamics of the system under scrutiny. Note that the former is commonly computed via hydrodynamic codes for the case of marine structures (see also the discussion provided in Section 1). With the definition of such a set, standard frequency-domain system identification techniques [7] are applied to compute a statespace description (A, B, C , D) describing the dynamics of the system. Remark 9. If experimental data is considered for the characterisation of the system, post-processing techniques are commonly employed before performing system identification, with FIGURE 1 Algorithmic overview of the proposed passivation approach the objective of, for instance, filtering any undesired effects, such as measurement noise (see [7,25]).
Though highly unlikely, if the identified model (A, B, C , D) is determined to be effectively passive, no passivation technique is required, and the algorithm naturally stops. If, in contrast, (A, B, C , D) is non-passive (which is normally the case), we introduce a perturbation ΔC in the output matrix of the target non-passive system, that is, we build a perturbed model (A, B, C + ΔC , D). If the target non-passive system is biproper, that is, D ≠ 0, we solve Problem 1 to compute the minimum norm ΔC to produce a passive system (A, B, C + ΔC , D), and FIGURE 2 Schematic of a generic marine system, including a description of the corresponding set of input-output variables the algorithm stops. However, if the target non-passive system is strictly proper, that is, D = 0, Problem 2 is considered analogously.

PASSIVATION OF MARINE STRUCTURES
Marine structures with zero forward speed can be generally represented in terms of the feedback interconnection between two different LTI subsystems, as depicted in Figure 2. We provide a brief description of the nature of these systems in the following paragraphs. The interested reader is referred to, for instance, [5], for a thorough treatment of this topic.
In particular, considering so-called linear potential flow theory (see [5]), the equation of motion of a marine structure with zero forward speed can be generally expressed in terms of the following linear continuous-time system: where z ∶ ℝ + → ℝ N , t ↦ z(t ), denotes the displacement vector of the marine system, with N ∈ ℕ the corresponding number of degrees-of-freedom (i.e. modes of motion) associated with the structure, f rad ∶ ℝ + → ℝ N , t ↦ f rad (t ), the so-called radiation force, which accounts for the fluid memory effects acting on the structure, , the viscous force, and f re ∶ ℝ + → ℝ N , t ↦ f re (t ), the restoring force. The input mapping u ∶ ℝ + → ℝ N , t ↦ u(t ), represents the sum of the external forces acting on the marine structure, including uncontrollable effects (i.e. the so-called wave excitation force), and controllable signals (i.e. user-designed control forces), while  ∈ ℝ N ×N is the so-called generalised WEC mass matrix [5]. Finally, note that the output mapping y ∶ ℝ + → ℝ N , t ↦ y(t ), is given by the velocity vector characterising the marine structure under analysis.
Remark 10. We set the output of system (11) to be the velocity of the structure due to the underlying importance of this variable in the design and synthesis of several families of controllers (see [3,26]). Nonetheless, note that we do this without any loss of generality, since either displacement, or any other linear combination of the marine structure motion variables, can also be chosen as the output of system .
Following linear potential flow theory, the hydrostatic restoring force f re is expressed in terms of a linear relation involving the system displacement, that is where s h ∈ ℝ N ×N is the so-called hydrostatic stiffness. Similarly, viscous effects are represented in terms of a relation proportional to the velocity of the marine system, that is with s v ∈ ℝ N ×N a viscous coefficient arising from a variety of different linearisation procedures (see, e.g. [27,28]). Finally, the radiation force is characterised in terms of the so-called Cummins' equation (see [29]), by means of the following convolution mapping: where the mapping k rad represents the radiation impulse response function.
Since the impulse response k rad fully characterises an LTI system, one can decompose system  in (11) in terms of the feedback interconnection of two main subsystems, S 1 and S 2 , that is with u S 1 = u − y S 2 , and Note that system S 1 represents inertial effects, while system S 2 represents fluid memory effects that incorporate the energy dissipation due to waves radiated as a consequence of the motion of the structure, that is, radiation forces.
Remark 11. Both systems S 1 and S 2 should be passive (see [2,9]). Note that this inherently guarantees that the closed-loop (inputoutput) behaviour  in (11) is passive and, hence, is automatically internally stable [20].
While system S 1 depends upon a finite number of fixed and static system parameters (e.g. mass and hydrostatic stifness of the structure), S 2 is built upon the specific impulse response FIGURE 3 (a) 3D-render and (b) low-order mesh (utilised in NEMOH to compute the frequency-domain characterisation of the associated radiation system S 2 ) for the analysed offshore platform function k rad , which is virtually always computed by means of hydrodynamic codes. Since these codes provide a frequencydomain characterisation of such a radiation system, a parametric (state-space) form for system S 2 , that is, a structure more suitable for time-domain motion simulation/control/estimation, is commonly obtained in terms of black-box system identification procedures (see Section 1). This, together with the discussion provided in Remark 11, highlights the importance of having suitable passivation techniques for S 2 , hence guaranteeing physically representative input-output dynamical models for any marine structure. This is specifically illustrated in Section 4.1, in terms of an offshore platform.
Another potential path is to apply a system identification procedure to the input-output system  directly, that is, using the frequency-response mapping associated with the closedloop behaviour. While one can guarantee input-output stability relatively straightforwardly with this approach, passivation is virtually always still required to have a physically consistent model. This is specifically considered in Section 4.2, where the response of a wave energy converter, computed directly from input-output experimental data, is forced to be passive by means of the presented technique.

Passivation of S 2 : An offshore platform
We consider model passivation for the radiation subsystem S 2 of an offshore platform, as schematically illustrated in Figure 3a. In particular, the open-source hydrodynamic solver NEMOH [8] is used to compute the non-parametric frequency-response associated with system S 2 , considering that the platform is allowed to move in two different degrees-of-freedom (DoFs): surge and heave. The low-order mesh, utilised by NEMOH to characterise the frequency-response of the platform, is illustrated in Figure 3b. Based on the frequency-domain results computed with NEMOH, an internally stable, strictly proper (i.e. with D = 0), non-passive, radiation system S 2 is computed using momentmatching [3,30], defined over a state-space of dimension (order) n = 26. Passivation is then applied to this system, by computing an optimal output perturbation ΔC , as detailed throughout Section 3. Figure 4 presents the sigma plot (i.e. the singular values associated with the frequency-response mapping) for the nonpassive identified model (dashed-black), and the model arising from subsequently applying the passivation technique presented in this note (solid-green). Since, as can be appreciated from Figure 4, these models are virtually identical in terms of their FIGURE 4 Sigma plot for the non-passive identified model (dashed-black) for the analysed offshore platform, and the model arising from subsequently applying the passivation technique presented in this note (solid-green)

FIGURE 5
Passivity indicator  for the non-passive (dashed-black), and passive (green-solid) models, for the analysed offshore platform. The passivity boundary is indicated with red color associated gains, we define the following real-valued operator to highlight any passivity violations: Let  ∶ ℝ → ℝ, ↦  ( ), be defined as where, clearly, if  ( ) < 0, for any ∈ ℝ, then the transfer function W is not PR, and hence its associated system Σ is non-passive. Figure 5 presents the indicator (17) for the non-passive system obtained using system identification (dashed-black), and its corresponding passive counterpart (solid-green), computed via the passivation technique proposed Set of eigenvalues associated with the closed-system , for both the non-passive model S 2 (black), and its passive counterpart (green), computed by means of the presented passivation strategy in this paper. Clearly, the former presents negative values for , hence directly highlighting its non-passivity. It is noteworthy that ‖ΔC ‖ 2 ∕‖C ‖ 2 = 0.0036 for this case, that is, only a very small (in terms of the spectral norm) perturbation is required to achieve passivation with the presented strategy.
Finally, and to highlight the importance of guaranteeing passivity in the parameterisation of S 2 for marine structures, the set of eigenvalues associated with the closed-system  is presented in Figure 6, for both the non-passive model (black), and its passive counterpart (green). Note that, even though a there is an almost negligible difference in the frequency-domain response of the non-passive system S 2 , and its passive counterpart (computed by means of the presented method), the input-output response arising from feedback interconnection with S 2 renders, for this case, an internally unstable system, hence losing any practical value in terms of analysis/prediction of the platform motion.

Passivation of : A wave energy converter
We consider, in this section, model passivation for the so-called Wavestar WEC [31]. This device consists of a hemispherical hull with a single operational DoF in pitch. On the full scale WEC,  Figure adapted from [18]. The reader is referred to [18] for further details the hydraulic power take-off (PTO) system consists of a cylinder, pumping fluid through a generator, with a rated power of 500 kW, for a device with 20 floaters [32]. Here, a single 1:20th scale model of the full scale device, as shown in Figure 7, is considered, with an electrical, direct drive, actuator PTO, fully described in [18].
In particular, following the procedure described in [18], a model for this device is obtained via black-box system identification, using input-output experimental data arising from chirp experiments performed on the physical system, which is located in a wave basin in Aalborg University, Denmark 3 . In particular, subspace-based techniques [25] are utilised for identification, rendering the internally stable, strictly proper, nonpassive, input-output, system , defined over a state-space of dimension (order) n = 6, described in terms of the set of matrices (A, B, C , 0) in (18). We apply the passivation method presented in this paper to calculate a perturbation ΔC as detailed in (18). Note that ‖ΔC ‖ 2 ∕‖C ‖ 2 = 0.0369, so that the perturbation computed to passivise system (A, B, C , 0) is, as in the case of Section 4.1, very small in terms of spectral norm.
The Bode plot associated with the input-output response for the corresponding WEC system is shown in Figure 8, including the set of experimental empirical transfer function estimates (solid-grey tones), the non-passive model identified via subspace methods (dashed-black), and the passive model computed via subsequently applying the passivation technique proposed in this note (solid-green). Furthermore, the Nyquist plot associated with non-passive (dashed-black) and passive (solid-green) systems is presented in Figure 9, clearly illustrating the passivity violation of the former, which has a Nyquist plot taking values beyond the positive real semiplane.

CONCLUSION
We present, in this paper, an LMI-based passivation technique for LTI systems, with specific application to physically representative mathematical modelling of marine structures. The proposed technique is based upon finding a small perturbation,

FIGURE 8
Bode plot for the input-output system  for the experimental WEC. Solid lines in grey tones indicate empirical transfer function estimates, for different combination of input-output signals. The identified state-space model is denoted with dashed-black line, while its passive counterpart, obtained via the proposed passivation method, is denoted with solid-green

FIGURE 9
Nyquist plot for the input-output system  for the experimental WEC. The identified state-space model is denoted with dashed-black line, while its passive counterpart is denoted with solid-green written in terms of the output matrix of a given state-space system, such that the conditions of the KYP lemma hold, that is, the transfer function associated with the perturbed system is PR (or SPR, if strict passivity is to be imposed). The performance of the strategy is illustrated in terms of two different marine structures with zero forward speed, highlighting the importance of such a technique in physically consistent modelling in marine applications.