Tutorial - Active Magnetic Bearings (AMB)

Tutorial - Active Magnetic Bearings (AMB)#

This tutorial provides a comprehensive guide on how to model, simulate, and analyze Active Magnetic Bearings (AMBs) using ROSS.

Section 1: Introduction and Working Principle#

Active Magnetic Bearings (AMBs) are support mechanisms that levitate a rotating shaft using electromagnetic forces, eliminating physical contact between the rotor and the stator.

1.1 Working Principle#

An AMB system operates as a closed-loop feedback control system consisting of three main components:

Sensors: Measure the radial displacement of the rotor ($x$).
Controller: Processes the displacement signal and calculates the necessary current correction ($i$) to maintain the rotor’s position (setpoint).
Actuators (Electromagnets): Receive the current and generate the magnetic force required to pull the rotor back to the center.

1.2 Mathematical Model#

In ROSS, the electromagnetic force $F$ is linearized around a nominal operating point (bias current $i_0$ and nominal gap $g_0$). The force equation is generally described by:

\[ F = K_i \; i + K_s \; x \]

Where:

$K_i$ is the current stiffness (Force/Current factor).
$K_s$ is the negative electromagnetic stiffness (Force/Displacement factor), representing the inherent instability of the magnetic attraction.

Section 2: Declaring a Magnetic Bearing in ROSS#

To use an AMB in ROSS, you must instantiate the MagneticBearingElement class. You need to provide physical parameters of the electromagnet and the gains for the PID controller.

import ross as rs
import numpy as np

# Physical parameters of the actuator
n_node = 4  # Node where the bearing is located
g0 = 1e-3  # Nominal air gap (m)
i0 = 1.0  # Bias current (A)
ag = 1e-4  # Pole area (m²)
nw = 200  # Number of windings
alpha = np.pi / 8  # Pole angle (rad)

# PID Controller gains
Kp = 1500  # Proportional gain
Kd = 10  # Derivative gain
Ki = 100  # Integral gain

# Instantiation
amb = rs.MagneticBearingElement(
    n=n_node,
    g0=g0,
    i0=i0,
    ag=ag,
    nw=nw,
    alpha=alpha,
    kp_pid=Kp,
    kd_pid=Kd,
    ki_pid=Ki,
    tag="AMB_1",
)

Section 3: Evaluating Time Response#

Active Magnetic Bearing (AMB) systems are inherently coupled, closed-loop systems where the control forces directly depend on the dynamic states of the rotor. To simulate this behavior accurately and efficiently, ROSS utilizes the AmbTimeResponse class.

Instead of relying solely on step-by-step numerical integration at every time step, this class constructs a comprehensive, unified State-Space Model in the modal domain. This global model fully encapsulates the rotor structural dynamics, the actuator magnetic forces, and the feedback controller transfer functions.

3.1 State-Space Formulation#

The coupled system is formulated by combining three fundamental sets of equations:

Rotor Dynamics (Physical Domain): $$M \ddot{q}(t) + C(\Omega) \dot{q}(t) + K q(t) = F_{ext}(t) + F_{amb}(t)$$
Magnetic Actuator Force Linearization: $$F_{amb}(t) = K_s \cdot q_{amb}(t) + K_i \cdot i(t)$$
Controller Dynamics (State-Space Representation): $$\dot{x}_c(t) = A_c x_c(t) + B_c e(t)$$i(t) = C_c x_c(t) + D_c e(t)$$

By projecting the rotor equations into the modal domain using the modal matrix $\Phi$ ($q = \Phi \eta$) and coupling them with the controller states, ROSS builds a single consolidated closed-loop state vector $z = [\eta^T, \dot{\eta}^T, x_c^T]^T$. The entire system is described by the following linear state-space equations:

\[\dot{z}(t) = A_{cl} z(t) + B_{cl} u(t)\]

\[y(t) = C_{cl} z(t) + D_{cl} u(t)\]

The global block-structured closed-loop matrices are defined as follows:

\[\begin{split}A_{cl} = \begin{bmatrix} 0_{m \times m} & I_{m \times m} & 0_{m \times n} \\ M_m^{-1}\Phi^T\left(\phi R^T (K_x - K_i D_c)R\phi^T - K\right)\Phi & -M_m^{-1}C_m & M_m^{-1}\Phi^T\phi R^T K_i C_c \\ -B_c R \phi^T \Phi & 0_{n \times m} & A_c \end{bmatrix}\end{split}\]

\[\begin{split}B_{cl} = \begin{bmatrix} 0_{m \times n_c} & 0_{m \times N} \\ 0_{m \times n_c} & M_m^{-1}\Phi^T \\ -B_c R \phi^T \phi R^T & 0_{n \times N} \end{bmatrix}\end{split}\]

\[\begin{split}C_{cl} = \begin{bmatrix} \Phi & 0_{N \times m} & 0_{N \times n} \\ \phi^T \Phi & 0_{n_c \times m} & 0_{n_c \times n} \\ R \phi^T \Phi & 0_{n_c \times m} & 0_{n_c \times n} \\ R^T(K_x - K_i D_c)R \phi^T \Phi & 0_{n_c \times m} & R^T K_i C_c \\ (K_x - K_i D_c)R \phi^T \Phi & 0_{n_c \times m} & K_i C_c \\ -D_c R \phi^T \Phi & 0_{n_c \times m} & C_c \end{bmatrix}\end{split}\]

\[\begin{split}D_{cl} = \begin{bmatrix} 0_{N \times n_c} & 0_{N \times N} \\ 0_{n_c \times n_c} & 0_{n_c \times N} \\ 0_{n_c \times n_c} & 0_{n_c \times N} \\ 0_{n_c \times n_c} & 0_{n_c \times N} \\ 0_{n_c \times n_c} & 0_{n_c \times N} \\ 0_{n_c \times n_c} & 0_{n_c \times N} \end{bmatrix}\end{split}\]

Where:

$m$ represents the number of dynamic modes retained after modal reduction.
$n$ indicates the total dimension of the controller’s internal state vector.
$N$ determines the total number of physical degrees of freedom (DoF) of the rotor.
$n_c$ expresses the number of coordinated control channels for the magnetic bearings (equal to twice the number of magnetic bearings on the rotor).
$M_m$ and $C_m$ are the modal mass and damping matrices, calculated respectively by $\Phi^T M \Phi$ and $\Phi^T C \Phi$.
$K$ is the structural physical stiffness matrix of the rotor.
$\phi$ defines the geometric coupling matrix that locates the bearing nodes within the physical degrees of freedom of the global system.
$R$ establishes the rotation matrix corresponding to the tilt and angular orientation of the sensor and actuator coordinate axes.
$K_x$ and $K_i$ represent diagonal matrices built from the linear position stiffness ($K_s$) and current stiffness ($K_i$) coefficients of the magnetic forces.
$A_c, B_c, C_c, D_c$ constitute the dynamic state-space matrices associated with the transfer function of the designed controller.

3.2 The `weight` Parameter#

When running a time response simulation via rotor.run_time_response(...), you can specify the option weight=True.

How it works: Setting weight=True automatically calculates and adds the static gravitational forces ($M \cdot g$) acting on the rotor elements to the external force vector ($F_{ext}$).
Engineering relevance: This is crucial for AMB systems, as it simulates the static load that the magnetic bearings must overcome to levitate the shaft. It allows you to observe the steady-state bias currents and position offsets required to maintain the rotor centered under gravity.

Note on Backwards Compatibility: The legacy step-by-step numerical integration via the Newmark method can still be explicitly invoked by setting method="newmark". Legacy scripts and workflows utilizing this method remain fully supported and functional.

import ross as rs
import numpy as np

# 1. Create the rotor with the AMB element
rotor = rs.rotor_example_amb_complex_controllers()
rotor.plot_rotor(nodes=999).show()

# Define simulation parameters
speed = 500.0  # rad/s
t = np.arange(0, 10, 0.001)  # Time vector
F = np.zeros((len(t), rotor.ndof))  # External force vector

# 2. Run time response using the new state-space default approach with gravity included
response = rotor.run_time_response(speed, F, t, weight=True)

# 3. Plot time response for a physical node using a Probe
obs_node = 28
probe_y = rs.Probe(node=obs_node, angle=np.pi / 2)
fig = response.plot_1d(probe=[probe_y])
fig.show()

# 4. Use the specialized automatic AMB plotting tool for displacements
response.plot_amb_disps().show()

Section 4: Collecting Forces and Control Currents#

The run_time_response method returns an instance of the AmbTimeResponseResults class. This dedicated object acts as a central repository for all internal tracking data related to the magnetic bearings.

4.1 Automated Plotting Methods#

The results object provides high-level automated plotting methods designed to visualize key control loop variables across all bearings simultaneously:

response.plot_amb_disps(): Plots the displacement profiles.
response.plot_amb_currents(): Plots the control currents over time.
response.plot_amb_forces(): Plots the total magnetic forces applied to the rotor.

4.2 Accessing Raw Arrays and Matrix Ordering Rules#

If you need to perform custom post-processing, optimization, or alternative plotting configurations, you can access the raw NumPy arrays stored as attributes inside the results object:

response.x_amb: Bearing displacements matrix.
response.I: Control currents matrix.
response.F_x: Total magnetic bearing forces matrix along the coordinate axes.

Data Structure and Indexing Conventions:#

When inspecting these raw matrices, the columns follow a strict geometric and layout convention:

Axis Arrangement: Within the data columns for any given bearing, the $x$-axis data always precedes the $y$-axis data. For instance, column 2*i holds the $x$ response, and column 2*i + 1 holds the $y$ response.
Node Ordering: Bearings appear along the columns in the exact spatial sequence they occupy on the rotor shaft. Bearings positioned at lower node numbers appear first (from left to right in the matrix columns).

For example, in a system with two magnetic bearings located at Node 12 and Node 43:

Column 0: Node 12 ($x$-axis)
Column 1: Node 12 ($y$-axis)
Column 2: Node 43 ($x$-axis)
Column 3: Node 43 ($y$-axis)

# 1. Use the automatic plotting methods built into AmbTimeResponseResults
response.plot_amb_currents().show()
response.plot_amb_forces().show()

# 2. Extract and inspect raw data vectors directly
raw_currents = response.I
raw_forces = response.F_x
raw_displacements = response.x_amb

print("Shape of the raw current matrix (Time steps x Channels):", raw_currents.shape)
print(
    "\nFirst 5 time steps of raw currents for all AMB axes (ordered by node, x before y):"
)
print(raw_currents[:5, :])

print("\nFirst 5 time steps of raw magnetic forces:")
print(raw_forces[:5, :])

Section 5: Sensitivity Analysis (ISO 14839)#

5.1 Definition and ISO Standard#

The sensitivity function $S(j\omega)$ determines the robustness of the control system. According to ISO 14839-3, the sensitivity is the ratio of the sensor output to a disturbance added to the sensor signal. It represents how much a disturbance is amplified by the feedback loop.

\[ S(j\omega) = \frac{1}{1 + G(j\omega)H(j\omega)} \]

The ISO standard classifies the stability margin based on the peak sensitivity value ($S_{max}$):

Zone A ($S_{max} < 3.0$): Unrestricted operation (approx 9.5 dB).
Zone B ($3.0 < S_{max} < 4.0$): Restricted operation.
Zone C ($S_{max} > 4.0$): Evaluation required, potential instability.

5.2 Simulation Procedure in ROSS#

To calculate the sensitivity function, ROSS performs a specific time-domain simulation:

A Logarithmic Chirp Signal (sweeping sine wave) is injected as a disturbance into the control loop at the AMB sensor location.
The system response (displacement) is recorded.
A Frequency Response Function (FRF) is computed between the Output (Disturbed Signal) and the Input (Excitation Signal).

# Run sensitivity analysis
# This automatically performs the chirp injection and FFT computation
rotor = rs.rotor_example_amb_simple()
sensitivity_results = rotor.run_amb_sensitivity(
    speed=0,
    t_max=10,  # Duration must be long enough for the chirp
    dt=1e-5,  # Time step
    disturbance_amplitude=1e-5,
    disturbance_min_frequency=0.001,  # Hz
    disturbance_max_frequency=150,
    amb_tags=["AMB_1"],
    verbose=0,
);

5.3 The `SensitivityResults` Object#

The run_amb_sensitivity method returns a SensitivityResults object. Below is a summary of its attributes and methods.

Attributes Table#

Attribute	Type	Description
`sensitivities`	`dict`	Dictionary containing the complex Sensitivity FRF ($S(j\omega)$) for each AMB and axis ($x$, $y$).
`sensitivities_abs`	`dict`	The magnitude (absolute value) of the sensitivity function.
`sensitivities_phase`	`dict`	The phase angle of the sensitivity function.
`sensitivities_frequencies`	`np.ndarray`	The frequency vector (in Hz) corresponding to the sensitivity arrays.
`max_abs_sensitivities`	`dict`	The peak sensitivity value ($S_{max}$) for each AMB/axis. Used for ISO classification.
`sensitivity_run_time_results`	`dict`	Contains raw time-domain arrays: `excitation_signal` (chirp), `disturbed_signal` (input to controller), and `sensor_signal`.
`sensitivity_compute_dofs`	`dict`	Mapping of AMB tags to their corresponding Degree of Freedom indices.

Accessing Results Data#

The attributes within the SensitivityResults object (such as max_abs_sensitivities, sensitivities, sensitivities_abs, and sensitivities_phase) are organized as nested dictionaries. This structure allows you to easily retrieve data for a specific magnetic bearing and a specific axis.

Structure Hierarchy:

First Level Key: The AMB tag (string), exactly as defined when creating the MagneticBearingElement.
Second Level Key: The axis (string), either "x" or "y".
Value: The corresponding data (scalar for max sensitivity, or array for FRFs).

Visual Representation:

sensitivity_results.max_abs_sensitivities = {
    "AMB_Tag_1": {
        "x": <scalar_value>,
        "y": <scalar_value>
    },
    "AMB_Tag_2": {
        "x": <scalar_value>,
        "y": <scalar_value>
    }
    # ...
}

Example: Retrieving Maximum Sensitivity

The following example demonstrates how to programmatically access the maximum absolute sensitivity ($S_{max}$) for a specific axis of a specific bearing to check against ISO limits.

# Assume 'sensitivity_results' is the object returned by rotor.run_amb_sensitivity()

# 1. Define the target Bearing Tag and Axis
target_amb_tag = "AMB_1"  # Must match the 'tag' used in MagneticBearingElement
target_axis = "x"

# 2. Access the nested dictionary
s_max = sensitivity_results.max_abs_sensitivities[target_amb_tag][target_axis]

print(f"Peak Sensitivity for {target_amb_tag} ({target_axis}-axis): {s_max:.4f}")

# 3. Check against ISO 14839-3 Zone A limit
if s_max < 3.0:
    print("Result: Zone A (Unrestricted Operation)")
else:
    print("Result: Zone B or C (Evaluation required)")

Plotting Methods#

.plot()

Displays the Bode plot (Magnitude and Phase) of the sensitivity function. This is the primary tool to check compliance with ISO 14839.

# Plot Bode diagram of Sensitivity
fig_bode = sensitivity_results.plot(
    frequency_units="Hz",
    magnitude_scale="decibel",  # Useful to check margins in dB
    xaxis_type="log",
)
fig_bode.show()

.plot_time_results()

Displays the raw time-domain signals used to calculate the sensitivity. This is useful for debugging to ensure the chirp signal was applied correctly and the system remained stable during the test.

# Plot the Chirp injection and system response
fig_time = sensitivity_results.plot_time_results()
fig_time.show()

Section 6: Implementing Complex Controllers (Cascade)#

ROSS supports arbitrary transfer functions using the python-control library. You can combine controllers (e.g., in series/cascade) by multiplying their transfer functions.

Example: A PID controller cascaded with a Lead Compensator (to improve phase margin).

import ross as rs

# 1. Define the Laplace variable 's'
s = rs.MagneticBearingElement.s

# 2. Define PID Controller
Kp, Ki, Kd = 1500, 100, 10
# Standard PID formulation with a derivative filter (N) usually handled internally or explicitly:
# Here we use a pure PID for simplicity, or construct manual TF
TF_PID = Kp + Ki / s + Kd * s

# 3. Define Lead Controller
# Lead: (alpha * T * s + 1) / (T * s + 1), where alpha > 1
alpha = 3.0
T = 0.001
TF_Lead = (alpha * T * s + 1) / (T * s + 1)

# 4. Cascade (Series) Connection
# In frequency domain, series connection is multiplication
TF_Combined = TF_PID * TF_Lead

# 5. Instantiate AMB with the complex controller
amb_cascade = rs.MagneticBearingElement(
    n=4,
    g0=1e-3,
    i0=1.0,
    ag=1e-4,
    nw=200,
    controller_transfer_function=TF_Combined,  # Pass the object here
    tag="AMB_Cascade",
)

ROSS handles the discretization and state-space conversion of TF_Combined automatically during the simulation.

It is worth noting that, if necessary, the MagneticBearingElement can still be defined by specifying only the kp_pid, ki_pid, and kd_pid parameters.

6.1 Auxiliary Methods for Defining and Evaluating Transfer Functions#

Conventions and returned types#

All “controller builder” functions return python-control transfer functions (control.TransferFunction), except where noted:

pid(...) → TransferFunction
lead_lag(...) → TransferFunction
second_order(...) → TransferFunction
low_pass_filter(...) → TransferFunction
notch_filter(...) → TransferFunction
lqg(...) → TransferFunction (converted from state-space to TF internally)
combine(*args) → product of transfer functions (series connection)

from ross.bearings.magnetic.amb_controllers import *

Below, each of the auxiliary methods introduced above is described in detail.

pid(k_p, k_i, k_d, n_f=10000)

Builds a PID controller with a filtered derivative term.

k_p (float): proportional gain
k_i (float): integral gain
k_d (float): derivative gain
n_f (float, optional): derivative filter “corner” scaling (higher → closer to ideal derivative without filtering)

C = pid(k_p=3.0, k_i=5.0, k_d=0.02, n_f=300)

lqg(A, B, C, Q_lqr, R_lqr, Q_kalman, R_kalman)

Builds an LQG controller (LQR state feedback + Kalman filter) and returns it as a transfer function.

A, B, C: system matrices (array-like). Converted to float numpy arrays.
Q_lqr, R_lqr: LQR weighting matrices.
Q_kalman, R_kalman: process/measurement noise covariances for the Kalman filter design.

Notes / tips

Dimensions must match python-control expectations, where n, m, and p are, respectively, the number of states, controls, and outputs.:
- A is (n, n), B is (n, m), C is (p, n)
- Q_lqr is (n, n), R_lqr is (m, m)
- Q_kalman is (n, n), R_kalman is (p, p)

A = [[0, 1], [-2, -0.5]]
B = [[0], [1]]
C = [[1, 0]]

Q_lqr = np.diag([10, 1])
R_lqr = np.array([[1]])

Q_kalman = np.diag([0.1, 0.1])
R_kalman = np.array([[0.5]])

C_lqg = lqg(A, B, C, Q_lqr, R_lqr, Q_kalman, R_kalman)

lead_lag(tau, alpha, k=1.0)

Creates a first-order lead/lag compensator:

\[C(s)=\frac{k(\tau s + 1)}{\alpha \tau s + 1}\]

tau (float): time constant
alpha (float): pole/zero separation factor
- Lead is typically 0 < alpha < 1 (pole farther left than zero)
- Lag is typically alpha > 1 (pole closer to origin than zero)
k (float, optional): gain multiplier

C_lead = lead_lag(tau=0.02, alpha=0.1, k=2.0)
C_lag = lead_lag(tau=0.5, alpha=5.0, k=1.0)

second_order(b2, b1, b0, a1, a0)

Creates a generic second-order transfer function:

\[H(s)=\frac{b_2 s^2 + b_1 s + b_0}{s^2 + a_1 s + a_0}\]

H = second_order(b2=1, b1=0.2, b0=10, a1=0.5, a0=25)

low_pass_filter(w_c, k=1.0)

First-order low-pass filter:

\[F(s)=\frac{k\,\omega_c}{s + \omega_c}\]

w_c (float): cutoff frequency in rad/s
k (float, optional): DC gain factor

F = low_pass_filter(w_c=200, k=1.0)

notch_filter(w_n, zeta_z, zeta_p, k=1.0)

Second-order notch filter (zeros and poles around the same natural frequency):

\[N(s)=k\frac{s^2 + 2\zeta_z\omega_n s + \omega_n^2}{s^2 + 2\zeta_p\omega_n s + \omega_n^2}\]

w_n (float): notch center frequency (rad/s)
zeta_z (float): zero damping (controls notch depth/shape)
zeta_p (float): pole damping (controls bandwidth/sharpness)
k (float, optional): gain multiplier

N = notch_filter(w_n=120, zeta_z=0.01, zeta_p=0.2, k=1.0)

combine(*args)

Multiplies transfer functions in series.

*args: any number of control.TransferFunction objects

C = combine(
    pid(2, 1, 0.02, n_f=200),
    lead_lag(0.03, 0.1),
    notch_filter(120, 0.02, 0.2),
    low_pass_filter(400),
)

plot_frequency_response(*systems, **kwargs)

Plots magnitude (dB) and phase (degrees) for one or more systems using Plotly.

*systems: one or more LTI systems compatible with ct.frequency_response(...) (control.TransferFunction or control.StateSpace)
w_min (float, default=1e-2): minimum frequency (rad/s)
w_max (float, default=1e3): maximum frequency (rad/s)
n_points (int, default=1000): number of log-spaced frequency points
title (str, default=”Frequency Response”): plot title
legends (list[str] | None): legend labels; must match number of systems

plot_frequency_response(
    C,
    N,
    C * N,
    legends=["Controller", "Filter", "Combined"],
    w_min=1e-1,
    w_max=1e5,
    n_points=1500,
    title="Bode (Controller / Filter / Combined)",
)

Section 7: Equivalent Stiffness Calculation#

For linear frequency-domain analyses (e.g., run_modal, run_campbell, run_ucs), ROSS cannot evaluate the discrete time-domain controller states directly. Instead, the MagneticBearingElement computes frequency-dependent equivalent coefficients that represent the closed-loop AMB behavior in the form of an equivalent stiffness and damping.

With the updated implementation, the displacement sensor gain (k_sense, in V/m) and the power amplifier gain (k_amp, in V/A) are explicitly included in the loop gain used to build these equivalent coefficients.

Internally, the element performs:

Build (or retrieve) the continuous-time controller transfer function $C(s)$ (PID with derivative filter or a custom transfer function).
Evaluate its frequency response $C(j\omega)$ over a frequency grid $\omega$ (rad/s).
Split the response into real and imaginary parts: $ C(j\omega)=\Re\{C(j\omega)\} + j \, \Im\{C(j\omega)\} $
Map these parts into equivalent stiffness and damping using the electromagnetic constants:
- $K_s$: negative electromagnetic stiffness (from force linearization)
- $K_i$: current-to-force gain (from force linearization)

Equivalent stiffness#

The equivalent stiffness includes the open-loop negative stiffness $K_s$ plus the real part of the controller contribution scaled by the sensor and amplifier gains:

\[K_{eq}(\omega) = K_s + K_i , k_{amp}, k_{sense},\Re\{C(j\omega)\}\]

Equivalent damping#

The equivalent damping is formed from the imaginary part of the controller contribution, scaled by the same gains and divided by $\omega$:

\[ C_{eq}(\omega) = \frac{K_i , k_{amp}, k_{sense}}{\omega},\Im\{C(j\omega\} \]

Notes on implementation details#

In the code, $C(j\omega)$ is obtained from control.frequency_response(), stored as Hjw, and then:
- C_real = Hjw.real
- C_imag = Hjw.imag
The arrays are computed as: $$ k_{eq} = K_s + K_i,k_{amp},k_{sense},C_{real} \qquad c_{eq} = K_i,k_{amp},k_{sense},C_{imag},\frac{1}{\omega} $$
After that, the element transforms the isotropic equivalent coefficients ${k_{eq}, c_{eq}}$ into the rotor $x\text{–}y$ coordinates using the rotation defined by sensors_axis_rotation, producing $k_{xx}, k_{xy}, k_{yx}, k_{yy}$ and $c_{xx}, c_{xy}, c_{yx}, c_{yy}$.

These frequency-dependent arrays are stored in the element. During calls like run_modal (at a given speed), ROSS interpolates the coefficients at the required excitation frequencies to assemble the global stiffness and damping matrices for the rotor model.