Sensitivity computation when using Magnetic Bearings

[ArthurIasbeck]

Theoretical Background

The API STANDARD 617 document (produced by the American Petroleum Institute) requires that rotors supported by active magnetic bearings (AMBs) have their stability assessed based on a series of criteria. One of these criteria involves analyzing the sensitivity of the closed-loop control system that enables the AMB operation. Given the transfer functions of the controller $G_r(s)$ and the plant $G_p(s)$ (which in this case describes the relationship between the force applied by the AMB and the displacement of the shaft at the bearing location), it is possible to calculate the sensitivity $G_s(s)$ as follows:

$$ G_s(s) = \frac{1}{1 + G_r(s) G_p(s)} $$

Since the closed-loop response is given by

$$ G_c(s) = \frac{G_r(s) G_p(s)}{1 + G_r(s) G_p(s)} $$

it can be concluded that $G_s(s) + G_c(s) = 1$. Therefore,

$$ G_s(s) = 1 - G_c(s) $$

The block diagram below graphically represents the control system associated with the AMB. The force that the AMB applies to the shaft and the displacement of the shaft at the bearing location are represented by $F$ and $x$, respectively.

According to the ROSS documentation, the frequency response obtained from the execution of the run_freq_response method is computed based on the application of an excitation force at a given node of the shaft. This force, $F_w$, represents an excitation to the control system, as shown in the figure below.

Therefore, since ROSS produces a frequency response that relates $x$ and $F_d$, it can be said that this response represents, in the frequency domain, the transfer function $G_d(s)$, whose definition can be obtained based on the analysis of the block diagram introduced above.

$$ \frac{F_d(s)}{X(s)} = G_d(s) = \frac{G_p(s)}{1 + G_r(s) G_p(s)} $$

where $F_d(s)$ and $X(s)$ are the Laplace transforms of the disturbance and output, respectively.

Given the definition of $T(s)$ previously presented, it follows that

$$ T(s) = G_r(s) G_d(s) $$

Since the gains of the controller $G_r(s)$ (in this case, a PID controller) are known, it is possible to obtain the magnitude and phase of the controller's frequency response (using the Python Control Systems). Since the magnitude and phase of $G_d(s)$ are provided by ROSS (from the execution of the run_freq_response method), it is possible to compute the magnitude and phase of $T(s)$ as follows:

$$ |T(s)| = |G_r(s)||G_d(s)| $$

$$ \angle T(s) = \angle G_r(s) + \angle G_d(s) $$

Having obtained the magnitude and phase of $T(s)$, the sensitivity can finally be computed, given that $S(s) = 1 - T(s)$.

Method employed in the computation of sensitivity

The proposed method for computing the sensitivity is presented below. As previously discussed, the computation of the sensitivity depends on the gains of the controller employed and the frequency response provided by ROSS, which represents the relationship between the output and a disturbance applied at the plant input.

import control as ct
import cmath

def compute_sensitivity(kp_pid, ki_pid, kd_pid, mag_W, phase_W, speed_range):
    # Controller frequency response computation
    s = ct.tf("s")
    C = kp_pid + ki_pid / s + kd_pid * s
    mag_C, phase_C, _ = ct.frequency_response(C, speed_range)

    # Close-loop frequency response computation
    mag_T = mag_W * mag_C
    phase_T = phase_W + phase_C

    complex_T = [
        complex(
            mag_phase[0] * np.cos(mag_phase[1]), mag_phase[0] * np.sin(mag_phase[1])
        )
        for mag_phase in zip(mag_T, phase_T)
    ]

    # Sensitivity computation
    complex_S = [1 - z for z in complex_T]

    mag_S = [abs(z) for z in complex_S]
    phase_S = [cmath.phase(z) for z in complex_S]

    return mag_S, phase_S, speed_range

The script below was used to validate the aforementioned method. Both the rotor utilized, and the results obtained can be verified below, having been acquired from an execution that took 25.43 seconds. The computation of the sensitivity, itself, took only 1.09 milliseconds.

• Complete code

  import ross as rs
  import numpy as np
  import matplotlib.pyplot as plt
  import cmath
  from loguru import logger
  import control as ct
  import time
  
  max_freq = 100  # Hz
  x_label = "w [rad/s]"
  
  def get_dof(node, local_dof):
      local_dof_dict = {"x": 0, "y": 1, "alpha": 2, "beta": 3}
      return node * 4 + local_dof_dict[local_dof]
  
  def build_rotor(show_rotor=False):
      logger.info("Initiating rotor build.")
  
      # --- SHAFT MATERIAL DEFINITION ---
      steel = rs.Material(name="Steel", rho=7850, E=211e9, G_s=81.2e9)
  
      # --- SHAFT ELEMENTS DEFINITION ---
      L = 0.1  # Length of each element (m)
      i_d = 0  # Internal Diameter (m) - assuming solid shaft
      o_d = 0.05  # External Diameter (m)
  
      shaft_elements = [
          rs.ShaftElement(
              L=L,
              idl=i_d,
              odl=o_d,
              material=steel,
              shear_effects=True,
              rotary_inertia=True,
              gyroscopic=True,
          )
          for _ in range(10)
      ]
  
      # --- DISC DEFINITION ---
      disk = rs.DiskElement.from_geometry(
          n=5, material=steel, width=0.07, i_d=0.05, o_d=0.28
      )
  
      # --- MAGNETIC BEARING DEFINITION ---
      # Electromagnetic parameters (see ROSS documentation for details)
      g0 = 1e-3  # Air gap
      i0 = 1.0  # Bias current
      ag = 1e-4  # Pole area
      nw = 200  # Winding turns
      alpha = 0.392  # Half of the angle between two poles
  
      # PID gains
      kp_pid = 1000000  # Kp gain
      kd_pid = 1000000  # Kd gain
      k_amp = 1.0  # Power amplifier gain
      k_sense = 1.0  # Sensor gain
  
      # Magnetic bearing at node 2
      bearing1 = rs.MagneticBearingElement(
          n=2,
          g0=g0,
          i0=i0,
          ag=ag,
          nw=nw,
          alpha=alpha,
          kp_pid=kp_pid,
          ki_pid=0,
          kd_pid=kd_pid,
          k_amp=k_amp,
          k_sense=k_sense,
      )
  
      # Simple support at node 8
      kxx = 1e6  # Stiffness in x direction (N/m)
      kyy = 1e6  # Stiffness in y direction (N/m)
      cxx = 1e3  # Damping in x direction (N*s/m)
  
      bearing2 = rs.BearingElement(n=8, kxx=kxx, kyy=kyy, cxx=cxx)
  
      # --- ROTOR ASSEMBLY AND VISUALIZATION---
      rotor = rs.Rotor(
          shaft_elements=shaft_elements,
          disk_elements=[disk],
          bearing_elements=[bearing1, bearing2],
      )
      if show_rotor:
          rotor.plot_rotor().show()
  
      logger.info("Rotor build completed successfully.")
      return rotor
  
  def compute_freq_resp(rotor):
      logger.info("Initiating frequency response computation.")
  
      speed_range = np.linspace(0, max_freq * 2 * np.pi, 2 * max_freq)
      freq_resp = rotor.run_freq_response(speed_range=speed_range)
      np.save(
          "output/freq_resp",
          freq_resp.freq_resp,
      )
      np.save(
          "output/speed_range",
          np.array(freq_resp.speed_range),
      )
  
      logger.info("Frequency response computation completed successfully.")
  
  def get_freq_response_at_mma():
      logger.info("Collecting frequency response at AMB.")
  
      freq_response = np.load("output/freq_resp.npy")
      speed_range = np.load("output/speed_range.npy")
  
      left_bearing_node = 2
      disk_node = 5
  
      freq_response_at_mma = freq_response[
          get_dof(disk_node, "y"), get_dof(left_bearing_node, "y"), :
      ]
  
      mag_W = [abs(z) for z in freq_response_at_mma]
      phase_W = [cmath.phase(z) for z in freq_response_at_mma]
  
      logger.info("Frequency response at AMB collected successfully.")
      return mag_W, phase_W, speed_range
  
  def plot_freq_resp():
      logger.info("Initiating frequency response plot.")
  
      mag_W, phase_W, speed_range = get_freq_response_at_mma()
  
      # Displacement
      # Magnitude
      plt.figure(figsize=(8, 6), dpi=130)
      plt.subplot(2, 1, 1)
      plt.plot(speed_range, mag_W)
      plt.legend(['$\mathrm{G_w(s)}$'])
      plt.xlabel(x_label)
      plt.ylabel("Magnitude [m/N]")
      plt.title("Frequency Response")
      plt.yscale("log")
      plt.xlim([np.min(speed_range), np.max(speed_range)])
      plt.grid()
      plt.tight_layout()
  
      # Fase
      plt.subplot(2, 1, 2)
      plt.plot(speed_range, phase_W)
      plt.legend(['$\mathrm{G_w(s)}$'])
      plt.xlabel(x_label)
      plt.ylabel("Phase [rad]")
      plt.xlim([np.min(speed_range), np.max(speed_range)])
      plt.grid()
      plt.tight_layout()
  
      logger.info("Frequency response plot completed successfully.")
  
  def compute_sensitivity_call():
      logger.info("Initiating sensitivity computation.")
      mag_W, phase_W, speed_range = get_freq_response_at_mma()
      start_time = time.time()
      mag_S, phase_S, speed_range = compute_sensitivity(
          kp_pid=1000000,
          ki_pid=0,
          kd_pid=1000000,
          mag_W=mag_W,
          phase_W=phase_W,
          speed_range=speed_range,
      )
      end_time = time.time()
      logger.info(f'Sensitivity computation time: {end_time - start_time} seconds.')
      return mag_S, phase_S, speed_range
  
  def compute_sensitivity(kp_pid, ki_pid, kd_pid, mag_W, phase_W, speed_range):
      # Controller frequency response computation
      s = ct.tf("s")
      C = kp_pid + ki_pid / s + kd_pid * s
      mag_C, phase_C, _ = ct.frequency_response(C, speed_range)
  
      # Close-loop frequency response computation
      mag_T = mag_W * mag_C
      phase_T = phase_W + phase_C
  
      complex_T = [
          complex(
              mag_phase[0] * np.cos(mag_phase[1]), mag_phase[0] * np.sin(mag_phase[1])
          )
          for mag_phase in zip(mag_T, phase_T)
      ]
  
      # Sensitivity computation
      complex_S = [1 - z for z in complex_T]
  
      mag_S = [abs(z) for z in complex_S]
      phase_S = [cmath.phase(z) for z in complex_S]
  
      return mag_S, phase_S, speed_range
  
  def plot_sensitivity(mag_S, phase_S, speed_range):
      logger.info("Initiating sensitivity response plot.")
  
      # Displacement
      # Magnitude
      plt.figure(figsize=(8, 6), dpi=130)
      plt.subplot(2, 1, 1)
      plt.plot(speed_range, mag_S)
      plt.legend(['S(s)'])
      plt.xlabel(x_label)
      plt.ylabel("Magnitude [m/N]")
      plt.title("Frequency Response")
      plt.yscale("log")
      plt.xlim([np.min(speed_range), np.max(speed_range)])
      plt.grid()
      plt.tight_layout()
  
      # Fase
      plt.subplot(2, 1, 2)
      plt.plot(speed_range, phase_S)
      plt.legend(['S(s)'])
      plt.xlabel(x_label)
      plt.ylabel("Phase [rad]")
      plt.xlim([np.min(speed_range), np.max(speed_range)])
      plt.grid()
      plt.tight_layout()
  
      logger.info("sensitivity plot completed successfully.")
  
  def main():
      start_time = time.time()
      rotor = build_rotor(show_rotor=True)
      compute_freq_resp(rotor)
      plot_freq_resp()
      mag_S, phase_S, speed_range = compute_sensitivity_call()
      plot_sensitivity(mag_S, phase_S, speed_range)
      plt.show()
  
      logger.info("Execution completed successfully.")
      end_time = time.time()
  
      logger.info(f"Execution time: {end_time - start_time}.")
  
  if __name__ == "__main__":
      main()

Implementation of sensitivity computation in ROSS

The results presented thus far indicate that the proposed method does indeed allow the sensitivity to be correctly computed. Therefore, it is proposed that such computation should take place during the execution of the run_freq_response method. Once the frequency response (freq_resp) has been properly computed by this method, it will be necessary to verify if there is any AMB associated with the rotor under analysis. If so, the sensitivity should be calculated based on the frequency response produced during the execution of the run_freq_response method. The FrequencyResponseResult object will then need to be modified to receive information regarding sensitivity in its constructor. In addition, a new method should be added to this class, allowing the user to construct a plot that represents both the magnitude and phase of the sensitivity.

The diagrams below illustrate the current structure of the run_freq_response method and the FrequencyResponseResult class, while the blocks highlighted in pink indicate the modifications that will be implemented in both.

[ArthurIasbeck]

During the sensitivity computing implementation, it was necessary to slightly modify the initially proposed structure. The changes implemented in the ROSS fork (specifically at: https://github.com/ArthurIasbeck/ross/tree/compute-sensitivity) are described below.

1. Update of the run_freq_response method

   This method, responsible for computing the rotor's frequency response, will now be able to compute the sensitivity associated with any magnetic bearing present in the rotor. For this purpose, a new optional parameter has been added to the method: compute_sensitivite_at. This parameter is a dictionary containing the information required to compute the sensitivities of the magnetic bearings. Each key should be a bearing tag, and the corresponding value should be a dictionary defining the input and output degrees of freedom that determine the sensitivity computation, as shown in the example below. If the compute_sensitivite_at parameter receives a non-null value, the compute_sensitivity() method will be invoked.

   compute_sensitivite_at = {
       "Bearing 0": {"inp": 2, "out": 2},
       "Bearing 1": {"inp": 8, "out": 8},
   }

   When the sensitivity calculation is requested by the user, the FrequencyResponseResults object returned by the run_freq_response method will contain an attribute called sensitivity_results, which will hold an instance of the SensitivityResults class. In the example below, it is possible to check which attributes are associated with sensitivity_results.

   sensitivity_results = {SensitivityResults} <ross.results.SensitivityResults object at 0x7281ed7d4070>
   	compute_sensitivite_at = {dict: 2} {'Bearing 0': {'inp': 9, 'out': 9}, 'Bearing 1': {'inp': 33, 'out': 33}}
   	max_abs_sensitivities = {dict: 2} {'Bearing 0': 0.8329731089141272, 'Bearing 1': 0.8329731089184789}
   	sensitivities = {dict: 2} {'Bearing 0': [0.78451171-2.57357183e-17j 0.78470679+6.42813212e-05j, ...

2. Implementation of the SensitivityResults class

   This class is responsible for storing the results from the sensitivity calculation (not only the sensitivity values themselves but also the degrees of freedom that were used as the basis for the computation and the maximum absolute sensitivity values for each bearing). The class also has a method, named plot(), which returns a Plotly figure containing two subplots, in which the magnitude and phase of the sensitivity associated with each bearing are represented (as indicated below).

3. Update of the FrequencyResponseResults class

   The constructor of this class now receives the optional parameter sensitivity_results (an instance of the SensitivityResults class) which holds the results obtained from the magnetic bearing sensitivity computation. These results will only be generated if the rotor model includes at least one magnetic bearing and the user has requested sensitivity computation. The sensitivity_results parameter contains information about the computed sensitivities for each magnetic bearing. This includes a list of the associated magnetic bearings, the maximum absolute sensitivity for each bearing, and the individual sensitivity values for each bearing. It is worth noting that the FrequencyResponseResults instance, constructed by the run_freq_response() method will only contain the sensitivity_results field if the sensitivity has actually been calculated.

   The FrequencyResponseResults class also received the addition of the plot_sensitivity() method, which returns a Plotly figure containing two subplots, in which the magnitude and phase of the sensitivity associated with each bearing are represented. This method calls the plot() method of the SensitivityResults class, providing it with the frequency range over which the sensitivity should be plotted. The definition of the plot_sensitivity() method allows the Bode plots for sensitivity to be constructed directly from the object returned by the run_freq_response() method.

4. Definition of the compute_sensitivity() method

   This method computes the sensitivity for each magnetic bearing whose tag is present in the compute_sensitivity_at dictionary, given the frequency response of the rotor.

5. Creation of the amb_rotor_example method

   This function returns an instance of a simple rotor model. The model consists of ten shaft elements, one disk, one simple bearing, and one magnetic bearing. It provides a basic model for writing doctests. This method was used in defining the examples present in the description of the run_freq_response() method, which can now compute the sensitivity of the magnetic bearings associated with the rotor.

---

In summary, the user executes the run_freq_response method, passing the compute_sensitivity_at argument, which indicates the degrees of freedom where the sensitivity should be calculated. The run_freq_response then calls the compute_sensitivity method, which is responsible for actually computing the sensitivity values. Back in the run_freq_response method, an instance of the SensitivityResults class is created and used in the definition of an instance of the FrequencyResponseResults class, which is finally returned to the user.

[ArthurIasbeck]

Below are the latest changes implemented in the fork https://github.com/ArthurIasbeck/ross/tree/compute-sensitivity.

1. Update of the plot_sensitivity and plot methods of the FrequencyResponseResults and SensitivityResults classes, respectively, to allow the user to adjust the units of the sensitivity plots.
2. Implementation of the test_compute_sensitivity test method in the test_rotor_assembly module, where sensitivity computation validation is performed.
3. Inclusion of the fig_kwargs argument in the plot_sensitivity and plot methods to enable user configuration of sensitivity plots.

[ArthurIasbeck]

The changes that enable ROSS to perform the computation of sensitivities associated with magnetic bearings present in the rotor have been duly implemented in the fork https://github.com/ArthurIasbeck/ross/tree/compute-sensitivity and the pull request https://github.com/petrobras/ross/pull/1108 has been opened, requesting that these changes be incorporated into ROSS.

[ArthurIasbeck]

Recently, the following changes were made (commit ID 021e8e9 from fork https://github.com/ArthurIasbeck/ross/tree/compute-sensitivity):

• Creation of the SensitivityResults class
• Implementation of an exclusive method for calculating sensitivity (run_amb_sensitivity)
• Update of the rotor_amb_example method to include a more complex model
• Implementation of the test_save_load_sensitivityresponse test
• Update of the test_compute_sensitivity method

[ArthurIasbeck]

I recently had access to ISO 14839-3, which defines how the sensitivity of active magnetic bearings should be measured. This ISO standard assumes that the control system associated with the AMB can be represented as follows:

where $G_p$ is the rotor's transfer function, $E$ is a disturbance associated with the sensor signal (which can represent white noise or a harmonic signal), and $F_d$ is a force that represents a disturbance applied directly to the rotor.

The system above can be simplified so that:

Based on these representations, the ISO standard proposes that the open-loop, closed-loop, and sensitivity transfer functions, represented by $G_o$, $G_c$, and $G_s$, respectively, be computed as follows:

Also, according to the ISO standard, it is possible to define the open-loop transfer function as $G_o = G_c / (1 - G_c)$. Therefore, it is possible to verify that $G_s = 1 - G_c$.

Starting from the simplified diagram introduced previously, it is possible to verify that

$$ X = \frac{G_r G_p}{1 - G_r G_p}E + \frac{G_p}{1 - G_rG_p} F_d $$

Assuming that the frequency response provided by ROSS is obtained by varying $F_d$ and adopting $E = 0$, it follows that the transfer function associated with the frequency response in question is given by:

$$ G_d = \frac{X}{F_d} = \frac{G_p}{1 - G_r G_p} $$

Therefore, it follows that

$$ G_c = -G_d G_r = (-1 + 0j)G_dG_r $$

$$ |G_c| = |G_d||G_r| $$

$$ \angle{G_c} = \pi + \angle{G_d} + \angle{G_r} $$

Since the controller transfer function $G_r$ is known and the transfer function relating the disturbance to the displacement $G_d$ is provided (indirectly) by ROSS, it is possible to compute the closed-loop transfer function $G_c$ and use it in the calculation of the sensitivity $G_s$.

The recent changes in the methods responsible for calculating the sensitivity aimed to comply with ISO 14839-3. There were few significant changes, among them, the addition of $\pi$ in the calculation of the phase of $G_c$.

[ArthurIasbeck]

As mentioned in previous comments, ISO 14839-3 describes the experimental procedure associated with measuring the sensitivity of the Active Magnetic Bearing (AMB) control systems. This procedure is based on the sampling of some signals and the application of noise $E$ at the system's output, as indicated in the diagram below.

The time-domain data that enabled the computation of the sensitivity were obtained through the execution of the run_time_response method, as indicated below:

def generate_time_response():
    np.random.seed(0)

    rotor = self.rotor_assembly()
    speed = 1200  # [rad/s] = [Hz] * 2 * np.pi

    time_samples = int(self.t_max / self.dt + 1)
    t = np.linspace(0, self.t_max, time_samples)

    f = np.zeros((len(t), rotor.ndof))
    
    response = rotor.run_time_response(speed, f, t, method="newmark")

The method responsible for implementing the controllers associated with each AMBs (magnetic_bearing_controller) is summarized below. The changes that enabled the calculation of the sensitivity are highlighted. It is worth noting that these changes are temporary and only allow the information necessary for calculating the sensitivity to be written to the terminal so that it can be collected and subsequently processed. This approach was used only for testing purposes, and the data needed for calculating the sensitivity will be collected later in a different way. The magnetic_bearing_controller method was only implemented in the ROSS version found at https://github.com/mcarolalb/ross/tree/dev/magnetic-controller.

def magnetic_bearing_controller(self, magnetic_bearings, time_step, disp_resp):
    ...
    
    for elm in magnetic_bearings:
        dofs = [self.number_dof * elm.n, self.number_dof * elm.n + 1]
        for idx_dofs, value_dofs in enumerate(dofs):
            # CODE ADDED - START
            if value_dofs == 49:
                # Impulse
                 if elm.e0[idx_dofs] == 1e-06:
                     e = 1
                 else:
                     e = 0

                # White noise (uniform)
                # e = np.random.rand() * 1e-6

                # White noise (Gaussian)
                # e = np.random.normal(scale=1e-6)

                v_2 = disp_resp[value_dofs]
                v_1 = v_2 + e
                data_dict = {
                    "e": e,
                    "v_1": v_1,
                    "v_2": v_2
                }
                err = setpoint - v_1
            else:
                err = setpoint - disp_resp[value_dofs]
						# CODE ADDED - END
						
            P = elm.kp_pid * err
            elm.integral[idx_dofs] += elm.ki_pid * err * dt
            D = elm.kd_pid * (err - elm.e0[idx_dofs]) / dt
            ...

The data collected from the terminal after the execution of the run_time_response method were saved in the file input/00_temp_data_to_compute_sensitivity.txt and subsequently processed using the methods indicated below, which enabled the calculation of the sensitivity.

def compute_sensitivity(
    data_file_name="input/00_temp_data_to_compute_sensitivity.txt",
):
    with open(data_file_name, "r") as file:
        data = file.readlines()
        data = data[2:-1]
        e_list = []
        v_1_list = []
        for line in data:
            values = json.loads(line)
            e_list.append(values["e"])
            v_1_list.append(values["v_1"])

        e = np.array(e_list, dtype=np.float64)
        v_1 = np.array(v_1_list, dtype=np.float64)
        t = np.linspace(0, self.t_max, len(e))

        f, abs_G_s, phase_G_s, _ = compute_freq_resp_from_temp(
            t, v_1, e, plot_result=plot_result
        )
        omega = 2 * np.pi * f

        return omega, abs_G_s, phase_G_s
        

def compute_freq_resp_from_temp(t, y, u, plot_result=False, ft_gain=1):
    Y, _, _, _ = compute_fft(y, t)
    U, f, _, _ = compute_fft(u, t)

    Gyy = np.conj(Y) * Y
    Guu = np.conj(U) * U
    Guy = np.conj(U) * Y
    H = ft_gain * (Guy / Guu)
    coe = np.absolute(Guy) ** 2 / (Gyy * Guu)

    abs_H, phase_H = compute_abs_phase(H)

    return f, abs_H, phase_H, H
    

def compute_abs_phase(X):
    X = np.round(X, 10)
    abs_X = np.absolute(X)
    phase_X = np.array([0 if np.absolute(z) == 0 else np.angle(z) for z in X])
    return abs_X, phase_X

The rotor used to obtain the evaluated data is represented below. The sensitivity calculation, in this case, was based on the data associated with the $y$ displacement of node 12, where the first AMB is located.

In the graphs below, it is possible to verify the sensitivity values obtained considering different types of disturbance $E$. The types of disturbances used are described in the legend. The sensitivity values obtained based on the frequency response are shown in red. It is possible to verify that the best result was obtained when an impulsive disturbance was applied and that the results obtained based on the frequency response are definitely not reliable.