Estimation of the force-displacement transfer function for Magnetic Bearings

As we can see in the excerpt below, the API (American Petroleum Institute) STANDARD 617 requires that the bearing transfer functions relating force and displacement be provided. The main objective of the changes proposed in this issue is to enable ROSS to provide the user with such transfer functions.

As can be seen in the diagram below, the PID controller produces a current signal that is injected into the electromagnetic actuator which, in turn, produces a force that is applied to the shaft and generates a displacement in the bearing plane. To determine the transfer function that relates such displacement to the magnetic force produced by the actuator, the least squares method was used. Since the relationship between force and displacement is usually represented by a second-order transfer function, given that most mechanical systems can be modeled using Newton's second law, it was considered that the relationship between $F$ and $y$ can be represented by the discrete model $G(z) = (b_0 z + b_1)/(z^2 + a_0 z + a_1)$.

The transfer function $G(z)$ can be represented in the time domain by the following difference equation,

$$ y(k) = -a_0 y(k-1) - a_1 y(k-2) + b_0 u(k-1) + b_1 u(k-2) $$

in which $y(k)$ and $u(k)$ represent the system output and input at instant $k$, respectively.

The difference equation in question can be rewritten based on the definition of $\varphi$ and $\theta$, as follows,

$$ y(k) = \left[\begin{matrix} y(k-1) & y(k-2) & u(k-1) & u(k - 2)\end{matrix} \right] \left[ \begin{matrix} -a_0 \ -a_1 \ b_0 \ b_1 \end{matrix} \right]^T = \varphi^T(k-1) \theta $$

The least squares method is based on minimizing the difference between the actual output $y_r$, in this case produced by ROSS, and the estimated output $y$ produced by the model. The residue to be minimized can be computed based on the following equation.

$$ V(\theta,t)=\frac12\sum_{i=1}^t\left(y_r(i)-\varphi^T(i) \theta\right)^2 $$

Assuming that

$$ \Phi = \left[ \begin{matrix} \varphi^T(1) \ \dots \ \varphi^T(k) \end{matrix} \right]^T, \quad Y = \left[ \begin{matrix} y(1) \ \dots \ y(k) \end{matrix} \right]^T $$

it is possible to analytically determine the values of $\theta$ that minimize $V$ by using the equation introduced below, the demonstration of which can be found in Astrom, Karl Johan, and Bjorn Wittenmark. Adaptive Control. 2nd ed. USA: Addison-Wesley Longman Publishing Co., Inc., 1994.

$$ \theta = (\Phi^T \Phi)^{-1} \Phi^T Y $$

In order for the equation above to be used to determine the parameters $\theta$ associated with $G(z)$, several values of magnetic force and displacement were produced by executing the rotor.run_time_response() method. $y_w$ disturbances were generated at the system output by applying external forces to the rotor disc (associated with node 27). The version of ROSS used allows obtaining the magnetic force produced by the AMB and is available at https://github.com/mcarolalb/ross.git@dev/magnetic-controller. The rotor model used to generate the results presented below is shown in the next figure.

The graph below shows the temporal evolution of the magnitude of the force applied to the rotor disk, responsible for producing disturbances in the output to be measured.

The current signal $i$ produced by the controller, the magnetic force $F$ and the displacement $y$ of the shaft in the AMB plane are shown in the graph below.

The application of the least squares method resulted in the following model,

$$ \frac{Y(z)}{F(z)}=\frac{-0.00008067z+0.00007740}{z^2-1.4500z+0.4560}, \quad dt = 0.005 \text{ sec} $$

The graph below allows a visual comparison between the output produced by ROSS $y_r$ and the output $y$ obtained by using the model introduced above.

It is noted that the curves practically overlap, which indicates an excellent fit of the model to the data produced by ROSS. In fact, the fit of the model can be computed as follows.

$$ fit=1-\frac{||y-y_{est}||}{||y-\overline{y}||}=0.9443 $$

The method used to obtain the above result can be seen below:

def estimate_second_order_model(U, Y):
    """
    Method for estimating a model based on applied inputs and collected outputs.
    The model used was: (b_0 z + b_1) / (z^2 + a_0 z + a_1)

    :param U: One-dimensional array containing the inputs applied over time.
    :param Y: One-dimensional array containing the outputs measured over time.
    :return: A tuple containing the coefficients of the estimated model (a_0, a_1, b_0, b_1)
    """
    y_1 = 0
    y_2 = 0
    u_1 = 0
    u_2 = 0

    Phi = []
    for k in range(0, U.shape[0]):
        phi_T = [-y_1, -y_2, u_1, u_2]
        Phi.append(phi_T)
        y_2 = y_1
        y_1 = Y[k]
        u_2 = u_1
        u_1 = U[k]

    Y = np.matrix(Y).T
    Phi = np.matrix(Phi)

    theta = (Phi.T * Phi).I * Phi.T * Y

    a_0 = theta[0, 0]
    a_1 = theta[1, 0]
    b_0 = theta[2, 0]
    b_1 = theta[3, 0]

    Y_est = []
    y_1 = 0
    y_2 = 0
    u_1 = 0
    u_2 = 0
    for k in range(0, U.shape[0]):
        phi_T = [-y_1, -y_2, u_1, u_2]
        y_est = (phi_T * theta)[0, 0]
        Y_est.append(y_est)
        y_2 = y_1
        y_1 = y_est
        u_2 = u_1
        u_1 = U[k]

    return a_0, a_1, b_0, b_1, Y_est

The estimation of the transfer functions that relate force and displacement is still a work in progress. Some of the activities yet to be carried out are listed below.

Evaluate the result of the estimate for other input types
Evaluate the quality of the estimate in the frequency domain
Propose a way to include the computation of the $F/x$ transfer function in ROSS
Allow the user to adjust the model (through the number of poles and zeros)
Compute the $i/x$ transfer function (actually used in the control system design)

petrobras / ross

Estimation of the force-displacement transfer function for Magnetic Bearings #1115