Association of generic controllers with active magnetic bearings

[ArthurIasbeck]

Active Magnetic Bearings (AMBs) operate through a control system responsible for determining the magnitude of the electric currents applied to their actuators’ coils. These currents define the intensity of the forces generated by the AMB. Consequently, the dynamics of this type of bearing are directly influenced by the architecture of the controller associated with it.

Currently, ROSS enables the implementation of a simple PID controller. The user defines the constants $K_p$, $K_i$, and $K_d$, and the control currents are computed when the run_time_response() method is invoked, as illustrated in the structure below (provided that the Newmark integration method is employed).

However, real AMBs can hardly be stabilized using only a PID controller, since their dynamics are complex and, under any operating condition, they are subject to external disturbances that cannot be predicted.

Thus, to ensure that simulations performed in ROSS reproduce the behavior of real systems more faithfully, the MagneticBearingElement class was adjusted to allow the user to employ not only the PID controller but any controller defined from a transfer function.

It was decided to represent the controller in the state space. In summary, the user provides the transfer function $C(s)$ in the domain of the complex variable $s$, which represents the controller. This function is then discretized to obtain $C(z)$ and subsequently converted into its state-space representation,

$$ x_e(k+1) = A_e x_e(k) + B_e e(k) $$ $$ u(k) = C_e x_e(k) + D_e e(k) $$

in which $x_e$ represents the internal states of the controller, $e$ corresponds to the error—defined as the difference between the desired position of the shaft supported by the AMB and its actual position—and $u$ denotes the control action, that is, the electric current applied to the actuator.

It is noted that, through this representation, the control action at the current instant $k$ can be directly obtained from the error at the same instant, provided that the internal states of the controller are updated at each iteration.

To ensure compatibility with previously developed scripts, it is still possible to instantiate the MagneticBearingElement class using the parameters kp_pid, ki_pid, and kd_pid. In this case, the code automatically implements a simple PID controller, whose transfer function is presented below,

$$ C(s) = K_p + \frac{K_i}{s} + K_d \frac{N_f s }{N_f s + 1} $$

The figure below presents a comparison between the results obtained before and after the update of the MagneticBearingElement class. In both cases, the class was instantiated using only the parameters kp_pid, ki_pid, and kd_pid.

With the aim of validating the implementation of a generic controller, the use of a lead controller in cascade with the PID controller was proposed, as indicated below.

$$ C_{\text{lead}}(s) = K \frac{1 + Ts}{1 + \alpha Ts} $$

To evaluate the effect of including the lead controller, the run_time_response method was used to obtain the shaft displacement after the application of an impulsive force to the rotor disk. The rotor employed in this validation is presented below, followed by the code used for its construction.

def build_simple_rotor():
    # Material
    steel = rs.Material(name="Steel", rho=7850, E=211e9, G_s=81.2e9)

    # Shaft elements
    L = 0.2 
    i_d = 0.0
    o_d = 0.05
    shaft_elements = [
        rs.ShaftElement(
            L=L,
            idl=i_d,
            odl=o_d,
            material=steel,
            shear_effects=True,
            rotary_inertia=True,
            gyroscopic=True,
            alpha=5.0,  
            beta=1e-4,  
        )
        for _ in range(4)
    ]
		
    # Disk
    disk = rs.DiskElement.from_geometry(
        n=2, material=steel, width=0.05, i_d=0.0, o_d=0.07
    )

    # Active magnetic bearings at nodes 1 and 5, with electromagnetic parameters
    # Example values (adjust according to your real system)
    g0 = 1e-3  # air gap in meters
    i0 = 1.0  # bias current in A
    ag = 1e-4  # pole area in m²
    nw = 200  # number of turns
    alpha = np.deg2rad(45)  # pole angle in radians
    kp = 100.0  # PID proportional gain
    ki = 0  # PID integral gain
    kd = 10.0  # PID derivative gain
    n_f = 10_000  # derivative filter cutoff frequency
    k_amp = 1.0  # amplifier gain
    k_sense = 1.0  # sensor gain

    s = rs.MagneticBearingElement.s
    pid_controller = kp + ki / s + kd * s * (1 / (1 + (1 / n_f) * s))

    k_lead = 1
    T_lead = 0.5
    alpha_lead = 0.1
    lead_controller = k_lead * (T_lead * s + 1) / (alpha_lead * T_lead * s + 1)

    controller_transfer_function = pid_controller * lead_controller

    mb1 = rs.MagneticBearingElement(
        n=0,
        g0=g0,
        i0=i0,
        ag=ag,
        nw=nw,
        alpha=alpha,
        controller_transfer_function=controller_transfer_function,
        k_amp=k_amp,
        k_sense=k_sense,
        tag="AMB_left",
    )

    mb2 = rs.MagneticBearingElement(
        n=4,
        g0=g0,
        i0=i0,
        ag=ag,
        nw=nw,
        alpha=alpha,
        controller_transfer_function=controller_transfer_function,
        k_amp=k_amp,
        k_sense=k_sense,
        tag="AMB_right",
    )

    # Rotor assembly
    rotor = rs.Rotor(
        shaft_elements=shaft_elements, disk_elements=[disk], bearing_elements=[mb1, mb2]
    )

    return rotor

The figure below presents the displacements obtained from the execution of the run_time_response method before and after the implementation of the lead controller, whose transfer function is shown next. As expected, the inclusion of the lead controller improves the steady-state performance of the system.

k_lead = 1
T_lead = 0.5
alpha_lead = 0.1
lead_controller = k_lead * (T_lead * s + 1) / (alpha_lead * T_lead * s + 1)

Finally, it was necessary to reformulate the calculation of the equivalent damping and stiffness coefficients of the AMB, as indicated in [1]:

$$ k_e = k_z + k_i , \text{Re} { C(j\omega) } $$

$$ c_e = \frac{k_i}{\omega} , \text{Im} { C(j\omega) } $$

Once the equivalent stiffness ($k_{eq}$) and damping ($c_{eq}$) coefficients have been computed, it is necessary to consider that these values are associated with the axes along which the electromagnetic actuators operate. These axes, denoted $x'y'$, differ from the $xy$ axes used in ROSS, as illustrated in the figure below. Therefore, a conversion between the two coordinate systems must be performed.

The stiffness matrix $K$ associated with the $xy$ axes can be defined as follows:

$$ \begin{bmatrix} F_x \\ F_y \end{bmatrix} = K \begin{bmatrix} x \\ y \end{bmatrix} $$

where

$$ K = \begin{bmatrix} k_{xx} & k_{xy} \\ k_{yx} & k_{yy} \end{bmatrix} $$

By employing the rotation matrix $R$:

$$ R(\alpha) = \begin{bmatrix} \cos\alpha & \sin\alpha \\ -\sin\alpha & \cos\alpha \end{bmatrix} $$

it follows that

$$ \begin{aligned} x' &= R(\alpha), x \\ F' &= R(\alpha), F \end{aligned} $$

Based on the definitions presented thus far, it is possible to conclude that

$$ K = R^{-1}(\alpha), K', R(\alpha) $$

Since

$$ K' = \begin{bmatrix} k_{eq} & 0 \\ 0 & k_{eq} \end{bmatrix} $$

the stiffness matrix $K$ can be determined from its counterpart $K'$. The same procedure applies to the damping matrix $C$. Thus, after computing the equivalent stiffness ($k_{eq}$) and damping ($c_{eq}$) coefficients, the matrices $K'$ and $C'$ are defined so that, based on the relations presented thus far, the final matrices $K$ and $C$ can be obtained.

Modifications implemented in the code

The implementation of the modifications listed below was requested through the opening of Pull Request #1237.

1. Dependencies

• Added control library: The requirements.txt has been updated to include the control package, which is now used for transfer function operations, frequency response analysis, controller synthesis, and system discretization.

2. MagneticBearingElement Refactoring (ross/bearing_seal_element.py)

• General Transfer Functions: Users can now pass a custom control.TransferFunction object via the controller_transfer_function parameter. If not provided, a PID controller (with a derivative filter) is constructed automatically.
• Frequency-Dependent Coefficients: The element now calculates equivalent stiffness ($K_{eq}$) and damping ($C_{eq}$) matrices as a function of frequency based on the closed-loop response of the controller, rather than using static PID gains.
• Coordinate Transformation: Added the sensors_axis_rotation parameter. The class now computes the rotation matrices to transform the control forces from the sensor frame ($v, w$) to the rotor frame ($x, y$), resulting in fully populated stiffness and damping matrices (including cross-coupling terms) when axes are rotated.
• Controller Discretization: Added get_analog_controller() and build_controller(dt). The latter discretizes the continuous transfer function using the Tustin method to generate State-Space matrices ($A_c, B_c, C_c, D_c$) for time-domain simulations.
• Updated compute_pid_amb: Refactored to calculate forces using the discrete state-space matrices, updating internal controller states ($x_c$) at each time step.

3. Rotor Assembly & Time Response (ross/rotor_assembly.py)

• Integration with Time Response: The magnetic_bearing_controller method in the Rotor class was updated to handle the coordinate transformation logic (projecting displacements to sensor axes and forces back to rotor axes).
• Initialization: The _init_ambs_for_integrate method now calls build_controller(dt) to ensure the magnetic bearings are discretized with the correct time step before simulation.

4. Controller Utilities Module (controller.py)

• New Module: Introduced a dedicated controller utilities module to standardize controller definition and analysis across the codebase.
• Controller Builders: Implemented reusable helper functions for common control structures, including:
  • pid (PID with filtered derivative),
  • lead_lag,
  • second_order,
  • low_pass_filter,
  • notch_filter,
  • lqg (LQR + Kalman filter, returned as a transfer function).
• Controller Composition: Added combine(*args) to support clean and explicit series (cascade) connections of multiple transfer functions.
• Frequency-Domain Visualization: Added plot_frequency_response(...), a Plotly-based utility for comparing magnitude and phase responses of one or more controllers or filters.
• All functions return control.TransferFunction objects (or compatible LTI systems), ensuring seamless compatibility with the MagneticBearingElement API and frequency-domain analyses.

5. Utilities (ross/utils.py)

• Added helper validation functions:
  • is_scalar_or_list
  • is_transfer_function_or_none
  • is_list_or_none
  • is_scalar

6. Documentation

• New Section 6.1 – Auxiliary Methods for Defining and Evaluating Transfer Functions:

  Added a dedicated subsection describing the controller utilities implemented in controller.py.

  This section documents controller builder functions, composition patterns, expected return types, and frequency-response plotting utilities, providing a clear and reusable workflow for defining complex AMB controllers compatible with the new generalized control architecture.

7. Tests

• New Test Case: Added test_magnetic_bearing_with_lead_controller_matches_frequency_response in test_bearing_seal_element.py to verify that the calculated $K(\omega)$ and $C(\omega)$ match the analytical expectation of a Lead controller.
• Updated Existing Tests: Updated test_magnetic_bearing_element and test_amb_controller (in test_rotor_assembly.py) to align with the new implementation and numerical precision.

References

[1] Xu, Yuanping, Jin Zhou, e Chaowu Jin. “Identification of Dynamic Stiffness and Damping in Active Magnetic Bearings Using Transfer Functions of Electrical Control System”. Journal of Mechanical Science and Technology 33, n. 2 (2019): 571–77. https://doi.org/10.1007/s12206-019-0110-y.