Sensitivity Computation

[ArthurIasbeck]

Pull Request Summary

This pull request is associated with the changes that enable ROSS to compute the sensitivities associated with magnetic bearings present in the rotor (as required by API 617).

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 $\mathrm{C(s)}$ and the plant $\mathrm{G(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 $\mathrm{S(s)}$ as follows:

$$ \mathrm{S(s)} = \frac{1}{1 + \mathrm{G(s)} \mathrm{C(s)}} $$

Since the closed-loop response is given by

$$ \mathrm{T(s)} = \frac{\mathrm{G(s) C(s)}}{1 + \mathrm{G(s)} \mathrm{C(s)}} $$

it can be concluded that $\mathrm{S(s)} + \mathrm{T(s)} = 1$. Therefore,

$$ \mathrm{S(s)} = 1 - \mathrm{T(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_w$, it can be said that this response represents, in the frequency domain, the transfer function $\mathbf{G_w(s)}$, whose definition can be obtained based on the analysis of the block diagram introduced above.

$$ \frac{\mathbf{W(s)}}{\mathbf{Y(s)}} = \mathbf{G_w(s)} = \frac{\mathbf{G(s)}}{1 + \mathbf{G(s)C(s)}} $$

where $\mathbf{W(s)}$ and $\mathbf{Y(s)}$ are the Laplace transforms of the disturbance and output, respectively.

Given the definition of $\mathbf{T(s)}$ previously presented, it follows that

$$ \mathbf{T(s)} = \mathbf{C(s)} \mathbf{G_w(s)} $$

Since the gains of the controller $\mathbf{C(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 package). Since the magnitude and phase of $\mathbf{G_w(s)}$ are provided by ROSS (from the execution of the run_freq_response method), it is possible to compute the magnitude and phase of $\mathbf{T(s)}$ as follows:

$$ |\mathbf{T(s)}| = |\mathbf{C(s)}| |\mathbf{G_w(s)}| $$

$$ \angle\mathbf{T(s)} = \angle\mathbf{C(s)} + \angle\mathbf{G_w(s)} $$

Having obtained the magnitude and phase of $\mathbf{T(s)}$, the sensitivity can finally be computed, given that $\mathrm{S(s)} = 1 - \mathrm{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.

def compute_sensitivity(speed_range, freq_resp, compute_sensitivity_at, ambs):
    """This method computes the sensitivity associated with the magnetic bearings in the rotor.

    Parameters
    ----------
    speed_range : array
        Array with the frequencies at which the sensitivity is computed.
    freq_resp : array
        Frequency response of the rotor (FrequencyResponseResults.freq_resp).
    compute_sensitivity_at : dict
        Dictionary defining the input and output degrees of freedom (DOFs) for
        sensitivity computation. The keys should be the tags of the magnetic bearings,
        and the values should be dictionaries with "inp" and "out" keys indicating
        the input and output DOFs.
    ambs : dict
        Dictionary of magnetic bearing elements associated with the rotor. The keys should be the tags of the
        magnetic bearings, and the values should be the magnetic bearing elements.

    Returns
    -------
    tuple
        A tuple containing two dictionaries:

        - max_abs_sensitivities: Dictionary containing the maximum absolute
          sensitivity value for each magnetic bearing. The keys are the magnetic
          bearing tags, and the values are the maximum sensitivity values.
        - sensitivities: Dictionary containing the sensitivity values for each
          magnetic bearing as a function of frequency. The keys are the magnetic
          bearing tags, and the values are complex arrays representing the
          sensitivity at each frequency point.

    Examples
    --------
    >>> rotor = rotor_amb_example()
    >>> compute_sensitivity_at = { "Bearing 0": {"inp": 27 * 4 + 1, "out": 12 * 4 + 1}, "Bearing 1": {"inp": 27 * 4 + 1, "out": 43 * 4 + 1} }
    >>> ambs = { bearing.tag: bearing for bearing in rotor.bearing_elements }
    >>> speed_range = np.linspace(0, 1000, 10)
    >>> frequency_response_result = rotor.run_freq_response(speed_range)
    >>> freq_resp = frequency_response_result.freq_resp
    >>> max_abs_sensitivities, sensitivities = compute_sensitivity(speed_range, freq_resp, compute_sensitivity_at, ambs)
    >>> max_abs_sensitivities # doctest: +ELLIPSIS
    {'Bearing 0': 1.0003602...
    >>> sensitivities # doctest: +ELLIPSIS
    {'Bearing 0': array([1.0002301...
    """
    sensitivities = {}
    max_abs_sensitivities = {}
    for amb_tag in compute_sensitivity_at.keys():
        amb = ambs[amb_tag]

        # Get disturbance frequency response
        freq_response_at_mma = freq_resp[
            compute_sensitivity_at[amb_tag]["inp"],
            compute_sensitivity_at[amb_tag]["out"],
            :,
        ]
        mag_w = [abs(z) for z in freq_response_at_mma]
        phase_w = [cmath.phase(z) for z in freq_response_at_mma]

        # Controller frequency response computation
        s = ct.tf("s")
        c = amb.kp_pid + amb.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]

        sensitivities[amb_tag] = np.array(complex_s)
        max_abs_sensitivities[amb_tag] = np.max(mag_s)

    return max_abs_sensitivities, sensitivities

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

The rotor used to validate the above method is depicted below and has two magnetic bearings, at nodes 12 and 27.

The script below enables the computation of sensitivity as well as its graphical representation.

from ross import rotor_amb_example
import numpy as np

rotor = rotor_amb_example()
speed_range = np.linspace(0, 1000 * 2 * np.pi, 2 * 1000)
compute_sensitivity_at = {"Bearing 0": {"inp": 27 * 4 + 1, "out": 12 * 4 + 1}}
sensitivity_result = rotor.run_amb_sensitivity(compute_sensitivity_at=compute_sensitivity_at, speed_range=speed_range)
sensitivity_result.plot(frequency_units="Hz", amplitude_units='m/N', phase_units='deg').show()

Modified files

• requirements.txt
• ross/results.py
• ross/rotor_assembly.py
• ross/tests/test_rotor_assembly.py
• ross/tests/test_results.py

Summary of changes made

1. Update of the run_amb_sensitivity method

   This method computes the sensitivity associated with the magnetic bearings in the rotor. Esse método possui como parâmetro obrigatório somente o compute_sensitivite_at, 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.

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

   Since the sensitivity is computed based on the frequency response, one of the other two optional parameters of the method must be provided: frequency_response_result or speed_range. The frequency_response_result parameter consists of a FrequencyResponseResults type object that represents the frequency response obtained by executing the run_freq_response method. The speed_range parameter, on the other hand, consists of an array containing the frequencies at which the sensitivities will be evaluated. If only the speed_range is provided, the run_freq_response method is invoked by the run_amb_sensitivity method and the frequency response is computed.

   This method returns a SensitivityResults type object, which contains the information listed below, regarding the sensitivity calculation.

   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. 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.

4. Creation of the rotor_amb_example method

   This function returns an instance of a rotor model. This rotor is supported by two magnetic bearings, and has a central disk, as shown in the figure below. It provides a basic model for writing doctests. This method was used in defining the examples present in the description of the run_amb_sensitivity method.

5. Implementation of the test_compute_sensitivity test method

   This method, implemented in the test_rotor_assembly.py module, performs a series of validations associated with sensitivity computation. The evaluated cases were as follows.

   1. Successful case where the sensitivity is computed for only one magnetic bearing
   2. Successful case where the sensitivity is computed for only one magnetic bearing but the frequency response is not provided
   3. Failing case where neither speed range nor frequency response are provided.
   4. Failing case where an invalid tag is provided in the sensitivity calculation.
   5. Failing case when there are no magnetic bearings on the rotor

6. Implementation of the test_save_load_sensitivityresponse

   This method, implemented in the test_rotor_assembly.py module, performs validations related to writing and reading TOML files that store the information present in a SensitivityResults object.

7. Addition of the Python Control Systems Library to requirements.txt

   The Python Control Systems Library was employed in the computation of the frequency response of the PID controllers associated with the magnetic bearings. Thus, the line below was added to the requirements.txt file.

   control
   

[codecov-commenter]

⚠️ Please install the to ensure uploads and comments are reliably processed by Codecov.

Codecov Report

Attention: Patch coverage is 98.54015% with 4 lines in your changes missing coverage. Please review.

Project coverage is 83.25%. Comparing base (f1bba84) to head (ab3a160).
Report is 73 commits behind head on main.

Files with missing lines	Patch %	Lines
ross/results.py	95.78%	4 Missing ⚠️

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1108      +/-   ##
==========================================
+ Coverage   81.67%   83.25%   +1.58%     
==========================================
  Files          38       38              
  Lines        8310     8272      -38     
==========================================
+ Hits         6787     6887     +100     
+ Misses       1523     1385     -138     

Files with missing lines	Coverage Δ
ross/bearing_seal_element.py	95.28% <100.00%> (+2.38%)	⬆️
ross/rotor_assembly.py	94.71% <100.00%> (+0.62%)	⬆️
ross/utils.py	95.51% <100.00%> (+0.55%)	⬆️
ross/results.py	80.16% <95.78%> (+1.55%)	⬆️

... and 3 files with indirect coverage changes

---

Continue to review full report in Codecov by Sentry.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 8df3e1c...ab3a160. Read the comment docs.

[ross-bott]

Hi there!
I have marked this issue as stale because it has not had activity for 45 days.
Consider the following options:

• If the issue refers to a large task, break it in smaller issues that can be solved in
  less than 45 days;
• Label the issue as wontfix or wontfix for now and close it.