Parameter Identiﬁcation in Nonlinear Mechanical Systems with Noisy Partial State Measurement Using PID-Controller Penalty Functions

: Dynamic models of physical systems often contain parameters that must be estimated from experimental data. In this work, we consider the identiﬁcation of parameters in nonlinear mechanical systems given noisy measurements of only some states. The resulting nonlinear optimization problem can be solved efﬁciently with a gradient-based optimizer, but convergence to a local optimum rather than the global optimum is common. We augment the dynamic equations with a morphing parameter and a proportional–integral–derivative (PID) controller to transform the objective function into a convex function; the global optimum can then be found using a gradient-based optimizer. The morphing parameter is used to gradually remove the PID controller in a sequence of steps, ultimately returning the model to its original form. An optimization problem is solved at each step, using the solution from the previous step as the initial guess. This strategy enables use of a gradient-based optimizer while avoiding convergence to a local optimum. The efﬁcacy of the proposed approach is demonstrated by identifying parameters in the van der Pol–Dufﬁng oscillator, a hydraulic engine mount system, and a magnetorheological damper system. Our method outperforms genetic algorithm and particle swarm optimization strategies, and demonstrates robustness to measurement noise. Pol–Dufﬁng oscillator using proportional–derivative controller (PD), genetic algorithm (GA), and particle swarm optimization (PSO) approaches.


Introduction
To generate a simulation of a mechanical system, one first requires a mathematical model of the system's dynamics. The parameters for dynamic system models are often estimated from experimental measurements, a process called parameter estimation or parameter identification. This task is difficult in general, particularly when the mathematical model is nonlinear. Additional challenges are presented when only some of the system's states are measured and when experimental data are corrupted by measurement noise. Parameter identification with partial state measurement is highly relevant in many areas of the biological and physical sciences [1].
Parameter identification is typically approached as a least-squares optimization problem, wherein the model parameters must be determined so as to minimize the difference between the measured response of the physical system and the simulated response of the analogous mathematical model [2,3]. Classical optimization techniques include the Gauss-Newton and Levenberg-Marquardt algorithms.
These gradient-based methods are very effective provided the objective function is convex; if it is not, they often converge to a local optimum rather than the global optimum [4]. In many mathematical models, the parameters appear nonlinearly with respect to the system states and the objective function is non-convex (i.e., it contains many local optima). Stochastic optimization techniques such as simulated annealing, the genetic algorithm (GA), and particle swarm optimization (PSO) are capable of finding the globally optimal solution to a non-convex problem; however, if a large number of parameters must be identified, these methods tend to become computationally expensive [5,6]. In general, parameter identification problems can be solved more efficiently if a gradient-based algorithm can be used even for non-convex optimization problems.
Homotopy methods have led to effective global optimization algorithms [7]. Vyasarayani et al. explored the application of homotopy optimization to identifying parameters in nonlinear systems, including those governed by differential algebraic equations [8]. In the homotopy optimization method, a high-gain observer and a morphing parameter are introduced into the dynamic equations (and, thus, into the objective function implicitly). This transformation makes the objective function convex [9] and enables use of a gradient-based optimization method. The dynamic equations are gradually morphed back to the original system by reducing the value of the morphing parameter in a sequence of steps. A new optimization problem is solved at each step, using the solution from the previous step as the initial guess. Provided the morphing parameter is reduced in sufficiently small steps, convergence to a local optimum is avoided despite using a simple gradient-based optimizer. This process is continued until the morphing parameter reaches zero, thereby recovering the original optimization problem, but now it is solved starting from a good initial guess.
In the present work, we use a proportional-integral-derivative (PID) controller to implement a homotopy optimization strategy and study the effectiveness of this approach at identifying parameters in the van der Pol-Duffing oscillator, a hydraulic engine mount system, and a magnetorheological (MR) damper system. Accurate mathematical models for engine mount and MR damper systems are important tools for engineering analysis, design, and control applications; thus, fast and accurate parameter identification methods are critical. The PID controller used in this work is a natural extension of the Luenberger observer [10]. We consider partial state measurement, wherein experimental measurements for only some of the system states are available, as well as the effect of measurement noise. We also compare the performance of the proposed approach to the performance of popular GA [11][12][13] and PSO [14][15][16] methods.
The remainder of the paper is organized as follows. The proposed parameter identification strategy is presented in Section 2. Numerical examples are presented and discussed in Section 3. Finally, conclusions are provided in Section 4.

Methodology
We assume the dynamics of the experimental system are governed by a system of ordinary differential equations, written here in first-order form: . , x ne (t)] T contains the time series of the n position-level states and p [p 1 , p 2 , . . . , p k ] is a vector of the k parameters in the experimental system. In the following, we assume that only one state variable (x 1e (t)) is measured; the other states in the experimental system are unobserved. The mathematical model representing the experimental system (Equation (1a), (1b)) is written in an analogous form:ẋ where x 1 (t) [x 1 (p, t), x 2 (p, t), . . . , x n (p, t)] T . Our notation here emphasizes that the trajectory of the mathematical model is a function of the unknown parametersp [p 1 ,p 2 , . . . ,p k ]. The objective is to identify the parameters in the mathematical model (p in (Equation (2a), (2b)) to minimize the error between the experimental measurement (x 1e (t)) and the corresponding predicted response (x 1 (p, t)).

Optimization Strategy
We identify model parametersp by minimizing the following objective function: In general, this objective function surface will contain many local minima; a naive optimization strategy may lead to any one of them. To avoid convergence to a local minimum, we introduce a morphing parameter (λ) and a proportional-integral-derivative (PID) controller into the mathematical model (Equation (2a), (2b)): where k p , k i and k d are the proportional, integral and derivative gains, respectively, and is the error we wish to minimize. The controller used here (with k i = 0) resembles the Luenberger observer [10] used in the state estimation of linear systems. The predicted response and objective function are now implicit functions of morphing parameter λ: The optimization process begins by providing an arbitrary initial guess for the parameters and setting morphing parameter λ = 1. Provided the proportional gain (k p ) is sufficiently high, the PID controller in Equation (4a), (4b) will make the objective function surface convex [9], thereby enabling use of any gradient-based optimizer. The optimizer uses the objective function J(p, 1) to seek values for the model parameters (p) such that the predicted response (x 1 (p, 1, t)) matches the measured response (x 1e (t)). We denote the solution from this step 1p * . We then reduce λ by a small amount δ and minimize J(p, 1 − δ), using 1p * as the initial guess. The solution from this step ( 1−δp * ) is then used as the initial guess for the subsequent step, where λ = 1 − 2δ. This process is repeated until λ = 0, whereupon the PID controller vanishes from the mathematical model (Equation (4)). In the final step, λ = 0 and the global minimum to the original problem ( 0p * ) is obtained.
The proposed approach closely resembles the penalty function method. Any constrained optimization problem can be transformed into an unconstrained optimization problem by treating each constraint violation as a penalty. Each penalty i contributes a new term to the objective function, scaled by a weighting parameter r i . Values are selected for each r i and the optimization problem is solved. If the violation of a constraint from the original problem is too large, the corresponding weighting parameter is increased and the optimization problem is solved again, using the previous solution as the initial guess. This process is continued until all constraint violations are acceptably low. Penalty functions and weighting parameters r i are analogous to the PID controller and morphing parameter λ in Equation (4). However, note that the proposed approach is able to find the global optimum, whereas the penalty function method arrives at only a near-optimal solution in general.

Error Convergence
To gain insight into the error dynamics of the proposed approach, we subtract Equation (4a), (4b) from Equation (1a), (1b) to obtain the following error equation: Equation (6) represents a damped second-order system with stiffness k p , damping k d , and an equilibrium point at e 1 = 0. The k p , k i , and k d terms affect the error dynamics as seen in a typical PID controller: the proportional gain determines sensitivity to position-level error, the integral gain determines sensitivity to steady-state error, and the derivative gain determines sensitivity to the rate at which the error is changing. (We note that, for the systems presented in this work, the integral term was found to provide only marginal benefit.) For sufficiently high values of k p , the proportional term dominates and will drive the error to zero regardless of discrepancies between the forcing functions on the right-hand side of Equation (6). A proportional gain should be selected at the outset (i.e., when λ = 1) so that the synchronization error [17] is small when the initial guess for the parameter vector is used in the mathematical model (Equation (4a), (4b)).

Results
In this section, we identify parameters in the van der Pol-Duffing oscillator, a hydraulic engine mount system, and a magnetorheological damper system. We explore the effect of measurement noise and provide comparisons with two popular stochastic optimizers: the genetic algorithm and particle swarm optimization strategies.

Van der Pol-Duffing Oscillator
The van der Pol-Duffing (VDPD) oscillator represents a family of classical nonlinear systems with rich dynamical properties, and they appear in many fields of science and engineering [18][19][20][21].
The experimental system we consider in this section is an externally-excited VDPD oscillator [22], given here in first-order form: where f (t) = g cos(ωt) is the external excitation, with amplitude g = 0.50 and frequency ω = 0.79. The behavior of the VDPD oscillator depends on three important variables: α, β, and µ. Parameter α determines the stiffness of the system, β determines the amount of nonlinearity, and µ > 0 is the damping coefficient. VDPD oscillators are typically classified into three categories: single-well (α and β both positive), double-well (α negative and β positive), and double-hump (α positive and β negative).
In this example, we seek to identify the system parameters α, β, and µ given the trajectory of x 1e (t).
We assume the initial conditions and the external excitation ( f (t)) are known. As stated in Section 2, we identify parameters for the mathematical model by minimizing Equation (5). In this example, the parameters arep α,β,μ . We assume bounds are known for each parameter; thus, p l ≤p ≤ p u where p l and p u are vectors of the lower and upper bounds, respectively. We solve a constrained optimization problem because the VDPD oscillator loses stability beyond particular ranges of parameter values. We consider one experimental system of each category (single-well, double-well, and double-hump); the system parameters and bounds are provided in Table 1.
We first consider the single-well case, where the true parameters are [α, β, µ] = [0.5, 0.5, 0.1]. The initial conditions for the experimental system and mathematical model are assumed to be [1,0]. We generate synthetic experimental data by numerically integrating Equation (7a), (7b) for 20 s [22] using a fourth-order Runge-Kutta method and storing only x 1e (t), as shown in Figure 1a. The shape of the objective function is shown in Figure 1b as a function of β, for four values of α and with µ = 0.1. There are clearly many peaks and valleys in J(β) and, thus, there are many local minima in the objective function J(p). We show the basins of attraction in Figure 1c and confirm that a simple gradient-based optimizer is very likely to converge to a local minimum. To avoid converging to a local minimum, we use the strategy described in Section 2 with a high-gain PD controller. We use initial parameter guesses at the upper bounds, as listed in Table 1 (i.e., α,β,μ = [10, 10,2]) and set δ = 0.2. The objective function is given by Equation (5) and is solved using the "fmincon" constrained optimization solver in MATLAB R ; we use the sequential quadratic programming (SQP) algorithm with a first-order optimality tolerance of 10 −6 as the termination criterion. The mathematical model is as follows: Proportional and derivative gains are set to k p = k d = 10, and these values are confirmed to achieve synchronization between the experimental and predicted responses when λ = 1. As shown in Figure 2a, all parameters are successfully identified. The performance of GA and PSO strategies is shown in Figure 2b,c for comparison. We used the "ga" and "pso" functions in MATLAB R with a population size of 30 in each case (i.e., 10n p where n p is the number of parameters to be identified [23]). Note that the GA and PSO strategies are also successful at identifying the parameters in the single-well case. We now repeat the process for double-well and double-hump VDPD oscillators; the experimental data, objective function shape, and basins of attraction for these systems are shown in Figure 1d-i. The double-hump oscillator has a highly unstable response, so we use a simulation duration of 2.5 s in this case. The performance of the PD controller, GA, and PSO strategies is shown in Figure 2d-f for the double-well case and in Figure 2g-i for the double-hump case. All three methods are again successful at identifying the parameters in the double-well case. However, notice that the GA and PSO both fail to find the global optimum in the double-hump case, converging to α,β,μ = [1.3818, −0.3617, 0.1581] (maximum relative error of 58%) and [1.4754, −0.4767, 0.0975] (maximum relative error of 5%), respectively. Notice the relatively small basin of attraction in Figure 1i, which suggests that the parameter identification problem is relatively challenging in the double-hump case. Even when the GA and PSO strategies are successful, they are substantially more computationally expensive than the proposed PD controller approach. As shown in Figures 3 and 4, the GA and PSO strategies require more iterations and more function evaluations (i.e., they must consider more candidate solutions) than the proposed approach. Table 2 summarizes the performance of the three strategies.   Finally, we consider the robustness of the proposed approach to uncertainty in the measured data, since experimental measurements are typically corrupted by some amount of noise. We apply Gaussian noise to each data point in the experimental response as follows: wherex 1e is the noise-corrupted measurement, NSR is the noise-to-signal ratio, and r is a zero-mean, normally-distributed random number. As shown in Table 3, the identified parameters are not substantially affected by zero-mean noise, as would be expected. As shown in Figure 5, the predicted response reliably tracks the experimental measurement with NSR = 5%, suggesting that the proposed approach is robust to measurement noise. Table 3. Performance of proportional-derivative controller optimization strategy for the van der Pol-Duffing oscillator with measurement noise (experimental parameter values shown in parentheses).

Hydraulic Engine Mount System
Engine mount systems are designed to isolate passengers in a vehicle from the mechanical vibration of the engine. Designers of these systems must consider two sources of excitation: the low-amplitude, high-frequency vibrations caused by the engine as well as the high-amplitude, low-frequency vibrations induced by the uneven road surface. The response of the engine mount system is critical in the development of engines, and an accurate mathematical model is essential. The mechanical details and operating principles of passive engine mount systems are described in the literature [24][25][26].
The dynamics of this system are governed by a set of differential algebraic equations [27]: where M(q, p) diag(m L , 0, m M , m H ) is the mass matrix, q [x 1e , x 2e , x 3e , x 4e ] T is the vector of generalized coordinates, p [d H2 , c E1 , c E2 , d E ] is the vector of parameters to be identified, C T q (q) is the constraint Jacobian, and L is the Lagrange multiplier. (For clarity of notation, we have not shown functional dependence of each state variable on time.) The generalized force vector is given as follows: In this example, the engine mount is subjected to a frequency-sweep excitation force of F(t) = 100 sin(4π · 25 t t). The elastomer force (F E ), membrane force (F M ), and hydro-damping force (F H ) depend on the stiffness and damping of linear and nonlinear components: The algebraic constraint equation that relates the generalized coordinates is as follows: which has the following Jacobian: Our objective is to identify the parameters p = [d H2 , c E1 , c E2 , d E ] from the experimental response. Once again, we assume that only x 1e (t) is measured and the other states are unobserved. In this example, we use 5-second simulations. The system parameters that are assumed to be known are provided in Table 4; the experimental values of the parameters to be identified are provided in Table 5 along with their bounds and initial guesses.  Table 5. Parameters to be identified in the hydraulic engine mount system model [27].

Parameter Units Experimental Value Initial Guess Bounds
To identify parameters p, we first express the dynamic (Equation (10a), (10b)) in first-order form and augment the system with a PD controller: is the Lagrange multiplier and we set gain k 1 = 10. We use the same methods as described in Section 3.1 and solve the optimization problem using the "fmincon" constrained optimization solver in MATLAB R , with a first-order optimality tolerance of 10 −10 as the termination criterion. We first attempt to identify the parameters without the PD controller by setting λ = 0; as shown in Figure 6a, the initial guess (i.e., all parameters at their lower bounds) leads to a local minimum. However, the parameters are correctly identified when the PD controller is used, as illustrated in Figure 6b. In this example, we observe a slight improvement in performance when a PID controller is used (i.e., when the integral term is included as shown in Equation (4a), (4b); we use integral gain k i = 0.02. As shown in Figure 6c and summarized in Table 6, the PID controller successfully identifies the unknown parameters and requires less computational effort than the PD controller. We also observe reasonable robustness of the proposed method to Gaussian measurement noise, as demonstrated by the results shown in Table 7. Note that parametersd H2 andd E are relatively small compared to the others and thus we observe greater sensitivity of these parameters to measurement noise. In practice, a low-pass filter or nondimensionalization may be used to reduce this relative sensitivity. Table 6. Performance of optimization strategies for the hydraulic engine mount system using proportional-derivative controller (PD), proportional-integral-derivative controller (PID), genetic algorithm (GA), and particle swarm optimization (PSO) approaches.  Table 7. Performance of proportional-derivative controller optimization strategy for the hydraulic engine mount system with measurement noise (experimental parameter values shown in parentheses).  Finally, we compare the performance of the proposed approach to that of the GA and PSO strategies. As before, we used the "ga" and "pso" functions in MATLAB R with a population size of 10n p in each case (where n p = 4 in this example). As shown in Figure 6d,e, the stochastic optimizers fail to correctly identify the parameters. Furthermore, as shown in Figure 7, they require more function evaluations than the PD and PID controller strategies. A summary of these comparisons can be found in Table 6. Thus, the PD and PID controllers demonstrate superior performance in terms of accuracy as well as computational expense.

Magnetorheological Damper System
Magnetorheological (MR) dampers are semi-active devices that use MR fluids to control vibration in structural and mechanical systems [28][29][30]. Adjusting the current delivered to an electromagnet changes its magnetic field, thereby changing the viscosity of the MR fluid. This simple control scheme makes MR dampers ideal for vibration control in vehicle systems [31,32]. MR dampers are also very reliable in practice as they simply act as passive dampers in the event of control system failure. However, for design and implementation in practical applications, one requires a precise mathematical model of the force-displacement behavior of the system. A suitable mathematical model must capture the nonlinear behavior of the system caused by hysteresis.
The Bouc-Wen model is widely used to model the hysteresis behavior of structural and mechanical components [33,34], and is often used to model MR damper systems. In this example, we consider an MR damper system that comprises three parallel force-generating components: where F e (t) is the total force applied by the MR damper system, x e (t) is the displacement of one end relative to the other (e.g., the relative vertical motion of the sprung and unsprung masses in a vehicle), c 0 and k 0 are damping and stiffness coefficients, and α represents the ratio between post-elastic stiffness to elastic stiffness of the damper. We do not consider the influence of the applied electrical current on the damping characteristics. The dynamics of the hysteresis component (z e (t)) are governed by the following differential equation: The parameters α, β, γ, A, and n determine the hysteresis behavior of the system. In this example, we identify the system parameters p [c 0 , k 0 , α, β, γ, A, n] assuming the total force (F e (t)) is measured. The system parameters and bounds are provided in Table 8. The system is excited by a sinusoidal input displacement of x e (t) = 1.5 sin(5πt), as shown in Figure 8. We use 1-second simulations in this example. To identify the system parameters, we rewrite Equations (18) and (19) in first-order form and add a PD controller: where x 1 (t) = 1.5 sin(5πt) and x 2 (t) is its time derivative. In this case, we set gain k 1 = 400 to achieve synchronization when λ = 1 and use the following objective function: The initial guess for each parameter (equal to its lower bound) is listed in Table 8.  The evolution of the optimization process using the PD controller is shown in Figure 9; the final identified parameters are listed in the last column of Table 8. Several of the identified parameters differ somewhat from the corresponding experimental values. However, when the identified parameters are used in the original mathematical model (Equations (18) and (19)), the response is in close agreement with the "experimental" response ( Figure 10). Thus, the identified parameter values nevertheless produce an accurate overall mathematical model. We again compare the performance of the PD controller with that of GA and PSO strategies [36][37][38], and we again observe that the PD controller obtains a lower objective function value in fewer function evaluations ( Figure 11). The comparison between the PD controller and stochastic optimization strategies is summarized in Table 9. Table 9. Performance of optimization strategies for the magnetorheological damper system using proportional-derivative controller (PD), genetic algorithm (GA), and particle swarm optimization (PSO) approaches.

Strategy Number of Iterations Number of Function Evaluations
Error (J)  Figure 9. Evolution of parameter estimates (normalized relative to experimental values) for magnetorheological damper system using proportional-derivative controller.

Figure 10.
Experimental data and estimated response using identified parameters for magnetorheological damper system: (a) damping force over time, (b) displacement-level hysteretic response, and (c) velocity-level hysteretic response. Figure 11. Performance of three strategies to identify parameters for magnetorheological damper system: (a) error convergence of objective function J(c 0 , k 0 , α, β, γ, A, n); (b) computational effort required (cumulative number of function evaluations).

Conclusions
In this paper, we have presented a homotopy optimization strategy that uses a proportional-integral-derivative (PID) controller as a penalty function. We have applied this strategy to identify parameters in the van der Pol-Duffing oscillator and two systems of practical importance in engineering applications: a hydraulic engine mount system and a magnetorheological damper system. Through detailed numerical examples, we have demonstrated the utility of the proposed approach to identify model parameters given noisy measurements of only a subset of the system states. We have also compared its performance to that of two popular stochastic optimizers, namely the genetic algorithm and particle swarm optimization strategies. In general, we found the proposed approach to be superior in terms of the solutions that were found as well as the computational effort required. For the hydraulic engine mount system presented in this work, the integral term was found to provide a marginal benefit over a proportional-derivative (PD) controller. We, therefore, recommend considering the PID controller optimization strategy for parameter identification tasks: it is robust, easy to implement, and efficient thanks to its use of a gradient-based optimizer. Future work includes investigating application of the proposed approach to more practical scenarios, such as with data obtained experimentally and approximate mathematical models. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.