# The Identification of Nonlinear Systems UsingC

Authors: Ming-Fei Chen; Dyi-Cheng Chen; Chung-Heng Yang; Shih-Feng Tseng
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

In order to instantaneously distinguish theC
t (coefficient of viscous damping) and K
t (coefficient of stiffness), which are both functions of time in an M.C.K. nonlinear system, a new identification method is proposed in this paper. The graphs of the C
t-K
t are analyzed and the dynamic behavior of M.C.K. systems in a C
t-K
t coordinate plane is discussed. This method calculates two adjacent sampling data, the displacement, velocity, and acceleration (which are obtained from the responses of a pulse response experiment) and then distinguishes C
t and K
t of an instantaneous system. Finally, this method is used to identify the aerostatic bearing dynamic parameters, C and K.


---

## Body

## 1. Introduction

In order to control an M.C.K. system, the stiffness and damping must be known. However, in an actual M.C.K. system, the stiffness and damping are essentially nonlinear. Even though several studies consider the frequency domain [1–6], these analyses do not allow real-time identification. References [7–10] identify the time domain. Huang et al. [11] modified the logarithmic-decrement method to avoid disturbance from measured noise, but the second peak is still required to analyze the damping. As a result, a system that has at least two peaks is an underdamped system. In order to control the system, the displacement response must be almost critically damped and the peaks of its displacement response must be less than one. Therefore, the use of the logarithmic-decrement method is limited. For real-time control systems, the current method requires improvement.A single degree-of-freedom linear oscillator with viscous damping has the equation of motion:(1)
M
x
¨
+
C
x
˙
+
K
x
=
f
(
t
)
.However, in reality, many such oscillators show some nonlinearity in stiffness and/or damping. This paper presents a simple method to identify this nonlinearity. If the mass is assumed to be known and constant, the motion can be measured, in response to a known external force,f. The displacement, the velocity, and the acceleration can be measured and calculated. At any particular time t, (1) has only two unknowns, K and C. The equation describes a straight line in the C
-
K plane. If the system is exactly described by (1), then the set of straight lines at successive times, t, must all pass through a single point in the plane, which then determines the values of C and K.If the system is nonlinear or departs from the assumptions of (1) in some other way, such as having nonviscous damping, the successive lines do not all pass through a single point. However, if the variation is reasonably smooth, they still define an envelope curve. This curve in the C
-
K plane gives the time-dependent values of the equivalent linearized C and K. This paper presents an algorithm to map this envelope and illustrates the method with numerical and experimental examples.


## 2. Theoretical Background and Methodology

### 2.1. TheC
(
t
)
-
K
(
t
) Coordinate Plane

IfM
=
1 kg, C
=
2 (N·sec/m), and K
=
4 (N/m), then (
2,4
) are the rectangular Cartesian coordinates of a point that has C as the abscissa and K as the ordinate, as shown in Figure 1.Figure 1
TheC
-
K coordinate plane.SinceC
(
t
), K
(
t
) for M.C.K. system values are both functions of time, C
(
t
), K
(
t
) form a trajectory on the C
(
t
)
-
K
(
t
) plane of coordinate. This study terms the characteristic trajectory of an M.C.K. system, as shown in Figure 2. Sections 5.2 and 5.3 show a simulation and experimental examination of this situation. Only C
(
t
) and K
(
t
) are constant in the time domain. The characteristic trajectory is shown as Figure 1.Figure 2
C
(
t
) and K
(
t
) are both functions of time on the C
(
t
)
-
K
(
t
) coordinate plane.


### 2.2. The Limits of the Critical Damping

For a critically damped system, the damping ratio is(2)
ξ
=
1
=
C
2
K
M
and the stiffness is derived as
(3)
K
(
t
)
=
C
(
t
)
2
4
M
.If the mass in Figure2 is 1 kg, then the critical damping limit is shown by the dotted line in Figure 3. If K
(
t
)
=
0.5
+
0.2
t (N/m) and C
(
t
)
=
1
+
t
0.5 (N·sec/m), then (3) gives t
<
0.22774 (sec) for underdamping, t
=
0.227744 (sec) for critical damping, and t
>
0.227744 (sec) for overdamping.Figure 3
The critical damping limits for Figure2, with M
=
1 kg.


## 2.1. TheC
(
t
)
-
K
(
t
) Coordinate Plane

IfM
=
1 kg, C
=
2 (N·sec/m), and K
=
4 (N/m), then (
2,4
) are the rectangular Cartesian coordinates of a point that has C as the abscissa and K as the ordinate, as shown in Figure 1.Figure 1
TheC
-
K coordinate plane.SinceC
(
t
), K
(
t
) for M.C.K. system values are both functions of time, C
(
t
), K
(
t
) form a trajectory on the C
(
t
)
-
K
(
t
) plane of coordinate. This study terms the characteristic trajectory of an M.C.K. system, as shown in Figure 2. Sections 5.2 and 5.3 show a simulation and experimental examination of this situation. Only C
(
t
) and K
(
t
) are constant in the time domain. The characteristic trajectory is shown as Figure 1.Figure 2
C
(
t
) and K
(
t
) are both functions of time on the C
(
t
)
-
K
(
t
) coordinate plane.


## 2.2. The Limits of the Critical Damping

For a critically damped system, the damping ratio is(2)
ξ
=
1
=
C
2
K
M
and the stiffness is derived as
(3)
K
(
t
)
=
C
(
t
)
2
4
M
.If the mass in Figure2 is 1 kg, then the critical damping limit is shown by the dotted line in Figure 3. If K
(
t
)
=
0.5
+
0.2
t (N/m) and C
(
t
)
=
1
+
t
0.5 (N·sec/m), then (3) gives t
<
0.22774 (sec) for underdamping, t
=
0.227744 (sec) for critical damping, and t
>
0.227744 (sec) for overdamping.Figure 3
The critical damping limits for Figure2, with M
=
1 kg.


## 3. Methodology

The general movement formula for an M.C.K. system is(4)
M
x
¨
+
C
x
˙
+
K
x
=
f
(
t
)
,
where M is a mass, C is a damping coefficient, K is a stiffness coefficient, and f
(
t
) is an external force. Equation (1) represents a two-dimensional equation (L
(
t
)) at time, t, on the C
(
t
)
-
K
(
t
) coordinate plane:
(5)
M
x
¨
+
C
(
t
)
x
˙
+
K
(
t
)
x
=
f
(
t
)
⋯
L
(
t
)
.In (5), C
(
t
), K
(
t
) are taken to be functions of time and are generally considered constants in the definition of the time domain. It is known that M, C
(
t
), K
(
t
), and f
(
t
) for M.C.K. systems, the displacement response, x
(
t
), the velocity, x
˙
(
t
), and the acceleration, x
¨
(
t
), at any time, t, are necessary for (5); namely,
(6)
M
x
¨
(
t
)
+
C
(
t
)
x
˙
(
t
)
+
K
(
t
)
x
(
t
)
=
f
(
t
)
⋯
L
(
t
)
.It is known that at time,t
1, the displacement response of the system is x
(
t
1
). After differentiation, x
˙
(
t
1
) and x
¨
(
t
1
) are known as are M and f
(
t
1
), but C
(
t
1
)and K
(
t
1
) are unknown. Equation (6) can be rewritten as the slope-intercept from an equation of the line:
(7)
K
(
t
1
)
=
-
x
˙
(
t
1
)
x
(
t
1
)
C
(
t
1
)
+
(
-
M
x
¨
(
t
1
)
+
f
(
t
1
)
x
(
t
1
)
)
⋯
L
(
t
1
)
.
At another time, t
2, L
(
t
2
) represents the linear equation of the M.C.K. motion formula on the C
(
t
)
-
K
(
t
) plane coordinates:
(8)
K
(
t
2
)
=
-
x
˙
(
t
2
)
x
(
t
2
)
C
(
t
2
)
+
(
-
M
x
¨
(
t
2
)
+
f
(
t
2
)
x
(
t
2
)
)
⋯
L
(
t
2
)
.Figure4 shows that from times, t
1 to t
2, C
(
t
) and K
(
t
) for an M.C.K. system have a linear binary equation, L
(
t
1
) and L
(
t
2
), at the locations of the C
(
t
)
-
K
(
t
) plane coordinates. In terms of the dynamics, the K
(
t
) intercept of the line is driven by the displacement response, by the acceleration, and by external forces from b
(
t
1
) to b
(
t
2
), and the slope is driven by the displacement response and the velocity from m
(
t
1
) to m
(
t
2
). C
(
t
), K
(
t
) for a linear binary equation for the motion equation of an M.C.K. system change the location of L
(
t
1
) to L
(
t
2
), which results in (
C
(
t
1
)
,
K
(
t
1
)
) and (
C
(
t
2
)
,
K
(
t
2
)
), respectively. If C
(
t
), K
(
t
) for the system are constant in the time domain, that is, (
C
(
t
1
)
,
K
(
t
1
)
)
=
(
C
(
t
2
)
,
K
(
t
2
)
)
=
(
C
,
K
), then (
C
,
K
) mark the intersection point with L
(
t
1
) and L
(
t
2
), as shown in Figure 4. C
(
t
), K
(
t
) for a linear binary equation for the motion equation of an M.C.K. system rotate on the C
(
t
)
-
K
(
t
) plane about the rotational center (
C
,
K
).Figure 4
The dynamic behavior of the slope-intercept from the equation of the line on theC
(
t
)
-
K
(
t
) plane coordinates, from times, t
1 to t
2.IfC
(
t
) and K
(
t
) for the system are both functions of time, t, in the continuous time domain, it is assumed that
(9)
t
1
=
t
,
t
2
=
t
+
Δ
t
.AsΔ
t
→
0, it is inevitable that (
C
(
t
)
,
K
(
t
)
)
=
(
C
(
t
→
t
+
Δ
t
)
,
K
(
t
→
t
+
Δ
t
)
)
=
(
C
(
t
+
Δ
t
)
,
K
(
t
+
Δ
t
)
) and satisfies L
(
t
) and L
(
t
+
Δ
t
), such that (
C
(
t
→
t
+
Δ
t
)
,
K
(
t
→
t
+
Δ
t
)
) is an intersection point with L
(
t
) and L
(
t
+
Δ
t
), as shown in Figure 5. (
C
(
t
→
t
+
Δ
t
)
,
K
(
t
→
t
+
Δ
t
)
) is the instantaneous center of rotation on the C
(
t
)
-
K
(
t
) plane coordinates for the motion formula for an M.C.K. system. The center also forms a trajectory that defines the coordinates of the characteristic trajectory. Equations (7) and (8) are rewritten as follows:(10)
K
(
t
→
t
+
Δ
t
)
=
-
x
˙
(
t
)
x
(
t
)
C
(
t
→
t
+
Δ
t
)
+
(
-
M
x
¨
(
t
)
+
f
(
t
)
x
(
t
)
)
⋯
L
(
t
)
,
K
(
t
→
t
+
Δ
t
)
=
-
x
˙
(
t
+
Δ
t
)
x
(
t
+
Δ
t
)
C
(
t
→
t
+
Δ
t
)
+
(
-
M
x
¨
(
t
+
Δ
t
)
+
f
(
t
+
Δ
t
)
x
(
t
+
Δ
t
)
)
⋯
L
(
t
+
Δ
t
)
.Figure 5
The dynamic behavior of the slope-intercept from the equation of the line on theC
(
t
)
-
K
(
t
) plane coordinates. C
(
t
) and K
(
t
) for the system are both functions of time, t, in the continuous time domain.(
C
(
t
→
t
+
Δ
t
)
,
K
(
t
→
t
+
Δ
t
)
) is solved as the simultaneous solution of (10), so the time-domain displacement response (the velocity and the acceleration calculated by differentiating the displacement) obtained by experiment can be used to solve C
(
t
) and K
(
t
) for an M.C.K. system. The entire method is described as follows:(11)
M
x
¨
+
x
˙
C
(
t
)
+
x
K
(
t
)
=
f
(
t
)
t
=
[
t
0
,
t
1
,
t
2
,
…
,
t
n
]
,
t
>
0
,
x
¨
=
[
a
t
0
a
t
1
a
t
2
⋮
a
t
n
]
=
[
a
(
t
)
]
,
x
˙
=
[
v
t
0
v
t
1
v
t
2
⋮
v
t
n
]
=
[
v
(
t
)
]
,
x
=
[
x
t
0
x
t
1
x
t
2
⋮
x
t
n
]
=
[
x
(
t
)
]
,
K
(
t
)
=
[
k
(
t
0
→
t
1
)
,
k
(
t
1
→
t
2
)
,
k
(
t
2
→
t
3
)
,
…
,
k
(
t
n
-
1
→
t
n
)
]
=
[
k
(
t
)
]
,
C
(
t
)
=
[
c
(
t
0
→
t
1
)
,
c
(
t
1
→
t
2
)
,
c
(
t
2
→
t
3
)
,
…
,
c
(
t
n
-
1
→
t
n
)
]
=
[
c
(
t
)
]
,
[
v
(
t
)
]
·
[
c
(
t
)
]
+
[
x
(
t
)
]
·
[
k
(
t
)
]
=
-
M
·
[
a
(
t
)
]
+
[
f
(
t
)
]
,
{
v
t
0
c
(
t
0
→
t
1
)
+
x
t
0
k
(
t
0
→
t
1
)
=
-
M
·
a
t
0
+
f
(
t
0
)
v
t
1
c
(
t
0
→
t
1
)
+
x
t
1
k
(
t
0
→
t
1
)
=
-
M
·
a
t
1
+
f
(
t
1
)
}
⟹
c
(
t
0
→
t
1
)
,
k
(
t
0
→
t
1
)
{
v
t
1
c
(
t
1
→
t
2
)
+
x
t
1
k
(
t
1
→
t
2
)
=
-
M
·
a
t
1
+
f
(
t
1
)
v
t
2
c
(
t
1
→
t
2
)
+
x
t
2
k
(
t
1
→
t
2
)
=
-
M
·
a
t
2
+
f
(
t
2
)
}
⟹
c
(
t
1
→
t
2
)
,
k
(
t
1
→
t
2
)
{
v
t
2
c
(
t
2
→
t
3
)
+
x
t
2
k
(
t
2
→
t
3
)
=
-
M
·
a
t
2
+
f
(
t
2
)
v
t
3
c
(
t
2
→
t
3
)
+
x
t
3
k
(
t
2
→
t
3
)
=
-
M
·
a
t
3
+
f
(
t
3
)
}
⟹
c
(
t
2
→
t
3
)
,
k
(
t
2
→
t
3
)
⋮
{
v
t
n
-
1
c
(
t
n
-
1
→
t
n
)
+
x
t
n
-
1
k
(
t
n
-
1
→
t
n
)
=
-
M
·
a
t
n
-
1
+
f
(
t
n
-
1
)
v
t
n
c
(
t
n
-
1
→
t
n
)
+
x
t
n
k
(
t
n
-
1
→
t
n
)
=
-
M
·
a
t
n
+
f
(
t
n
)
}
⟹
c
(
t
n
-
1
→
t
n
)
,
k
(
t
n
-
1
→
t
n
)
.As is seen, the values ofC
(
t
) and K
(
t
) from the last formula cannot be determined. If necessary, the number of samples can be easily adjusted to satisfy the demands of the analysis.In theory, this method is very sensitive to noise. Since the velocity and acceleration are obtained by differentiating the displacement, the displacement, the velocity, and the acceleration must be differentiable and continuous. In order to ensure the feasibility of this method, the measurement data must be smoothed.


## 4. Inference

For example, for the impact response of the experiment to momentumF
^, the initial conditions are x
(
0
)
=
0 and x
˙
(
0
+
)
=
F
^
/
M.Ifξ
=
0
(
C
=
0
) represents no damping, then x
(
t
), x
˙
(
t
), and x
¨
(
t
) are given as follows:
(12)
x
(
t
)
=
F
^
M
ω
n
sin
(
ω
n
t
)
,
x
˙
(
t
)
=
F
^
M
cos
⁡
(
ω
n
t
)
,
t
>
0
,
x
¨
(
t
)
=
-
F
^
ω
n
M
sin
(
ω
n
t
)
.Consequently, Figure5 gives
(13)
m
(
t
)
=
-
x
˙
(
t
)
x
(
t
)
=
ω
n
cot
(
ω
n
t
)
=
K
M
cot
(
ω
n
t
)
,
(14)
b
(
t
)
=
-
M
x
¨
(
t
)
x
(
t
)
=
M
ω
n
2
=
K
.Equation (14) defines the relationship between the stiffness coefficient, K, and the intercept, b
(
t
), on the K
(
t
) axis.If0
<
ξ
<
1, so the system is underdamped; then x
(
t
), x
˙
(
t
), and x
¨
(
t
) are given as follows:
(15)
x
(
t
)
=
F
^
M
ω
d
e
-
ξ
ω
n
t
sin
(
ω
d
t
)
,
x
˙
(
t
)
=
F
^
M
e
-
ξ
ω
n
t
cos
⁡
(
ω
d
t
)
-
ξ
ω
n
·
x
(
t
)
,
t
>
0
,
x
¨
(
t
)
=
-
ω
n
2
(
1
-
2
ξ
2
)
·
x
(
t
)
-
2
ξ
ω
n
F
^
M
e
-
ξ
ω
n
t
cos
⁡
(
ω
d
t
)
.Consequently, Figure5 gives
(16)
m
(
t
)
=
-
x
˙
(
t
)
x
(
t
)
=
ω
d
cot
(
ω
n
t
)
-
ξ
ω
n
,
(17)
b
(
t
)
=
-
M
x
¨
(
t
)
x
(
t
)
=
M
ω
n
2
(
1
-
2
ξ
2
)
+
2
M
ξ
ω
d
ω
n
cot
(
ω
d
t
)
=
K
(
1
-
2
ξ
2
)
+
C
ω
d
cot
(
ω
d
t
)
.Ifξ
=
0, (16) and (17) are, respectively, equal to (13) and (14).Ifξ
>
1, which means that the system is overdamped, then x
(
t
), x
˙
(
t
), and x
¨
(
t
) are given as follows:
(18)
x
(
t
)
=
F
^
M
2
ω
n
ξ
2
-
1
e
-
ξ
ω
n
t
(
e
ω
n
ξ
2
-
1
t
-
e
ω
n
ξ
2
-
1
t
)
,
x
˙
(
t
)
=
F
^
M
2
ω
n
ξ
2
-
1
e
-
ξ
ω
n
t
×
[
(
-
ξ
ω
n
+
ω
n
ξ
2
-
1
)
e
ω
n
ξ
2
-
1
t
-
(
-
ξ
ω
n
-
ω
n
ξ
2
-
1
)
e
-
ω
n
ξ
2
-
1
t
]
,
x
¨
(
t
)
=
F
^
M
2
ω
n
ξ
2
-
1
e
-
ξ
ω
n
t
×
[
(
-
ξ
ω
n
+
ω
n
ξ
2
-
1
)
2
e
ω
n
ξ
2
-
1
t
-
(
-
ξ
ω
n
-
ω
n
ξ
2
-
1
)
2
e
-
ω
n
ξ
2
-
1
t
]
.Consequently, Figure8 gives
(19)
m
(
t
)
=
-
(
∫
(
(
-
ξ
ω
n
+
ω
n
ξ
2
-
1
)
e
ω
n
ξ
2
-
1
t
-
(
-
ξ
ω
n
-
ω
n
ξ
2
-
1
)
e
-
ω
n
ξ
2
-
1
t
)
×
(
e
ω
n
ξ
2
-
1
t
-
e
-
ω
n
ξ
2
-
1
t
)
-
1
∫
)
,
b
(
t
)
=
-
M
(
∫
(
(
-
ξ
ω
n
+
ω
n
ξ
2
-
1
)
2
e
ω
n
ξ
2
-
1
t
-
(
-
ξ
ω
n
-
ω
n
ξ
2
-
1
)
2
e
-
ω
n
ξ
2
-
1
t
)
×
(
e
ω
n
ξ
2
-
1
t
-
e
-
ω
n
ξ
2
-
1
t
)
-
1
∫
)
.Although this is a linear system, the nonlinear time-variant systems observed within a small time interval are almost all linear, time-invariant systems. Using the reasoning form, b, the relationship between the damping ratio, ξ, and the behavior of the system on the C
-
K plane is known. This determines the system damping ratio on the C
-
K plane. In the special case of the absence of damping, (
C
(
t
)
=
0
), b directly indicates the K value on the C
-
K plane. This means that an M.C.K. system can be directly observed on the C
-
K plane (M is a fixed value).


## 5. Application

The simulation response of an M.C.K. system is the measurement data obtained from this experiment. This study inversely calculates the two adjacent data points of the simulation response to determine theC
(
t
) and K
(
t
) values for the system. In Section 5.3 an experiment involving an aerostatic bearing system is performed to validate this method.


### 5.1. Theoretical Analysis of Pulse Response Simulation

Assuming that the parameters of the MCK areM
=
1 (kg), C
=
2 (N·sec/m), K
=
4 (N/m), F
^
=
-
0.01 (kg·m/sec), and f
(
t
>
0
)
=
0 (N), the initial conditions are x
(
0
)
=
0 (m), x
˙
(
0
+
)
=
-
0.01 (m/sec), time = 7 s, the displacement response of 71 points is taken, and the time space is 0.1 s. Equation (15) gives solutions for the displacement, the velocity, and the acceleration, as shown in (20), (21), and (22), respectively. Consider
(20)
x
(
t
)
=
0.01
3
e
-
t
sin
(
3
t
)
,
(21)
x
˙
(
t
)
=
0.01
e
-
t
cos
⁡
(
3
t
)
-
x
(
t
)
,
(22)
x
¨
(
t
)
=
-
2
x
(
t
)
-
0.02
e
-
t
cos
⁡
(
3
t
)
.Figures6(a), 6(b), and 6(c) show the displacement, the velocity, and the acceleration, respectively. Equation (11) gives the dynamic behavior of the slope-intercept from an equation of the line on the C
(
t
)
-
K
(
t
) plane coordinates, as shown in Figure 6(d). Each line has an instantaneous center of rotation at the C
(
t
)
-
K
(
t
) values and L1–L8 are used to label this portion of the data only. Figure 6(a) shows that the C
(
t
) and K
(
t
) values for this system occur before the first peak and that the K
(
t
) value is solved in Figure 6(e). That is, the K
(
t
) value is constant in the time domain. The C
(
t
) value is solved in Figure 6(f). This means that the C
(
t
) value in the time domain is also a constant. Table 1 shows the portion of the data, C
(
t
), K
(
t
), for time functions and the slope of the M.C.K. system in Figure 6, in which L
(
0.2
)
=
L
1, L
(
0.4
)
=
L
2, L
(
0.6
)
=
L
3
,
…, and so forth.Table 1
The portion of the data,C
(
t
), K
(
t
), from Figure 6(d).
Lt
t (sec)
x (m)
v (m/sec)
a (m/sec2)
b
(
t
)
m
(
t
)
K
(
t
)
(N/m)
C
(
t
) (N·sec/m)
L1
0.2
−0.0016
−0.0061
0.018612
11.59676
−3.79838
4
2
L2
0.4
−0.00247
−0.00269
0.015259
6.173171
−1.08659
4
2
L3
0.6
−0.00273
-
5
E
-
05
0.011026
4.036969
−0.01848
4
2
L4
0.8
−0.00255
0.001723
0.006754
2.64883
0.675585
4
2
L5
1
−0.0021
0.002687
0.003011
1.436505
1.281747
4
2
L6
1.2
−0.00152
0.002984
0.000111
0.072932
1.963534
4
2
L7
1.4
−0.00094
0.002795
−0.00185
−1.97586
2.987929
4
2
L8
1.6
−0.00042
0.002304
−0.00292
−6.92301
5.461503
4
2Using the time-domain displacement response from the Pulse Response Simulation to determine theC
(
t
), K
(
t
) values. (a) The time-domain displacement response is x. (b) The velocity response is v, which is the differential of x. (c) The acceleration response is a, which is the differential of v. (d) The characteristic trajectory of the motion formula on the C
(
t
)
-
K
(
t
) plane. (e) The coefficient of stiffness is solved as graph K
(
t
). (f) The coefficient of viscous damping is solvedas graph C
(
t
).
(a)
(b)
(c)
(d)
(e)
(f)Figure7(a) shows that C
(
t
) and K
(
t
) are solved by inverse calculation. That is, C
(
t
) and K
(
t
) for this system are constant, similar to the parameter settings for the response simulation. The C
(
t
), K
(
t
) coordinate is (
2,4
), which is on the upper left of the critical damping limits, so this is an underdamped system. The calculated margin of error is less than 10
-
14, so it can be ignored. Figures 7(b), 7(c), and 7(d) compare C
(
t
) and K
(
t
) with C
=
2 and K
=
4. In the plot, the last item is neglected because it cannot be determined at the end of the C
(
t
), K
(
t
) value.The theoretical analysis of the Pulse Response Simulation at theC
(
t
)
-
K
(
t
) coordinate plane. (a) The characteristic trajectory of the motion formula and the critical damping limit for the M.C.K. system on the C
(
t
)
-
K
(
t
) plane. (b) The error in the C
(
t
)
-
K
(
t
) values distributed on the C
(
t
)
-
K
(
t
) coordinate plane. (c) The time-domain error in K
(
t
). (d) The time-domain error in C
(
t
).
(a)
(b)
(c)
(d)Using the time-domain displacement response for the simulation experiment to solve theC
(
t
), K
(
t
) values, which are both functions of time. (a) The time-domain displacement response is x. (b) The velocity response is v, which is the numerical differential of x. (c) The acceleration response is a, which is the numerical differential of v. (d) The characteristic trajectory of the motion formula on the C
(
t
)
-
K
(
t
) plane. (e) The coefficient of stiffness is solved as graph, K
(
t
), K
(
t
), and is not constant. (f) The coefficient of viscous damping is solvedas graph, C
(
t
), C
(
t
), and is not constant.
(a)
(b)
(c)
(d)
(e)
(f)


### 5.2. Numerical Analysis of a Simulated Pulse Response

In order to confirm the feasibility of this method in identifying nonlinear systems, it is assumed that a curve imitates an impact test process, as shown in Figure8(a). This virtual curve must represent the trajectories of a nonlinear system, in order to verify the accuracy of the proposed method. For these reasons, a displacement theoretical equation for the virtual nonlinear system is presented in (23). The displacement trajectories for the virtual nonlinear systems are obtained from (23) and are shown in Figure 8(a). Using the finite difference method to obtain the speed in Figure 8(b) from the displacement curve in Figure 8(a), and the acceleration in Figure 8(c) from the speed in Figure 8(b). The methods already mentioned are then used to obtain the parameters, C
(
t
) and K
(
t
). Using (23), (24), and (25), another set of theoretical solutions for the displacement, the velocity, and the acceleration are then calculated. The theoretical finite difference solutions are compared with the theoretical solutions to verify the accuracy of the method. Because this is an envelope approximation, the smaller the time intervals, the more accurate the solutions. The accuracy of the calculation process depends on the number of iterations.Assuming thatM
=
1 (kg) and C
(
t
), K
(
t
) are unknown, the M.C.K. system receives a single pulse of force, with no other external forces. The simulation of the displacement response, defined by (23), is shown in Figure 8(a):
(23)
x
(
t
)
=
-
3
e
-
3
t
1.3
sin
(
2
t
)
,
(24)
x
˙
(
t
)
=
-
117
100
t
0.3
x
t
-
6
e
-
3
t
1.3
sin
(
2
t
)
,
(25)
x
¨
(
t
)
=
(
-
351
300
t
-
0.7
+
4563
300
t
0.3
-
4
)
x
t
+
234
5
t
0.3
e
-
3
t
1.3
cos
⁡
(
2
t
)
.The simulation takes 4 s and the time space is 0.001 s. Although (23) can be deduced as the velocity of (24) and as the acceleration of (25), the displacement response is used by the numerical differential to obtain the velocity and the acceleration during the actual experiment, as shown in Figures 8(b) and 8(c). The initial conditions are obtained from the displacement response and the velocity, so the size of the pulse can be ignored. The status of other forces, which is f
(
t
>
0
)
=
0 (N), is the only other consideration. Figure 8(d) shows that (11) solves the C
(
t
), K
(
t
) values as a characteristic trajectory on the C
(
t
)
-
K
(
t
) plane. Table 2 shows the portion of the C
(
t
), K
(
t
) values that represent time functions and the slope of the M.C.K. system in Figure 8, in which L
(
0.01
)
=
L
1, L
(
0.03
)
=
L
2, L
(
0.09
)
=
L
3
,
…
,
L
(
0.99
)
=
L
8. Figures 8(e) and 8(f) show that the K
(
t
) and C
(
t
) values can change in the time domain. In the plot, the last value for C
(
t
), K
(
t
) cannot be directly solved. However, by using the Central Difference Method, the velocity can be determined on the first application. The acceleration can be calculated from this velocity data, using the same method, so the first two and last three items are not considered. When the two adjacent sampling points have been calculated, the system can determine instantaneous C
(
t
), K
(
t
) values. Data 2 in Figure 9(a) shows that (23) derives the velocity for (24) and the acceleration for (25) by using 4001 points to determine the C
(
t
), K
(
t
) values of the system, and data 1 is obtained by the numerical differentiation of the velocity and the acceleration. Figures 9(b), 9(c), and 9(d) show that the error is reduced after 10
-
2 s of the initial response, which is seen by comparing data 1 with data 2. Data 3 in Figure 9(a) shows that M
=
1 (kg) at the critical damping limit for the M.C.K. system. The characteristic trajectory of C
(
t
), K
(
t
) increases with time and approaches the critical damping limits. That is, the damping ratio increases with time and approaches the critical damping ratio and the coordinate of the characteristic trajectory for C
(
t
), K
(
t
) in the upper left boundary, so the M.C.K. system is underdamped at the time of the simulation. This is similar to the design of a shock absorber system that includes variable stiffness and damping. Firstly, the response of the system is determined and then the differentiation method previously discussed is used to determine the necessary stiffness and damping.Table 2
The portion of the data,C
(
t
)
,
K
(
t
), for the time functions from Figure 8(d).
Lt
t (sec)
x (m)
v (m/sec)
a (m/sec2)
b
(
t
)
m
(
t
)
K
(
t
) (N/m)
C
(
t
) (N·sec/m)
L1
0.01
−0.05955
−5.89536
13.58601
228.1615
−99.0059
70.05914
1.596898
L2
0.03
−0.17433
−5.5664
18.55715
106.4511
−31.9311
35.84346
2.211251
L3
0.09
−0.47109
−4.28554
22.77945
48.35459
−9.09703
20.07025
3.109183
L4
0.19
−0.78703
−2.07592
20.35073
25.85776
−2.63768
15.25057
4.021417
L5
0.33
−0.90441
0.198556
11.87931
13.1349
0.219543
14.23512
5.011415
L6
0.51
−0.73229
1.434002
2.598733
3.548783
1.958248
15.49573
6.100836
L7
0.73
−0.40647
1.351987
−2.25806
−5.55523
3.326129
18.51146
7.235646
L8
0.99
−0.14247
0.67754
−2.37718
−16.6854
4.755633
22.61141
8.263207The simulated pulse response on theC
(
t
)
-
K
(
t
) plane. (a) The characteristic trajectory of the motion formula and the critical damping limit for the M.C.K. system on the C
(
t
)
-
K
(
t
) plane. (b) The error in the C
(
t
), K
(
t
) values on the C
(
t
)
-
K
(
t
) coordinate plane. (c) The time-domain error in K
(
t
). (d)The time-domain error in C
(
t
).
(a)
(b)
(c)
(d)


### 5.3. An Aerostatic Bearing System

During the transient response interval, an aerostatic bearing has nonlinear stiffness and damping characteristics. Therefore an experiment to determine the dynamic parameters for an aerostatic bearing is proposed. The aerostatic bearing is a closed loop feedback system. When it is located at a balance point, the gas film gap remains constant. When the load changes, the aerostatic bearing air gap has a new equilibrium position.In general, an aerostatic bearing’s dynamic performance depends on the dynamic parameters of the bearing, such as the stiffness and damping values. Since finding an exact solution for these values involves complex fluid mechanics problems, it is difficult to determine them. In order to verify the nonlinear dynamic parameter identification method presented in this paper, an experiment to measure stiffness and damping values of the aerostatic bearing was performed. Figures10 and 11 show a schematic diagram and a photograph of the experimental system. A free-falling steel ball is the pulse force, which acts on the aerostatic bearing and excites transient responses. Using the transient responses and the presented method, the bearing dynamic characteristics, such as the stiffness and damping, can be determined. The experimental results are shown in Figure 12.Figure 10
The impact response experiment for the aerostatic bearing system.Figure 11
A photograph of the impact response experiment.TheC, K values for the aerostatic bearing system on the C
-
K plane.
(a)
(b)
(c)
(d)
(e)
(f)For the experimental results of the displacement, the smoothing process is performed using the MATLAB software function, “polyfit,” as shown by the red line in Figure12(a). On the time axis, the smooth curve deviates from the actual measured curve after 0.025 seconds, because there is insufficient data. As it is close to the steady state, it does not affect the determination of the transient response.Figure12(a) section [A, E] shows a complete waveform, which marks the major part of the transient response of this system. In section [A, E] the fitted curve and the actual measured values are in good agreement and this is the interval of interest. Intervals [A, B] and [D, E] are the range of the gas film squeeze process and interval [B, D] is the range of the gas film gap expansion process. Points A, C, and E are the aerostatic bearing equilibrium positions where the displacement is zero.Figure12(b) shows that interval [A, E] is a continuous curve whose track is formed by an envelope. Figures 12(c) and 12(d) show the change in the values of K and C in the time domain. Intervals [A, B] and [B, D] are smooth curves, which show that the change in K
(
t
) is inversely proportional to the air gap shown in Figure 12(e). In these sections, the stiffness of the aerostatic bearing decreases. At point B, the air gap has the smallest value, but the stiffness, K
(
t
), has a maximum value. However, at point C, the air gap has the largest value, but the stiffness K
(
t
) is a minimum.The curve at interval [A, E], shown in Figure12(f), shows the relationship between C
(
t
) and the air gap, which is continuous and nonlinear. Interval [B, D] shows the rising state of the aerostatic bearing. During interval [B, C] it is seen that C
(
t
) is inversely proportional to the absolute value of the air gap. However, during the interval [C, D] it is seen that C
(
t
) is proportional to the absolute value of the air gap. Intervals [A, B] and [D, E] show the relationship between the air gas and C
(
t
) during compression of the bearing.


## 5.1. Theoretical Analysis of Pulse Response Simulation

Assuming that the parameters of the MCK areM
=
1 (kg), C
=
2 (N·sec/m), K
=
4 (N/m), F
^
=
-
0.01 (kg·m/sec), and f
(
t
>
0
)
=
0 (N), the initial conditions are x
(
0
)
=
0 (m), x
˙
(
0
+
)
=
-
0.01 (m/sec), time = 7 s, the displacement response of 71 points is taken, and the time space is 0.1 s. Equation (15) gives solutions for the displacement, the velocity, and the acceleration, as shown in (20), (21), and (22), respectively. Consider
(20)
x
(
t
)
=
0.01
3
e
-
t
sin
(
3
t
)
,
(21)
x
˙
(
t
)
=
0.01
e
-
t
cos
⁡
(
3
t
)
-
x
(
t
)
,
(22)
x
¨
(
t
)
=
-
2
x
(
t
)
-
0.02
e
-
t
cos
⁡
(
3
t
)
.Figures6(a), 6(b), and 6(c) show the displacement, the velocity, and the acceleration, respectively. Equation (11) gives the dynamic behavior of the slope-intercept from an equation of the line on the C
(
t
)
-
K
(
t
) plane coordinates, as shown in Figure 6(d). Each line has an instantaneous center of rotation at the C
(
t
)
-
K
(
t
) values and L1–L8 are used to label this portion of the data only. Figure 6(a) shows that the C
(
t
) and K
(
t
) values for this system occur before the first peak and that the K
(
t
) value is solved in Figure 6(e). That is, the K
(
t
) value is constant in the time domain. The C
(
t
) value is solved in Figure 6(f). This means that the C
(
t
) value in the time domain is also a constant. Table 1 shows the portion of the data, C
(
t
), K
(
t
), for time functions and the slope of the M.C.K. system in Figure 6, in which L
(
0.2
)
=
L
1, L
(
0.4
)
=
L
2, L
(
0.6
)
=
L
3
,
…, and so forth.Table 1
The portion of the data,C
(
t
), K
(
t
), from Figure 6(d).
Lt
t (sec)
x (m)
v (m/sec)
a (m/sec2)
b
(
t
)
m
(
t
)
K
(
t
)
(N/m)
C
(
t
) (N·sec/m)
L1
0.2
−0.0016
−0.0061
0.018612
11.59676
−3.79838
4
2
L2
0.4
−0.00247
−0.00269
0.015259
6.173171
−1.08659
4
2
L3
0.6
−0.00273
-
5
E
-
05
0.011026
4.036969
−0.01848
4
2
L4
0.8
−0.00255
0.001723
0.006754
2.64883
0.675585
4
2
L5
1
−0.0021
0.002687
0.003011
1.436505
1.281747
4
2
L6
1.2
−0.00152
0.002984
0.000111
0.072932
1.963534
4
2
L7
1.4
−0.00094
0.002795
−0.00185
−1.97586
2.987929
4
2
L8
1.6
−0.00042
0.002304
−0.00292
−6.92301
5.461503
4
2Using the time-domain displacement response from the Pulse Response Simulation to determine theC
(
t
), K
(
t
) values. (a) The time-domain displacement response is x. (b) The velocity response is v, which is the differential of x. (c) The acceleration response is a, which is the differential of v. (d) The characteristic trajectory of the motion formula on the C
(
t
)
-
K
(
t
) plane. (e) The coefficient of stiffness is solved as graph K
(
t
). (f) The coefficient of viscous damping is solvedas graph C
(
t
).
(a)
(b)
(c)
(d)
(e)
(f)Figure7(a) shows that C
(
t
) and K
(
t
) are solved by inverse calculation. That is, C
(
t
) and K
(
t
) for this system are constant, similar to the parameter settings for the response simulation. The C
(
t
), K
(
t
) coordinate is (
2,4
), which is on the upper left of the critical damping limits, so this is an underdamped system. The calculated margin of error is less than 10
-
14, so it can be ignored. Figures 7(b), 7(c), and 7(d) compare C
(
t
) and K
(
t
) with C
=
2 and K
=
4. In the plot, the last item is neglected because it cannot be determined at the end of the C
(
t
), K
(
t
) value.The theoretical analysis of the Pulse Response Simulation at theC
(
t
)
-
K
(
t
) coordinate plane. (a) The characteristic trajectory of the motion formula and the critical damping limit for the M.C.K. system on the C
(
t
)
-
K
(
t
) plane. (b) The error in the C
(
t
)
-
K
(
t
) values distributed on the C
(
t
)
-
K
(
t
) coordinate plane. (c) The time-domain error in K
(
t
). (d) The time-domain error in C
(
t
).
(a)
(b)
(c)
(d)Using the time-domain displacement response for the simulation experiment to solve theC
(
t
), K
(
t
) values, which are both functions of time. (a) The time-domain displacement response is x. (b) The velocity response is v, which is the numerical differential of x. (c) The acceleration response is a, which is the numerical differential of v. (d) The characteristic trajectory of the motion formula on the C
(
t
)
-
K
(
t
) plane. (e) The coefficient of stiffness is solved as graph, K
(
t
), K
(
t
), and is not constant. (f) The coefficient of viscous damping is solvedas graph, C
(
t
), C
(
t
), and is not constant.
(a)
(b)
(c)
(d)
(e)
(f)


## 5.2. Numerical Analysis of a Simulated Pulse Response

In order to confirm the feasibility of this method in identifying nonlinear systems, it is assumed that a curve imitates an impact test process, as shown in Figure8(a). This virtual curve must represent the trajectories of a nonlinear system, in order to verify the accuracy of the proposed method. For these reasons, a displacement theoretical equation for the virtual nonlinear system is presented in (23). The displacement trajectories for the virtual nonlinear systems are obtained from (23) and are shown in Figure 8(a). Using the finite difference method to obtain the speed in Figure 8(b) from the displacement curve in Figure 8(a), and the acceleration in Figure 8(c) from the speed in Figure 8(b). The methods already mentioned are then used to obtain the parameters, C
(
t
) and K
(
t
). Using (23), (24), and (25), another set of theoretical solutions for the displacement, the velocity, and the acceleration are then calculated. The theoretical finite difference solutions are compared with the theoretical solutions to verify the accuracy of the method. Because this is an envelope approximation, the smaller the time intervals, the more accurate the solutions. The accuracy of the calculation process depends on the number of iterations.Assuming thatM
=
1 (kg) and C
(
t
), K
(
t
) are unknown, the M.C.K. system receives a single pulse of force, with no other external forces. The simulation of the displacement response, defined by (23), is shown in Figure 8(a):
(23)
x
(
t
)
=
-
3
e
-
3
t
1.3
sin
(
2
t
)
,
(24)
x
˙
(
t
)
=
-
117
100
t
0.3
x
t
-
6
e
-
3
t
1.3
sin
(
2
t
)
,
(25)
x
¨
(
t
)
=
(
-
351
300
t
-
0.7
+
4563
300
t
0.3
-
4
)
x
t
+
234
5
t
0.3
e
-
3
t
1.3
cos
⁡
(
2
t
)
.The simulation takes 4 s and the time space is 0.001 s. Although (23) can be deduced as the velocity of (24) and as the acceleration of (25), the displacement response is used by the numerical differential to obtain the velocity and the acceleration during the actual experiment, as shown in Figures 8(b) and 8(c). The initial conditions are obtained from the displacement response and the velocity, so the size of the pulse can be ignored. The status of other forces, which is f
(
t
>
0
)
=
0 (N), is the only other consideration. Figure 8(d) shows that (11) solves the C
(
t
), K
(
t
) values as a characteristic trajectory on the C
(
t
)
-
K
(
t
) plane. Table 2 shows the portion of the C
(
t
), K
(
t
) values that represent time functions and the slope of the M.C.K. system in Figure 8, in which L
(
0.01
)
=
L
1, L
(
0.03
)
=
L
2, L
(
0.09
)
=
L
3
,
…
,
L
(
0.99
)
=
L
8. Figures 8(e) and 8(f) show that the K
(
t
) and C
(
t
) values can change in the time domain. In the plot, the last value for C
(
t
), K
(
t
) cannot be directly solved. However, by using the Central Difference Method, the velocity can be determined on the first application. The acceleration can be calculated from this velocity data, using the same method, so the first two and last three items are not considered. When the two adjacent sampling points have been calculated, the system can determine instantaneous C
(
t
), K
(
t
) values. Data 2 in Figure 9(a) shows that (23) derives the velocity for (24) and the acceleration for (25) by using 4001 points to determine the C
(
t
), K
(
t
) values of the system, and data 1 is obtained by the numerical differentiation of the velocity and the acceleration. Figures 9(b), 9(c), and 9(d) show that the error is reduced after 10
-
2 s of the initial response, which is seen by comparing data 1 with data 2. Data 3 in Figure 9(a) shows that M
=
1 (kg) at the critical damping limit for the M.C.K. system. The characteristic trajectory of C
(
t
), K
(
t
) increases with time and approaches the critical damping limits. That is, the damping ratio increases with time and approaches the critical damping ratio and the coordinate of the characteristic trajectory for C
(
t
), K
(
t
) in the upper left boundary, so the M.C.K. system is underdamped at the time of the simulation. This is similar to the design of a shock absorber system that includes variable stiffness and damping. Firstly, the response of the system is determined and then the differentiation method previously discussed is used to determine the necessary stiffness and damping.Table 2
The portion of the data,C
(
t
)
,
K
(
t
), for the time functions from Figure 8(d).
Lt
t (sec)
x (m)
v (m/sec)
a (m/sec2)
b
(
t
)
m
(
t
)
K
(
t
) (N/m)
C
(
t
) (N·sec/m)
L1
0.01
−0.05955
−5.89536
13.58601
228.1615
−99.0059
70.05914
1.596898
L2
0.03
−0.17433
−5.5664
18.55715
106.4511
−31.9311
35.84346
2.211251
L3
0.09
−0.47109
−4.28554
22.77945
48.35459
−9.09703
20.07025
3.109183
L4
0.19
−0.78703
−2.07592
20.35073
25.85776
−2.63768
15.25057
4.021417
L5
0.33
−0.90441
0.198556
11.87931
13.1349
0.219543
14.23512
5.011415
L6
0.51
−0.73229
1.434002
2.598733
3.548783
1.958248
15.49573
6.100836
L7
0.73
−0.40647
1.351987
−2.25806
−5.55523
3.326129
18.51146
7.235646
L8
0.99
−0.14247
0.67754
−2.37718
−16.6854
4.755633
22.61141
8.263207The simulated pulse response on theC
(
t
)
-
K
(
t
) plane. (a) The characteristic trajectory of the motion formula and the critical damping limit for the M.C.K. system on the C
(
t
)
-
K
(
t
) plane. (b) The error in the C
(
t
), K
(
t
) values on the C
(
t
)
-
K
(
t
) coordinate plane. (c) The time-domain error in K
(
t
). (d)The time-domain error in C
(
t
).
(a)
(b)
(c)
(d)


## 5.3. An Aerostatic Bearing System

During the transient response interval, an aerostatic bearing has nonlinear stiffness and damping characteristics. Therefore an experiment to determine the dynamic parameters for an aerostatic bearing is proposed. The aerostatic bearing is a closed loop feedback system. When it is located at a balance point, the gas film gap remains constant. When the load changes, the aerostatic bearing air gap has a new equilibrium position.In general, an aerostatic bearing’s dynamic performance depends on the dynamic parameters of the bearing, such as the stiffness and damping values. Since finding an exact solution for these values involves complex fluid mechanics problems, it is difficult to determine them. In order to verify the nonlinear dynamic parameter identification method presented in this paper, an experiment to measure stiffness and damping values of the aerostatic bearing was performed. Figures10 and 11 show a schematic diagram and a photograph of the experimental system. A free-falling steel ball is the pulse force, which acts on the aerostatic bearing and excites transient responses. Using the transient responses and the presented method, the bearing dynamic characteristics, such as the stiffness and damping, can be determined. The experimental results are shown in Figure 12.Figure 10
The impact response experiment for the aerostatic bearing system.Figure 11
A photograph of the impact response experiment.TheC, K values for the aerostatic bearing system on the C
-
K plane.
(a)
(b)
(c)
(d)
(e)
(f)For the experimental results of the displacement, the smoothing process is performed using the MATLAB software function, “polyfit,” as shown by the red line in Figure12(a). On the time axis, the smooth curve deviates from the actual measured curve after 0.025 seconds, because there is insufficient data. As it is close to the steady state, it does not affect the determination of the transient response.Figure12(a) section [A, E] shows a complete waveform, which marks the major part of the transient response of this system. In section [A, E] the fitted curve and the actual measured values are in good agreement and this is the interval of interest. Intervals [A, B] and [D, E] are the range of the gas film squeeze process and interval [B, D] is the range of the gas film gap expansion process. Points A, C, and E are the aerostatic bearing equilibrium positions where the displacement is zero.Figure12(b) shows that interval [A, E] is a continuous curve whose track is formed by an envelope. Figures 12(c) and 12(d) show the change in the values of K and C in the time domain. Intervals [A, B] and [B, D] are smooth curves, which show that the change in K
(
t
) is inversely proportional to the air gap shown in Figure 12(e). In these sections, the stiffness of the aerostatic bearing decreases. At point B, the air gap has the smallest value, but the stiffness, K
(
t
), has a maximum value. However, at point C, the air gap has the largest value, but the stiffness K
(
t
) is a minimum.The curve at interval [A, E], shown in Figure12(f), shows the relationship between C
(
t
) and the air gap, which is continuous and nonlinear. Interval [B, D] shows the rising state of the aerostatic bearing. During interval [B, C] it is seen that C
(
t
) is inversely proportional to the absolute value of the air gap. However, during the interval [C, D] it is seen that C
(
t
) is proportional to the absolute value of the air gap. Intervals [A, B] and [D, E] show the relationship between the air gas and C
(
t
) during compression of the bearing.


## 6. Conclusions

This paper presents a novel and effective method for the identification of nonlinear systems, usingC
(
t
)
-
K
(
t
) plane coordinates. A new time-domain analysis is proposed. In the C
(
t
)
-
K
(
t
) plane coordinates, the dynamic behavior of an M.C.K. system provides useful information. In C
(
t
)
-
K
(
t
) plane coordinates, any motion formula for an M.C.K. system can be expressed as a line in a geometric plane (a two-dimensional space) and each line has an instantaneous center of rotation on the C
(
t
)
-
K
(
t
) coordinate plane in this system. The two adjacent sampling data are calculated and then the values of C
(
t
) and K
(
t
) for the instantaneous system are determined. C
(
t
) and K
(
t
) are both functions of time, so in a traditional M.C.K. system, C and K are fixed constants, but there may be exceptions. In other words, the design that allows variable stiffness and damping can be identified in the time domain.



---

*Source: 497872-2014-02-26.xml*