{'fraction_non_alphanumeric': 0.0294382257558322, 'fraction_numerical': 0.009980123443874882, 'mean_word_length': 4.469275660830759, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlinear state-space modelling of the kinematics of an oscillating circular cylinder in a fluid flow

The flow-induced vibration of bluff bodies is an important problem of many marine, civil, or mechanical engineers. In the design phase of such structures, it is vital to obtain good predictions of the fluid forces acting on the structure. Current methods rely on computational fluid dynamic simulations (CFD), with a too high computational cost to be effectively used in the design phase or for control applications. Alternative methods use heuristic mathematical models of the fluid forces, but these lack the accuracy (they often assume the system to be linear) or flexibility to be useful over a wide operating range. In this work we show that it is possible to build an accurate, flexible and low-computational-cost mathematical model using nonlinear system identification techniques. This model is data driven: it is trained over a user-defined region of interest using data obtained from experiments or simulations, or both. Here we use a Van der Pol oscillator as well as CFD simulations of an oscillating circular cylinder to generate the training data. Then a discrete-time polynomial nonlinear state-space model is fit to the data. This model relates the oscillation of the cylinder to the force that the fluid exerts on the cylinder. The model is finally validated over a wide range of oscillation frequencies and amplitudes, both inside and outside the so-called lock-in region. We show that forces simulated by the model are in good agreement with the data obtained from CFD.


Introduction
The kinematics of oscillating bluff bodies in a fluid flow have been a hot research topic for decades. Most important are the potentially harmful vortex-induced vibrations (VIV) where the structure is excited by alternating vortex shedding. Classical examples are the vibration of chimney stacks exposed to wind, pipelines on the sea-bed excited by the ocean currents or tubes in heat exchangers [1].
For this kind of fluid-structure interactions, typically an analytical solution cannot be found. However, accurate predictions of the kinematics and the resulting dynamics are vital for instance in the design process, in monitoring applications or for control. To obtain predictions with high fidelity, so far only two possibilities existed: solving the Navier-Stokes equation via computation fluid dynamic (CFD) simulations or performing experiments. Both approaches are time-consuming, expensive, and require dedicated lab or computing facilities. These drawbacks make the current methods too expensive for many intended applications, where only limited time and resources might be available.
To tackle the shortcomings of CFD and experiments, low-order models of VIV have been developed [2]. Most models are heuristic, with only a limited range of applicability. Different types exist, the most noteworthy can be classified as force decomposition models [3], single degree-of-freedom (sdof) models [4] or wake oscillator models, where the wake oscillation is generally formulated as a Van der Pol equation [5,6,7].
Although capable of qualitatively describing the characteristic behaviour of VIV, an efficient and powerful model, flexible enough to span a wide domain in parameter space with a single set of parameters has yet to be obtained. To identify such a model, system identification can be a very powerful tool. In recent years, the application of linear system identification techniques has already shown some promising results, e.g. the best linear approximation in least squares sense [8] or building an input-output relationship using delayed values of the input and the output (ARX model) [9].
Nevertheless, the challenges in modelling the system at hand are substantial. Fluid-structure interactions are, for one, inherently nonlinear [10]. This is strongly pronounced by the fact that VIV is a selfexcited yet self-limited oscillation, resulting in a stable limit cycle [11]. In addition, it has been shown that the vortex shedding behaviour is hysteretic [12]. Nonlinear modelling of VIV is currently a very active research topic. In [13], the auto-regressive moving averaging (ARMA) technique is combined with the modal analysis method. A nonlinear force representation is obtained by including higher harmonic terms resulting from the fluid-structure coupling. The coefficients in the model are determined to provide a best fit to the measured time series in maximum-likelihood sense. Reasonable predictions in terms of maximum amplitude are found, although only a limited amount of validation cases were investigated. In recent work [14], local linear models are used as a surrogate model for aerodynamic loads on wings. A general framework is however not yet in place when it comes to nonlinear identification of complex dynamical systems. A lot can still be gained in terms of quality and applicability of the model but also in terms of insight into the nonlinear behaviour of VIV.
In this work, we propose a novel approach to model flow problems that applies state of the art nonlinear identification techniques. The powerful framework of nonlinear state-space models is explored and put to the test. Nonlinear state-space models have proven to be very rich and flexible, providing solutions to a variety of nonlinear problems also including hysteresis [15,16].
The methodology of this work can be summarised as follows: 1. Use established methods, CFD simulations in this case, although the proposed method is equally suitable for experiments, to study the fluid forces within a user defined domain of parameter space. Then capture these observations under the form of time series. 2. Identify a nonlinear black-box model that is able to reproduce the observed phenomenon. 3. Exploit the benefits of the obtained model to simulate, at very low cost, a variety of possible regimes of the system.
The lay-out of the paper is as follows: in Section 2, the simulation set-up is introduced. Section 3 describes the discrete polynomial nonlinear state-space model structure and how such a model is identified. Section 4 consists of two parts. In the first part, the method is applied to a nonlinear oscillator of the Van der Pol type. In the second part the kinematics of an oscillating cylinder is addressed. Section 5 concludes the paper.

Simulation set-up
VIV are often studied on a 2D circular cylinder, suspended by a spring and a damper and restricted to move in the transverse direction, perpendicular to the incoming flow. Such a set-up is a self-excited oscillator, where the fluid forces generate the displacement. Common practice, however, is the use of imposed (or forced) oscillations [12,17,18] to study the interaction between the flow field and the cylinder. With total control of the displacement one can explore a wide range of frequencies and amplitudes, contrary to the limited response behaviour observed when the mass, stiffness and the damping are fixed to certain values. The set-up of both the free and the imposed oscillations can be seen in Fig. 1.

Input-output description of the system
Since the displacement is imposed, this is logically considered the input to the system. To quantify the behaviour of the wake one could consider any of the classical fluid variables: velocity, vorticity or pressure in some points downstream. Although not all variables are equally informative [19], they remain good candidates to build a model from. More obvious, from an application point of view, is a displacement-force relationship. Moreover, the fluid force on the cylinder can be seen as a combined fluid variable, concentrated in a single observable point, the centre of the cylinder. The system can thus be considered as a single input-single output system. To generalise the findings, both the input and the output are made dimensionless. The input is given by A/D with A the amplitude of the displacement and D the diameter of the cylinder. The output is given by the force coefficient c y , where F y is the resultant of the forces acting on the cylinder in the y-direction, perpendicular onto the incoming flow, U is the unperturbed flow speed, ρ is the density of the fluid.

The autonomous oscillator c y
The frequency of natural vortex shedding, f St , is given as Strouhal's relationship: where D is the diameter and U is the fluid velocity of the unperturbed flow [20]. S is the so-called Strouhal number, which is a slowly varying function of the Reynolds number (Re) and has a value of about 0.16 around Re = 100 [21]. Even in the case of a static cylinder, c y will oscillate at f St . In other words, c y behaves as an autonomous oscillator, tracking a stable limit cycle [11]. If however an excitation (imposed oscillation) of certain frequency and amplitude is applied, the frequency of vortex shedding can deviate from the natural vortex shedding frequency and synchronise to the excitation frequency. The region of parameter space for which this synchronisation phenomenon is observed is called the lock-in region [21]. The lock-in region is of special interest since within this region the maximum values of the fluid forces are observed [1].

CFD simulations of an oscillating circular cylinder
Using CFD simulations of the flow around cylinders with imposed oscillations, time series of the fluid variables are obtained. Broadband signals that span a frequency band around f St are chosen as excitation signals. The simulations are performed using the open source CFD package OpenFOAM [22]. A transient solver called pimpleDyMFoam with a variable time step was used to solve the Reynolds averaged Navier Stokes (RANS) equations for an incompressible flow in a finite volume discretisation scheme [23,24]. All simulations were performed at Re = 100 to ensure laminar, predominantly 2D [1], vortex formation. The unperturbed velocity U is set by the Strouhal relationship (Eq. 2), selecting f St = 3 Hz results in U = 1.7964 m/s for S = 0.167 (this is the uncorrected S , see below). The cylinder diameter D = 0.10 m. Consequently the kinematic viscosity ν = 1.7964 × 10 −3 m 2 s . Since the flow is laminar, no turbulence model was used. The time step was adjusted by the solver, satisfying at all times the Courant condition, Cr < 0.5 with Cr the Courant number. After transients have settled the time step generally converged to a value O(2 × 10 −5 ) s. From these raw data, time series are extracted using linear interpolation at a fixed sampling frequency of f s = 50 Hz. Given that the highest excited frequency is 4.5 Hz (see Section 4.2.2), output signals containing up to the fifth harmonic will be alias-free.
Since the vortex shedding is only 2 dimensional at Re = 100, a 2D mesh with the cylinder positioned 10 D from the inlet and centred in the y-direction was used. The computation domain is 40D long (L) and 30D high (H). The top and bottom of the domain are constrained by a slip boundary condition while on the cylinder surface itself, the velocity is set equal to the grid velocity, which follows from the cylinder motion (no-slip condition). At the inlet, a uniform velocity profile is used. The outlet is conditioned with a zero velocity gradient.
Validation of the CFD models was carried out based on vortex shedding related quantities, which are tabulated in Table 1. Two meshes were evaluated, a fine mesh containing a total of 55362 hexahedra cells and a coarser one of 35172 hexahedra cells. The tested quantities showed to be independent from the mesh used with deviations well below 1 %. All the data that were used for identification was produced using the finest mesh, although the coarser one would have sufficed given the accuracy when compared to literature.
For the aspect ratio of our computational domain (H/D = 30), we can expect a small effect of blockage on the force coefficients and Strouhal number. [25,26]. To compensate for this blockage effect, the free stream velocity was correct from the formula derived in [27] and quoted in [28,29,30], where U and U are the real and the corrected free stream velocity correspondingly, D is the diameter of the cylinder, H is the height if the computational domain and c d is the uncorrected mean drag coefficient. In [29] it was suggested that c d can be corrected from: with c d the uncorrected drag coefficient and c d the corrected one. It was however noted that the discrepancy would be very small if instead of utilising Eq. 4, the drag coefficient would be calculated using the corrected free stream velocity U . Using U also c y and the Strouhal number can be corrected. Note that a correction in U also implies a correction in Re. The true Reynolds number, given the present blockage ratio is approximately equal to 101. Fig. 2 illustrates the vortex formation in the wake behind an oscillating cylinder.  [26] 0.225 1.319 0.1648 Table 1: Validation results on stationary cylinder simulations at Re= 100 where c y denotes the rms value of the corrected force coefficient in the y-direction (perpendicular on the oncoming flow), c d is the mean value of the corrected drag coefficient and S is the corrected Strouhal number in the case of the present study and the reported values from literature found in [26].

Black-box identification using a polynomial nonlinear state-space model
System identification generally requires a four-step approach. First of all, it is essential to collect data of the user-defined problem (Section 3.1). To that end, an experiment, tailored to the problem at hand, is designed. Secondly a model structure is selected (Section 3.2). The model structure has to be flexible or generic enough to be able to capture the observed behaviour of the system. The overall goal will be to match the model as well as possible to the data. A part of this third step will be to evaluate the fit of the model. Therefore the notion of cost is introduced (Section 3.3). To obtain a model estimate, the objective will be to minimise the cost, given the data and the selected model structure (Section 3.4). Finally, the model is validated on a new dataset, not used during the identification, to avoid an over-dependence on stochastic contributions (Section 3.6). This step will also require estimates of the initial states and input (Section 3.5).

Experiment design
The fact that it is unlikely that the true model is found implies that the identified model becomes dependent on the provided data. It is therefore crucial that the coefficients of the model are tuned on data which span the operating domain as well as possible, in this way the errors are pushed into regions that are of little or no interest. When selecting the excitation signals of the training data, not only the frequency band and the amplitude, but also the type of signal must be taken into account. The intended operating conditions, richness and compactness of the data are all important factors. The selection of an adequate training set for identification of c y is discussed in Section 4.2.2.

Model structure
To be able to work with a data-based black-box modelling approach, the model structure has to be flexible enough to cover the problem of interest. The Polynomial Nonlinear State-Space (PNLSS) structure [37] is able to provide such flexibility. An additional advantage of the state-space representation is that it is very well suited to incorporate multiple-inputs and multiple-outputs (MIMO). For now, only single-input single-output (SISO) of a 2D structure is modelled. In future work it might be considered to exploit the MIMO abilities of the structure to relate the displacement of multiple spanwise cross sections of a 3D structure to their corresponding force coefficients. It should be noted that pressure fluctuations at certain cross sections may nonlinearly effect others. In recent work [38,39] this aspect has already been addressed using the aerodynamic strip theory. MIMO nonlinear state-space models could form an interesting alternative.
The nonlinear state-space model structure is defined as a classical linear state-space model where the state and the output equation are extended with some nonlinear functions f and g Eq. (5). The model is built in discrete time.
To maintain a nomenclature suited to the model to be built, c y is used as output variable and y is the input, which corresponds to the displacement of the cylinder in the y-direction.
When the linear part is of n x -th order (n x state variables), the state vector x(t) ∈ R n x , the input vector y(t) ∈ R n y with n y equal to the number of inputs and c y (t) ∈ R n Cy , which is the vector of all the n Cy outputs. From hereon only the SISO case is considered where n c y = n y = 1. Correspondingly, the dimensions for the matrices of coefficients become: In the PNLSS case, the nonlinear functions f and g are expanded in basis functions, leading to Eq. (6): where vectors ζ(t) ∈ R n ζ and η(t) ∈ R n η contain nonlinear monomials in x(t) and y(t). The monomials range from degree 0 up to a chosen degree p for the state equation Eq. (6a) and q for the output equation Eq. (6b). All degrees within this range are included, except for degree 1, since these contributions are already covered by the linear part 6 of the structure. The monomials ζ(t) and η(t) are formed by all possible products of the input and the state variables raised to a chosen degree, e.g. for an element in ζ(t): with the total degree of the monomial satisfying the condition: . Therefore E ∈ R n x ×n ζ and F ∈ R n x ×n η , with n ζ and n η the number of combinatorial products between the states and the input, satisfying the condition on the total degree of the monomial up and till p for E and q for F.

The cost function
In this work, the cost is defined as the squared error between the output of the model and the training data c y . A time-domain cost function was preferred to facilitate the use of signals of different types and with different frequency content. Since the cost is dominated by model errors and not by noise, no noise weighting is applied. The least-squares cost function is given by: where θ is the vector of parameters to be optimised. These parameters are all the coefficients in the A,B,C,D,E and F matrices: here vec denotes the operator which stacks all the elements of a matrix in a single column vector. c y,CFD denotes the training data while c y,mod is the simulated output of the model.

Minimising the cost function
The procedure of identifying a PNLSS model consists of 3 major steps [37]: 1. First, the best nonparametric linear approximation (FRF) of the system is estimated. The best linear approximation, denotedĜ BLA , of a nonlinear system is defined as the model G, belonging to the set of all linear models G such thatĜ where C Y ( jω) and Y( jω), denote the discrete Fourier transforms of c y , and y respectively. 2. Next, linear identification techniques are used to obtain a linear parametric estimate of the model,Ĝ BLA,par (θ lin ) (see Eq. 16 for a definition of θ lin ). 3. Finally, all coefficients are further tuned to data via nonlinear optimisation of the cost function Eq. (8). As initial values for the A,B,C and D matrices the parametric BLA is used, E and F are set to zero.

Initialisation of the nonlinear model via the best linear approximation
The best linear approximation, as defined in Eq. (10), is used to initialise the full nonlinear model. To obtain the best linear approximation,Ĝ BLA , random-phase multisine realisations are used [40]. The input signal, y(t), corresponds to a sum of N harmonically related sines with a user-defined amplitude spectrum. The phases φ n are random variables from a uniform distribution between 0 and 2π.
In general, M different random phase multisines are applied and P number of periods of the input and the output signal are measured. Per experiment m and period p the discrete Fourier transform (DFT) of the input signal (Y(k) [m,p] ∈ C) 7 and of the output signal (C Y (k) [m,p] ∈ C) are computed. The spectra are then averaged over the periods. At frequency line k we get:Ŷ For every experiment m, the transfer function estimate,Ĝ(z k ) [m] ∈ C, with z k = e j 2π fsk N , f s the sampling frequency and N the number of frequency lines, then becomes:Ĝ Finally, the transfer function estimates are averaged over the realisations to get the Best Linear Approximation: An advantage of using such periodic excitations is that the influence of nonlinear distortions can be separated from measurement noise [41]. Variations found over the P periods of the measured transfer functions are caused by measurement noise while variations over the transfer functions of M experiments (phase realisations) are due to the combined effect of measurement noise and nonlinear distortions (see Section 4.2).
From this nonparametric estimate and an estimate of the sample variance (σ 2 BLA ), a parametric linear state-space model,Ĝ BLA,par , is estimated using the frequency domain identification toolbox (FDIDENT 1 ) in MATLAB. Here the nonlinear cost function is minimised in the frequency domain. The cost function is weighted by the total variance which contains contributions of both the noise and the nonlinear distortions: where the vector of all the linear coefficients is given by:

Nonlinear optimisation of the full model
The full nonlinear model is obtained by minimising the least-squares cost function Eq. (8) with respect to θ. The optimisation is performed using the Levenberg-Marquardt (LM) algorithm. It is well known that gradient descent based methods, such as LM, are very sensitive to the choice of initial values. Given the starting values, the optimisation can be seen as a descending path along the cost function. Unfortunately, the optimisation might easily end up in a poor local minimum, close to the starting point. To increase the chances of ending up in a sufficiently good local minimum, one would like to descend the cost function over multiple paths, in different descending directions. A way to come to different descending paths, from a certain starting point, is by calculating the gradient based on a subset of the parameters, as opposed to all parameters simultaneously. Keeping a number of parameters constant during optimisation reduces the degree of freedom, resulting in a different path.
An example of a fictive, second degree of freedom cost function, is given in Fig. 3. By optimising both parameters, θ 1 and θ 2 , simultaneously, the red path is followed directly to the local minimum L 1 . Alternatingly optimising both but starting from θ 1 , keeping θ 2 constant to its initial value, results in the path indicated in black. Again ending up in L 1 . Changing only the order of optimisation, now starting from θ 2 , creates the green path, which leads to the more optimal minimum L 2 . Depending on what parameters and the order in which the parameters are selected for optimisation, a different path along the cost function is followed. If both θ 1 and θ 2 are optimised simultaneously, in a gradient descent based way, the red path is followed, ending up in the poor local minimum L 1 . If θ 1 and θ 2 are optimised alternatingly, starting by θ 1 , the optimisation would again end up in L 1 (black path). Starting by θ 2 , however, L 1 could be avoided, ending up in the more optimal minimum L 2 (green path).
In this work, the same strategy is applied to the optimisation of the PNLSS model. In that case, as will be explained in Section 4.2, a total of 2772 parameters need to be optimised. Since the true optimal subset of parameters leading to a good minimum is unknown, a number of predefined subsets are evaluated. A typical subset would include all the linear coefficients in the A, B, C and D matrices and only a selection of the coefficients in E and F. In total, we chose to evaluate 90 possible paths from the starting point. The possibilities are explained in detail in Appendix A. The one leading to the most optimal point is selected and the search is repeated. When none of the 90 options yields a better minimum, the optimisation is ended. A list of the tested subsets along with a schematic overview of the optimisation loop is given in Appendix A, Fig. A.1.

Estimating the initial state of the model
Given the fact that discrete PNLSS models are formulated recursively, x(t + 1) = f (x(t), u(t)), the initial state x 0 and initial input u 0 need to be known at t = 0, in order to be able to simulate the response of the model to a certain input sequence. To obtain values for x 0 and u 0 , they can be seen as parameters that need to be estimated from the data. The problem is again nonlinear in the parameters and the Levenberg-Marquardt algorithm is used. As initialisation x 0 and u 0 are set to zero. At this stage, no guarantee of finding the true initial state and input can be given and errors can be introduced.
When dealing with nonlinear systems, the impact of an error on x 0 and u 0 can be significant. This is due to the fact that there might exist multiple stable solutions (nonunicity) for a set of different initial conditions. Although a full sensitivity analysis revealing the bifurcation structure is out of the scope of this work, the impact of the initial state on the fluid force coefficient will be illustrated with an example in Section 4.2.4.

Validating the model
How to assess the quality of the model will strongly depend upon the user. Generally, there exists a trade off between the accuracy of the model and the extent of parameter space covered. In this work we will demonstrate that a single model, valid over a wide range of parameter space, can be obtained at the expense of accuracy in some isolated regions. See Section 4.2.5 for an in depth discussion on the validation.

Application of the method
In this section, the PNLSS identification procedure is applied to two systems: the nonlinear Van der Pol oscillator, and the problem of modelling c y . The Van der Pol oscillator is used as an illustration of the method since it can be shown that there exists an analytical PNLSS equivalent.

The Van der Pol oscillator
The Van der Pol oscillator is of special interest since such equations have been frequently used as wake oscillator models. One of their characteristics is that they are able to perform autonomous oscillations, and also describe a synchronisation behaviour when they are excited [5,6]. The most basic representation is given by [42]: where on the left hand side we have a nonlinear oscillator with an autonomous angular frequency Ω aut and on the right hand side we have a forcing term related to the derivative of the displacement (forced Van der Pol equation).
There is debate on what the best right hand side expression should be used to model VIV, displacement, velocity, or acceleration [7]. In this article we use the velocity coupling.
Here Ω aut = 2π f aut with f aut the autonomous frequency. We chose µ equal to 0.3, in accordance with [7], to assure the characteristic behaviour of frequency synchronisation, limit cycle oscillations and autonomous oscillations are present.

Visualisation of the lock-in range
We study the lock-in region by imposing a sinusoidal displacement, y(t) = A sin(2π f ex t), of various amplitude levels and ranging in frequency around f aut . In Fig. 4, the combinations of amplitude and relative frequency of y, for which the frequency of c y locks on to the excitation frequency ( f ex ) are indicated with a black marker. A distinct vshaped region, also called Arnold tongue [11], becomes visible. More information on the response and the bifurcation structure of the forced Van der Pol equation can be found in [43].

An analytical PNLSS expression of the Van der Pol equation
It can easily be shown that the Van der Pol equation is covered by the PNLSS model class. Using a simple Euler discretisation, Eq. (17) can be rewritten in a discrete PNLSS form. With the state variables being x(t) = [c y (t)ċ y (t)] T and T s indicating the time step the equivalent state-space formulation becomes: Since in the case of an oscillating cylinder in a fluid flow, the analytical expressions are unknown, we will also tackle the identification of the Van der Pol oscillator in a black-box way. We will show that a good estimate of Eq. (18) can be found working solely from data, and thus without making use of any of the coefficients from the analytical expression, nor from the type (second order) of nonlinearity present. To get an estimate of the PNLSS model, steps 1 to 3 as described in Section 3.4 are executed. Using the FDIDENT toolbox in MATLAB, a parametric model scan was conducted. In Fig. 6 the cost of the parametric fit is plotted against the model order. In this case a linear second-order model (1/2) was found adequate (clear drop in cost). Moreover it can be seen from the analytical expression in Eq. 18 that 2 state variables would suffice. A plot of the nonparametric and the parametric estimate ofĜ BLA can be seen in Fig. 7. Also the level of total distortion and measurement noise are indicated in red and green correspondingly. Even though the data was not intentionally corrupted with noise, still a certain amount of distortions were measured. One major reason is that the autonomous frequency was not exactly equal to 3 Hz, because of integration errors when solving Eq. 17 in a discrete way. As a result leakage could not be avoided.

Nonlinear optimisation
Starting fromĜ BLA as initial values for the A, B, C and D matrices, the nonlinear optimisation strategy as explained in Section 3.4.2 is applied to minimise the cost function: where c y,VdP is the data computed from Eq. (17) and c y,mod is the simulated output of the model. The 7 multisine realisations are used as training data. Since the objective is to model the Van der Pol equation in a black-box way, no precise selection of the nonlinear degrees was made. Instead, for both the state and the output equation all monomials ranging in degree from 2 up and till 5 were used (whereas judging from Eq. 17 a total degree 3 would be exact). In the output equation also the DC term (degree 0) was included. Theoretically, it should be possible to retrieve the exact model parameters (degree 3) by introducing the proper constraints on the A, B, C and D matrices. In practice, however, this is not straightforward, as there is no intuitive way to initialise the remaining coefficients in Eq. 18.
After the linear identification step, the model showed an rms error of 33.1 % when simulating the output of the  We have thus successfully retrieved a discrete PNLSS model that well describes the data obtained from Eq. (17) without imposing any prior knowledge on the PNLSS equations. Although an analytical PNLSS expression Eq. (18) exists, we have to note that the exact coefficients are not retrieved (error does not drop to zero). This is the price to be payed for allowing for a larger flexibility in the nonlinear terms (degrees 0 till 5, while degree 3 would suffice), combined with the strong dependence on the initialisation of the linear part.

Model estimation of c y from CFD data
In this section, a model describing the kinematics of an oscillating cylinder in a fluid flow is estimated from CFD data.

Initial values of the model: the Best Linear Approximation
As in the Van der Pol example in Section 4.1, the linear model is estimated using 7 random phase multisine realisations. From each realisation one period is regarded as transients and is not used for the identification. In this case, contrary to the Van der Pol example where the input and output had no physical meaning, each realisation corresponds to a trajectory of an imposed cylinder oscillation. The frequency band covered ranges from 0 to 1.5 f St , ( f St = 3 Hz), at a frequency resolution of 0.1 Hz and an amplitude level of A/D = 0.3. In Fig. 8, one period of the input and the output, as well as their amplitude spectra, averaged over the successive periods, are plotted. The DFT of the input and the output are denoted Y and C Y , respectively.
Given the nonparametric estimate (Ĝ BLA ), found through application of Eq. (12) to (14), a parametric FIR (Finite Impulse Response) model was estimated. A FIR model was preferred to ensure stability [44]. The selection of the model order is based on a trade off between the cost of the error of the fit (cost) and the number of state variables required. The higher the model order, the higher the number of state variables resulting in a combinatorial growth of monomials in ζ(t) and η(t). A plot of the cost as a function of the model order is given in Fig. 9. An FIR model of 5th order was selected. Fig. 10 shows the nonparametric and the parametric estimateĜ BLA andĜ BLA,par . The red line indicates the variation found over the different phase realisations that were applied. It is a measure for the combined effect of measurement noise and stochastic nonlinear influences (due to the random selection of the phases). Discrepancies between the nonparametric and the parametric estimates can be tolerated given the high uncertainty of the linear estimate due to nonlinear effects. The green line shows the noise floor. All components which are present in the signal but do not act at an integer multiple of f 0 will cause leakage. This will make up most of the power in what is considered to be noise. From Eq. 2 it follows that the Strouhal frequency is equal to 3 Hz. In practice we can however observe that this will not be exactly equal to 3. This will cause leakage to occur. Besides leakage also numerical averaging can cause fluctuations in the solution of the CFD simulations. It can be seen that especially in a band around f St the total level of distortion is of equivalent power as the function itself, emphasising the need for a nonlinear model. The parametric FIR estimate, written in state-space formulation, will serve as initial values for the A, B, C and D matrices in Eq. (6).

The estimation dataset
The envisaged application of the identified model is the simulation of the forces on real-life systems such as an oil riser or a chimney stack. Generally, apart from gusts, the operation regime would consist of a constant (or slowly varying) flow speed resulting in a constant frequency of excitation and thus oscillation. Therefore, accurate simulation of the response to single sine oscillations will be most important.
When composing the training set, a trade-off must be made between signals that lean towards the intended operating condition of the model, which is a user requirement, and the compactness and richness of the data, a system identification requirement. To that end, sine-sweep excitations with a slow sweep rate were found most adequate. Including multisine excitations has also been tested but did not yield improved validation results on single sine validation tests, and for that reason they were omitted. The final training dataset was obtained by gradually increasing the amplitude density of the sine sweeps in order to improve the validation results. The sweeps cover the same frequency band as the multisines that were used during the linear identification step: 0 to 1.5 f S t (0 -4.5 Hz with a period T 0 = 16 s for a single sweep direction) at amplitude levels from A/D = 0.05 to 0.3 in discrete steps of 0.025. Table 2 gives an overview of the type of excitation signal and the identification technique used throughout the identification process. Fig. 11 shows the imposed oscillation (on top) of one part of the training data (A/D = 0.25). In the middle, the corresponding c y is plotted and the bottom figure shows the relative instantaneous frequency (relative to f St ) regarding both quantities. Looking at the relative instantaneous frequency we can clearly observe synchronisation of   the frequency of c y between, roughly, 10 s and 26 s. Within this region, rapid variation in magnitude of c y occurs. The different behaviour of the force coefficient during sweep up and down also highlights the hysteretic nature of the phenomenon. Oscillations in the relative instantaneous frequency are due to the presence of multiple frequency components in the signal. All sweep experiments end with a period in which the cylinder is held stationary so that the vortex shedding can relax to its autonomous behaviour. The stable limit cycle period is needed to ensure the smooth concatenation of different data parts in the training set (Fig. 12).

Training of the model
Training of the model on the estimation data is done following the scheme in Fig. A.1 of Appendix A. Table 3 gives an overview of how the error drops during successive iterations of the nonlinear optimisation over subsets of the coefficients. The selected subsets of the state and the output equation are indicated under 'Optimisation choices'. Note that the iteration number in the left column refers to a cycle of the loop in Fig. A.1 and not to the underlying Levenberg-Marquardt algorithm which is executed in each iteration and typically ran for 1000 steps. The relative rms error, e rms , between the modelled output and the true output of the estimation set is given by: where N i is the number of points of the i th data part and N is the total number of points. If, instead of iteratively optimising over subsets of the coefficients, all coefficients would have been tuned simultaneously, the optimisation would have hit a local minimum corresponding to a final error of e rms = 49.9 %, pointing out the relevance of the   Table 3: List of the errors during training of the model on the estimation data.
iterative approach. In Fig. 12, the input (y-displacement of the cylinder), output (force coefficient) and the error on the modelled output are plotted.

4.2.4.
Estimating the initial state of the model As described in Section 3.5, the initial state and input, x 0 and u 0 also need to be estimated from the data. It is known that nonlinear systems can have a non-unique solution following a bifurcation point. In [45] it was shown that the initial conditions, i.e. the starting time of the imposed motion with respect to the zero crossing of the autonomous oscillating force coefficient, can have a tremendous influence on the transient time before a steady state solution is obtained. The transients can range from only 15 motion cycles up to 200 cycles during which the solution changes continuously and unpredictably between cycles with positive and negative energy transfer. Moreover a complex higher order bifurcation structure was revealed while sweeping though the lock-in range at a constant amplitude level. In [46,47] a stability analysis of the flow about imposed oscillating circular cylinders was conducted. It was found that different stable states existed, the distinction being based on the bifurcation that arises. One of the bifurcations is of the Neimark-Sacker type. This leads to quasi-periodic states. But also outside lock-in it can be clear that the initial condition influences the quasi-periodic force coefficient. Since within this zone the wake will not synchronise to the cylinder motion, the starting time of the imposed motion will have a direct impact on the obtained solution. This is visualised by plotting the output to an imposed sweep and hold experiment for different initial conditions, Fig. 13. The impact of errors on x 0 and u 0 will thus also be significant.
The test signal that is used is composed of two parts. In the first part the frequency sweeps up to f ex / f St = 0.5, from then on the frequency is held constant. A high amplitude level of A / D = 0.25 is applied. In the top part of Fig. 13, the output of the model for two different x 0 values is plotted. The bottom part shows the instantaneous frequency of the input. u 0 is set to zero since at all time steps prior to the start of the input, the cylinder was indeed stationary. Clearly, the influence of the initial state remains present throughout the entire signal length. We want to stress that although it cannot be guaranteed that the exact bifurcation structure of the underlying system is also present in the model, this influence of the initial state is clearly a system property and not due to numerics.

Model validation
To study the validity of the estimated model, it is subjected to various validation experiments. How to validate the model strongly depends on the intended application, and is therefore a user choice. In this work, two user requirements are put forward: 1. We want the model to perform well on signals which relate closely to how VIV is observed in practice. Therefore the validation experiments will be of the type of single sine oscillations. 2. As a second requirement, we want a single model to be valid over a wide range in frequency and amplitude of the input signal. The set of validation experiments will thus consist of a large number single sine oscillations of different frequencies and amplitudes.
As quality measures, four quantities are computed for each validation experiment: • The correlation coefficient R • The relative error on the maximum of the amplitude, This is an interesting measure since the maximum amplitude is a strong design criterium.
• We also want to introduce a quality measure which is phase independent, since we have seen that the model is very sensitive to errors on the estimate of the initial state (Section 4.2.4) and since we know that such errors will mainly manifest themselves as phase errors. This is obtained by transforming the error to the frequency domain. A frequency-domain error, based only on the magnitudes of the DFT of the true output, C Y , and the simulated output, C Y mod is given by: where ω s = 2π f s , f s being the sampling frequency.
An overview of the applied validation experiments and the corresponding errors is given in Table 4.  Table 4: Validation results of mono sine oscillations.
Each validation experiment starts with a sweep up, until the desired frequency is reached. For the remaining part of time the frequency is held constant. The quality measures are calculated only on the part where the desired frequency is reached and held constant.
By plotting the true output and the simulated output together as time series it was noticed that in practically all cases the model was able to accurately reproduce the maximum amplitude, the correct frequency and even the characteristic features of the signal. In Fig. 14a, the simulated output of validation experiment 21 is plotted. Although the case is outside of lock-in, a strong interaction between the natural vortex shedding and the imposed oscillation is visible. The quasi-periodicity becomes even more clear when looking at the PSD and the phase plot of c y (Fig. 14b). The model is able to reproduce this quasi-periodic oscillation with an accuracy of 95 %.
In a few cases, higher values of e rms were observed. This increased error appears to be due to small phase errors between the true and the simulated output, for which e rms is a very sensitive measure (an example of such case suffering from a phase error is given in Fig. 15a). Phase errors are easily introduced due to the sensitivity on the initial state (see Section 4.2.4). As a result we can have two solutions of the model, which track the same limit cycle, although they will never be equal to each other at a single time instant.
From a user point of view, having small phase errors on the output will in many cases be tolerable, given that a good estimate of the maxium amplitude, the frequency of the response and the feature of the signal are obtained. In these cases, the quality of the simulated output can best be quantified in terms of e DFT .
In Fig. 15b, the time domain output of the case which performs the worst in terms of e DFT is plotted (experiment 17). The oscillating behaviour of the instantaneous frequency of c y measured from CFD suggests that this case resides on the border of lock-in, exhibiting quasi-periodic oscillations (see Section 4.2.4).This region of parameter space is also known as the transition region (transition to lock-in) [48]. Although this very sensitive instability is also captured by the model, the entering and exiting of synchronisation is out of phase with respect to the true output. Altering the initial conditions of the CFD set up, i.e. the point at which the cylinder oscillation starts interacting with the natural vortex shedding, revealed that the fluid is also able to stabilise to multiple slightly different solutions. This highlights the sensitivity of this particular part of parameter space.
Since the correlation coefficient R has a mean value of 0.96 with only a limited spread we can conclude that the modelled output is in good relation with the true output. The estimated model is, in other words, able to reproduce the distinct behaviour of the force coefficient both inside and outside the lock-in region for cylinders oscillating at a single frequency. A contour plot of the relative rms error (Fig. 16a) and of R (Fig. 16b) visualises the results in the frequency-amplitude plane.
Considering potential applications of the model, e.g. for design purposes, the relative error on the maximum of the amplitude might well be the most crucial quality measure. A low mean value of 3.2 % and only a small spread can be reported here (Fig. 17).

Conclusions
In this work, a discrete PNLSS model that relates the displacement of an oscillating cylinder in a fluid flow to the fluid forces acting on the cylinder is identified. The model has been trained on a concatenation of imposed swept sine oscillations and it has been validated on single sine experiments, as single sines are considered to be the intended application of the model. Data, both for training and for validation were obtained using CFD simulations at the fixed free-stream condition corresponding to a Reynolds number of 100. The PNLSS model is initialised using the best linear approximation, obtained from exciting the system with random phase multisine realisations. Subsequently, an iterative optimisation strategy is used to avoid poor local minima. The identification approach has numerically been validated on data obtained from the Van der Pol equation. The quality of the fluid-force model was assessed based on: the relative rms error, the relative error on the maximum amplitude, the correlation coefficient and a phase independent error. It was shown that the output of the model was in good agreement with the CFD data over a wide range of oscillation frequencies and amplitudes. Sensitive regions in the frequency-amplitude plane, especially on the edge of the lock-in region, were observed to be hard to simulate in terms of phase. The feature and the maximum amplitude, however, remained accurate.
• All monomials with only contributions in the input.
• All monomials that are affine in the states. Monomials where only one state appears, and it appears linearly.
• All monomials without cross products.
• All monomials that are full state-affine. Monomials where at most one state appears, and if it appears, it appears linearly.
• All monomials with total degree up to 3.
• All monomials with only contributions in the states.
• Only odd nonlinear degrees.
• All monomials except for the DC component (all states as well as the input raised to the power 0).
Clearly, also other subsets could have been proposed. Given the fact that 10 possible subsets were proposed (of which the diagonal option is only valid for the state equation, the E matrix) and not necessarily the same choice has to be made for the state as well as for the output equation, we arrive at a total of 90 possibilities (paths on the cost function) to be tested during each optimisation run.