Simulation and Actual Examples

1 Simulation

Consider the embedded system shown in previous theory section:

In this setup, the plant \(P(s)\) is defined as:

\[P(s) = \frac{1}{Ms^2+Cs}\]

And the controller \(C(s)\) is a PD controller:

\[C(s) = K_P + \frac{K_Ds}{\tau_D s+1}\]

We will simulate this system using a Chirp signal \(\mathrm{Sig}(t)\) as the excitation input, where the frequency changes exponentially over time:

\[\mathrm{Sig}(t)=A\sin\left(2\pi f_{\mathrm{min}}\frac{k^t-1}{\ln(k)}\right)\]

Here, \(A\) is the amplitude of the Chirp signal, \(f_{\mathrm{min}}\) is the starting frequency at \(t=0\), and the frequency at any given time \(t\) is \(f_{\mathrm{min}}k^t\).

Note that the differentiator in the PD controller includes a time constant \(\tau_D\). This is because a pure differentiator \(s\) is mathematically impossible to implement in practice, so it is realized in this filtered form. Since this chapter focuses on system identification, the details of this implementation will be covered in the Controller Design chapter.

The Python execution module is provided below. This simulation implements the PD controller to stabilize the motor model, then applies a Chirp signal that sweeps exponentially from \(0.1[\mathrm{Hz}]\) to \(2.0[\mathrm{kHz}]\) over \(50[\mathrm{s}]\). The frequency response analyzed from the input-output data is then compared with the theoretical frequency response obtained via \(j\omega\) analysis. The PD controller is designed using a technique called Pole Placement, which will also be detailed in the Controller Design chapter.

Run_sysid.py

# Copyright (c) 2021 Koichi Sakata

import numpy as np
from scipy import signal
from control import matlab
from pylib_sakata import ctrl
from pylib_sakata import plot
from pylib_sakata import fft

print('Start simulation!')

# Common parameters
Ts = 1/8000
dataNum = 10000
freqrange = [1, 1000]
freq = np.logspace(np.log10(freqrange[0]), np.log10(freqrange[1]), dataNum, base=10)
s = ctrl.tf([1, 0], [1])
z = ctrl.tf([1, 0], [1], Ts)
print('Common parameters were set.')

# Plant model
M = 2.0
C = 10.0
K = 0.0
Pmechs = ctrl.tf([1.0], [M, C, K])
Dz = z**-3
Pnz = ctrl.c2d(Pmechs, Ts, method='zoh') * Dz
Pnz_frd = ctrl.sys2frd(Pnz, freq)
print('Plant model was set.')

# Design PD controller
freq1 = 10.0
freq2 = 10.0
zeta2 = 1.0
Cz = ctrl.pd(freq1, freq2, zeta2, M, C, K, Ts)
Cz_frd = ctrl.sys2frd(Cz, freq)
print('PID controller was designed.')

print('System identification simulation is running...')
Snz = ctrl.feedback(Pnz, Cz, sys='S')
SPnz = ctrl.feedback(Pnz, Cz, sys='SP')
t = np.linspace(0.0, 50, int(50/Ts))
chirp = signal.chirp(t, f0=0.1, f1=2000, t1=50, method='logarithmic', phi=-90)
u, tout, xout = matlab.lsim(ctrl.tf2ss(Snz), chirp)
y, tout, xout = matlab.lsim(ctrl.tf2ss(SPnz), chirp)

fft_axis, chirp_fft = fft.fft(chirp, Ts)
fft_axis, u_fft = fft.fft(u, Ts)
fft_axis, y_fft = fft.fft(y, Ts)

Pmeas_frd, coh = fft.tfestimate(u, y, freq, Ts)

print('Plotting figures...')
# Time response
fig = plot.makefig()
ax1 = fig.add_subplot(311)
ax2 = fig.add_subplot(312)
ax3 = fig.add_subplot(313)
plot.plot_xy(ax1, t, chirp, '-', 'm', 0.5, 1.0, [0, 50], ylabel='Sig [N]', title='Time response')
plot.plot_xy(ax2, tout, u, '-', 'b', 0.5, 1.0, [0, 50], ylabel='In [N]')
plot.plot_xy(ax3, tout, y*1.0e3, '-', 'b', 0.5, 1.0, [0, 50], xlabel='Time [s]', ylabel='Out [mm]')

# Plant
fig = plot.makefig()
ax_mag = fig.add_subplot(311)
ax_phase = fig.add_subplot(312)
ax_coh = fig.add_subplot(313)
plot.plot_tffrd(ax_mag, ax_phase, Pmeas_frd, '-', 'm', 1.5, 1.0, ax_coh=ax_coh, coh=coh, title='Frequency response of plant')
plot.plot_tffrd(ax_mag, ax_phase, Pnz_frd, '--', 'b', 1.5, 1.0, freqrange, legend=['Measurement', 'Model'])

print('Finished.')

Start simulation!
Common parameters were set.
Plant model was set.
PID controller was designed.
System identification simulation is running...
The common pole-zeros of the zpk model have been deleted.
Plotting figures...
Finished.

The first output shows the time response of the excitation signal \(\mathrm{Sig}(t)\), input signal \(\mathrm{In}(t)\), and output signal \(\mathrm{Out}(t)\). The second output compares the measured frequency response (analyzed from the time-series data) with the model’s frequency response. The results match almost perfectly, confirming that the system identification process based on input-output time responses is valid.

The Coherence shown in the bottom plot is an index indicating the degree of correlation between the input and output signals. It ranges from \(0\) to \(1\), where a value closer to \(1\) indicates a higher correlation. If the coherence is low, the reliability of the frequency response in that specific range is also low. In this case, since the plot is based on noise-free simulation data, the coherence stays consistently near \(1\).

2 Python Functions for Analysis

In actual system identification, rather than the simulations shown in previous sections, an excitation signal \(\mathrm{Sig}(t)\) is applied to the control system of the actual plant. To perform the analysis, we record the actual system’s input signal \(\mathrm{In}(t)\), output signal \(\mathrm{Out}(t)\), and a flag signal \(\mathrm{Flag}(t)\) that indicates the excitation period.

The Python package pylib-sakata provides functions to derive frequency response data from CSV files recorded and saved using the TwinCAT3 Measurement YT Scope Project.

pylib_sakata.meas.measdata2frd(filePath, inputName, outputName, flagName, freq, inputGain=1.0, outputGain=1.0, windivnum=4, overlap=0.5)

This function is for system identification from input and output time response data of measurement file.

• Parameters:
  • filePath: file path of measurement file (.csv, .txt, and .mat are supported.)
  • inputName: input data name in the measurement file
  • outputName: output data name in the measurement file
  • flagName: flag data name in the measurement file
  • freq: 1-D array frequency data [Hz]
  • inputGain: inputdata gain (Optional), Default: 1.0, unit of input can be fixed by this parameter.
  • outputGain: outputdata gain (Optional), Default: 1.0, unit of output can be fixed by this parameter.
  • dt: sampling time of the time response data
  • windivnum: number of windows to divide the time response data
  • overlap: overlap ratio divided time response data (0 <= overlap < 1)
• Returns:
  • freqresp: instance of FreqResp class
  • coh: 1-D array coherence data

class pylib_sakata.fft.FreqResp(freq, resp, dt=0)

• Parameters:
  • freq: 1-D array frequency data [Hz]
  • resp: 1-D array frequency response data [complex data]
  • dt: sampling time (Optional), Default: 0, set the value >= 0. If dt = 0, the system is continuous time system.

Examples

import numpy as np
from pylib_sakata import meas

freq = np.logspace(np.log10(1.), np.log10(1000.), 10000, base=10)
freqresp, coh = meas.measdata2frd('data/001-inject.csv', 'ServoOut', 'PosErrUm', 'FlagNoise', freq, 1., 1.e-6)

3 Actual Examples

Here are some examples of system identification performed using the Python functions introduced in in previous theory section.

3.1 System Identification of a Servomotor

Here is a practical example of system identification for a BECKHOFF servomotor, model AM8111-1F21.

3.1.1 Measurement during Excitation

The following video shows the data acquisition process during Chirp signal excitation using the YT chart function of TwinCAT3 Scope View (TE13xx).

3.1.2 Analysis Software

The following Python module is used to perform system identification using the actual data obtained from the Chirp excitation. The dataset used is freq_resp_PD_20250903.csv, which can be found Here.

Run_sysid_exp.py

# Copyright (c) 2021 Koichi Sakata

import numpy as np
from scipy import signal
from control import matlab
from pylib_sakata import ctrl
from pylib_sakata import plot
from pylib_sakata import fft
from pylib_sakata import meas

print('Start simulation!')

# Common parameters
Ts = 1/8000
dataNum = 10000
freqrange = [1, 1000]
freq = np.logspace(np.log10(freqrange[0]), np.log10(freqrange[1]), dataNum, base=10)
s = ctrl.tf([1, 0], [1])
z = ctrl.tf([1, 0], [1], Ts)
print('Common parameters were set.')

# Plant model
M = 0.027
C = 0.7
K = 0.0
Pmechs = ctrl.tf([1], [M, C, K])
numDelay, denDelay = matlab.pade(Ts*4, n=4)
Ds = ctrl.tf(numDelay, denDelay)
Dz = z**-3
Pns = Pmechs * Ds
Pnz = ctrl.c2d(Pmechs, Ts, method='zoh') * Dz
Pnz_frd = ctrl.sys2frd(Pnz, freq)

print('Getting measurement data...')
measfileName = 'data/freq_resp_PD_20250903.csv'
# Frequency response
Pmeas_frd, coh = meas.measdata2frd(measfileName, 'ServoOutN[0]', 'ActPosUm[0]', 'FlagInject', freq, 1., 1.e-6, 8, 0.8)

# Time response
measdata = meas.getdata(measfileName)
time = measdata.time
NoiseOut = measdata.value[meas.getdataindex(measdata, 'NoiseOut[0]')]
ServoOutN = measdata.value[meas.getdataindex(measdata, 'ServoOutN[0]')]
ActPosUm = measdata.value[meas.getdataindex(measdata, 'ActPosUm[0]')]
FlagInject = measdata.value[meas.getdataindex(measdata, 'FlagInject')]

print('Plotting figures...')
# Time response
fig = plot.makefig()
ax1 = fig.add_subplot(411)
ax2 = fig.add_subplot(412)
ax3 = fig.add_subplot(413)
ax4 = fig.add_subplot(414)
plot.plot_xy(ax1, time, NoiseOut, '-', 'b', 0.5, 1.0, ylabel='[Nm]', legend=['Sig'], loc='upper right', title='Time response')
plot.plot_xy(ax2, time, ServoOutN, '-', 'b', 0.5, 1.0, ylabel='[Nm]', legend=['In'], loc='upper right')
plot.plot_xy(ax3, time, ActPosUm, '-', 'b', 0.5, 1.0, ylabel='[um]', legend=['Out'], loc='upper right')
plot.plot_xy(ax4, time, FlagInject, '-', 'b', 0.5, 1.0, xlabel='Time [s]', ylabel='[.]', legend=['Flag'], loc='upper right')

# Plant
fig = plot.makefig()
ax_mag = fig.add_subplot(311)
ax_phase = fig.add_subplot(312)
ax_coh = fig.add_subplot(313)
plot.plot_tffrd(ax_mag, ax_phase, Pmeas_frd, '-', 'm', 1.5, 1.0, ax_coh=ax_coh, coh=coh, title='Frequency response of plant')
plot.plot_tffrd(ax_mag, ax_phase, Pnz_frd, '--', 'b', 1.5, 1.0, freqrange, legend=['Measurement', 'Model'])

print('Finished.')

Start simulation!
Common parameters were set.
Getting measurement data...
Plotting figures...
Finished.

The first output represents the time response of the excitation signal \(\mathrm{Sig}(t)\), the input signal \(\mathrm{In}(t)\), the output signal \(\mathrm{Out}(t)\), and a flag signal (indicating excitation with a value of \(1\)). By using the Python functions introduced in the previous section, the actual frequency response characteristics can be derived from these measured time-series data. The second output shows a comparison between the frequency response of the model and the actual system.

By tuning the model’s mass (or inertia \(J\) for rotational systems) \(M\) and damping coefficient \(C\), the model’s frequency response is adjusted to match the actual plant, as shown in the plot. The closer they align, the more accurately the model and its parameters represent the real-world system’s characteristics. Once a satisfactory model is obtained, you can proceed to the controller design phase based on this validated model.

3.2 System Identification of a Servomotor with a Load

A servomotor is an actuator used to drive mechanical mechanisms. In practical field applications, motors are rarely operated alone; they are almost always connected to some form of mechanical system.

Here, we will perform system identification of a servomotor with a load, using a simple mechanism created with 3D CAD and attached to the motor shaft. You can find the 3D CAD models and G-code HERE.

This load mechanism is structured such that the shaft mounting member and the outer peripheral member are connected by four leaf springs. This design results in low stiffness only in the rotational direction, while other directions are constrained. Consequently, the system with this load mechanism attached behaves as a two-inertia system.

4 Two-Inertia System Model

A two-inertia system can be represented by the mechanical configuration shown below, consisting of a driving-side mass \(M_1\), a load-side mass \(M_2\), a spring \(K_\mathrm{reso}\) and damper \(C_\mathrm{reso}\) connecting them, an applied force \(u(=f)\), and positions \(y_1\) (driving side) and \(y_2\) (load side). This is referred to as a two-mass system in translational terms. For rotational systems, it is called a two-inertia system, and the same model applies by replacing mass with inertia and force with torque.

The model of a two-inertia system can be expressed as the linear sum of a rigid-body model \(P_{\mathrm{rigid}}(s)\) and a resonance model \(P_{\mathrm{reso}}(s)\):

\[\begin{aligned} P_1(s) &= \frac{y_1}{u} = P_{\mathrm{rigid}}(s) + k_1 P_{\mathrm{reso}}(s)\\ P_2(s) &= \frac{y_2}{u} = P_{\mathrm{rigid}}(s) + k_2 P_{\mathrm{reso}}(s) \end{aligned}\]

Here, \(P_1(s)\) is the transfer function from force (torque) \(u\) to the driving-side position \(y_1\), and \(P_2(s)\) is the transfer function from \(u\) to the load-side position \(y_2\). The rigid-body and resonance models are defined as follows:

\[\begin{aligned} P_{\mathrm{rigid}}(s) &= \frac{1}{(M_1+M_2)s^2+Cs}\\ P_{\mathrm{reso}}(s) &= \frac{1}{s^2 + 2\zeta\omega s+\omega^2} \end{aligned}\]

The mode influence constants \(k_1\) and \(k_2\) are determined by the relationship between the driving-side and load-side masses (inertias):

\[\begin{aligned} k_1 &= \frac{M_2}{M_1(M_1+M_2)}\\ k_2 &= \frac{-1}{M_1+M_2} \end{aligned}\]

The resonant angular frequency \(\omega[\mathrm{rad/s}]\) and damping ratio \(\zeta\) are calculated as:

\[\begin{aligned} \omega &= \sqrt{\frac{K_\mathrm{reso}(M_1+M_2)}{M_1 M_2}}\\ \zeta &= \frac{C_\mathrm{reso}}{2}\sqrt{\frac{M_1+M_2}{K_\mathrm{reso} M_1 M_2}} \end{aligned}\]

4.0.1 Measurement during Excitation

As with the standalone motor, here is a video showing the data acquisition process during Chirp signal excitation with the load mechanism attached, using the YT chart function of TwinCAT3 Scope View (TE13xx).

4.0.2 Analysis Software

As with the standalone motor, the following Python module is used to perform system identification using actual data obtained from Chirp excitation. The dataset used is freq_resp_2mass_20250902.csv, which can be found HERE.

Run_sysid_2mass_exp.py

# Copyright (c) 2021 Koichi Sakata

import numpy as np
from scipy import signal
from control import matlab
from pylib_sakata import ctrl
from pylib_sakata import plot
from pylib_sakata import fft
from pylib_sakata import meas

print('Start simulation!')

# Common parameters
Ts = 1/8000
dataNum = 10000
freqrange = [1, 1000]
freq = np.logspace(np.log10(freqrange[0]), np.log10(freqrange[1]), dataNum, base=10)
s = ctrl.tf([1, 0], [1])
z = ctrl.tf([1, 0], [1], Ts)
print('Common parameters were set.')

# Plant model
M = 0.12
C = 0.7
K = 0.0
fanti = 168
freso = 380
Creso = 3.5
M1 = (fanti / freso) ** 2 * M
M2 = M - M1
Kreso = (2.0 * np.pi * fanti) ** 2 * M2
k1 = M2/(M1 * (M1 + M2))
k2 = -1.0/(M1 + M2)
omegaPreso = np.sqrt(Kreso * (M1 + M2)/(M1 * M2))
zetaPreso = 0.5 * Creso*np.sqrt((M1 + M2)/(Kreso * M1 * M2))
Pmechs1 = ctrl.tf([1], [M, C, K]) + k1 * ctrl.tf([1], [1, 2*zetaPreso*omegaPreso, omegaPreso**2])
Pmechs2 = ctrl.tf([1], [M, C, K]) + k2 * ctrl.tf([1], [1, 2*zetaPreso*omegaPreso, omegaPreso**2])
numDelay, denDelay = matlab.pade(Ts*4, n=4)
Ds = ctrl.tf(numDelay, denDelay)
Dz = z**-3
Pns1 = Pmechs1 * Ds
Pns2 = Pmechs2 * Ds
Pnz1 = ctrl.c2d(Pmechs1, Ts, method='zoh') * Dz
Pnz2 = ctrl.c2d(Pmechs2, Ts, method='zoh') * Dz
Pnz1_frd = ctrl.sys2frd(Pnz1, freq)
Pnz2_frd = ctrl.sys2frd(Pnz2, freq)
Pnz = Pnz1
Pnz_frd = Pnz1_frd
print('Plant model was set.')

print('Getting measurement data...')
measfileName = 'data/freq_resp_2mass_20250902.csv'
# Frequency response
Pmeas_frd, coh = meas.measdata2frd(measfileName, 'ServoOutN[0]', 'ActPosUm[0]', 'FlagInject', freq, 1., 1.e-6, 8, 0.8)

# Time response
measdata = meas.getdata(measfileName)
time = measdata.time
NoiseOut = measdata.value[meas.getdataindex(measdata, 'NoiseOut[0]')]
ServoOutN = measdata.value[meas.getdataindex(measdata, 'ServoOutN[0]')]
ActPosUm = measdata.value[meas.getdataindex(measdata, 'ActPosUm[0]')]
FlagInject = measdata.value[meas.getdataindex(measdata, 'FlagInject')]

print('Plotting figures...')
# Time response
fig = plot.makefig()
ax1 = fig.add_subplot(411)
ax2 = fig.add_subplot(412)
ax3 = fig.add_subplot(413)
ax4 = fig.add_subplot(414)
plot.plot_xy(ax1, time, NoiseOut, '-', 'b', 0.5, 1.0, ylabel='[Nm]', legend='Sig', loc='upper right', title='Time response')
plot.plot_xy(ax2, time, ServoOutN, '-', 'b', 0.5, 1.0, ylabel='[Nm]', legend='In', loc='upper right')
plot.plot_xy(ax3, time, ActPosUm, '-', 'b', 0.5, 1.0, ylabel='[um]', legend='Out', loc='upper right')
plot.plot_xy(ax4, time, FlagInject, '-', 'b', 0.5, 1.0, xlabel='Time [s]', ylabel='[.]', legend='Flag', loc='upper right')

# Plant
fig = plot.makefig()
ax_mag = fig.add_subplot(311)
ax_phase = fig.add_subplot(312)
ax_coh = fig.add_subplot(313)
plot.plot_tffrd(ax_mag, ax_phase, Pmeas_frd, '-', 'm', 1.5, 1.0, ax_coh=ax_coh, coh=coh, title='Frequency response of plant')
plot.plot_tffrd(ax_mag, ax_phase, Pnz_frd, '--', 'b', 1.5, 1.0, freqrange, legend=['Measurement', 'Model'])

print('Finished.')

Start simulation!
Common parameters were set.
Plant model was set.
Getting measurement data...
Plotting figures...
Finished.

4.0.3 Tips for System Identification in a Two-Inertia System

In the system identification of a standalone motor, the process is straightforward as it only involves fitting two parameters: mass (inertia) \(M\) and damping \(C\). In contrast, a two-inertia system has five parameters to fit, which can be overwhelming. However, there is a specific strategy for efficient identification.

First, by aligning the mass line of the rigid-body model \(P_{\mathrm{rigid}}(s)\), the total mass \(M_1 + M_2\) is determined. Next, identify the resonant frequency \(f_{\mathrm{reso}}[\mathrm{Hz}]\) and the anti-resonant frequency \(f_{\mathrm{anti}}[\mathrm{Hz}]\) from the plot.

Using the following equations:

\[\begin{aligned} f_{\mathrm{reso}} &= \frac{1}{2\pi}\sqrt{\frac{K_\mathrm{reso}(M_1+M_2)}{M_1 M_2}}\\ f_{\mathrm{anti}} &= \frac{1}{2\pi}\sqrt{\frac{K_\mathrm{reso}}{M_2}} \end{aligned}\]

The relationship can be expressed as:

\[f_{\mathrm{reso}} = f_{\mathrm{anti}}\sqrt{\frac{(M_1+M_2)}{M_1}}\]

By using the previously fitted total mass \(M_1 + M_2\) and the values of \(f_{\mathrm{reso}}\) and \(f_{\mathrm{anti}}\) read from the graph, the driving-side mass \(M_1\) is determined.

Subsequently, the load-side mass \(M_2\) and the spring constant \(K_\mathrm{reso}\) are derived from the resonance or anti-resonance equations.

Finally, adjust the resonance damping \(C_{\mathrm{reso}}\) to match the height of the gain peak at \(f_{\mathrm{reso}}\). The rigid-body damping \(C\) is adjusted in the same manner as with the standalone motor.

By following this logical sequence, even a complex two-inertia system can be identified accurately.