Main Utilities

Be sure to check out Tutorials for example usages of these utilities!

Complementary Filter Synthesis

Optimal complementary filter synthesis using \(\mathcal{H}_\infty\) methods.

class kontrol.ComplementaryFilter(noise1=None, noise2=None, weight1=None, weight2=None, filter1=None, filter2=None, f=None)

Bases: object

Complementary filter synthesis class.

Parameters:
  • noise1 (TransferFunction) – Sensor noise 1 transfer function model. A transfer function that has magnitude response matching the amplitude spectral density of noise 1.
  • noise2 (TransferFunction) – Sensor noise 2 transfer function model. A transfer function that has magnitude response matching the amplitude spectral density of noise 2.
  • weight1 (TransferFunction, optional) – Weighting function 1. Frequency dependent specification for noise 1. Defaults None.
  • weight2 (TransferFunction, optional) – Weighting function 2. Frequency dependent specification for noise 2. Defaults None.
  • filter1 (TransferFunction, optional) – The complementary filter for noise 1. Defaults None.
  • filter2 (TransferFunction, optional) – The complementary filter for noise 2. Defaults None.
  • f (array, optional) – The frequency axis in Hz for evaluating the super sensor noise. Defaults None.
h2synthesis()

Complementary filter synthesis that minimizes the \(\mathcal{H}_2\)-norm of the super sensor noise. The noise must be specified before using this method.

hinfsynthesis()

Complementary filter synthesis that minimizes the \(\mathcal{H}_\infty\)-norm of the super sensor noise. The noise must be specified before using this method.

Notes

This is a utility class for complementary filter synthesis, sensor noise estimation and analysis. To use synthesis methods, specify noise1 and noise2 as the transfer function modesl for the sensor noises, and specify weight1 and weight2 as the frequency dependent specification for the two sensor noises.

References

[1]T. T. L. Tsang, T. G. F. Li, T. Dehaeze, C. Collette. Optimal Sensor Fusion Method for Active Vibration Isolation Systems in Ground-Based Gravitational-Wave Detectors. https://arxiv.org/pdf/2111.14355.pdf
f

The frequency axis in Hz.

filter1

First complementary filter.

filter2

Second complementary filter.

h2synthesis(clean_filter=True, **kwargs)

Synthesize complementary filters using H2 synthesis.

Returns:
  • filter1 (TransferFunction) – The complementary filter filtering noise 1.
  • filter2 (TransferFunction) – The complementary filter filtering noise 2.
  • clean_filter (boolean, optional) – Remove small outlier coefficients from filters. Defaults True.
  • **kwargs – Keyword arguments passed to kontrol.TransferFunction.clean()
hinfsynthesis(clean_filter=True, **kwargs)

Synthesize complementary filters using H-inifinity synthesis.

Returns:
  • filter1 (TransferFunction) – The complementary filter filtering noise 1.
  • filter2 (TransferFunction) – The complementary filter filtering noise 2.
  • clean_filter (boolean, optional) – Remove small outlier coefficients from filters. Defaults True.
  • **kwargs – Keyword arguments passed to kontrol.TransferFunction.clean()
noise1

Transfer function model of sensor noise 1.

noise2

Transfer function model of sensor noise 2.

noise_super(f=None, noise1=None, noise2=None, filter1=None, filter2=None)

Compute and return predicted the ASD of the super sensor noise

f : array, optional
The frequency axis in Hz. Use self.f if None. Defaults None.
noise1 : array, optional
The amplitude spectral density of noise 1. Use self.noise1 and self.f to estimate if None. Defaults None.
noise2 : array, optional
The amplitude spectral density of noise 2. Use self.noise1 and self.f to estimate if None. Defaults None.
filter1 : TransferFunction, optional
The complementary filter for filtering noise1 Use self.filter1 if not specified. Defaults None
filter2 : TransferFunction, optional
The complementary filter for filtering noise1 Use self.filter2 if not specified. Defaults None
Returns:The amplitude spectral density of the super sensor noise.
Return type:array

Frequency Series Class

Frequency domain and transfer function data modeling

class kontrol.FrequencySeries(f, x)

Bases: object

Frequency series class

f

Frequency axis

Type:array
x

Frequency series data.

Type:array
x_empirical

Frequency series of the empirical model.

Type:array or None
model_empirical

The model function whose first argument is the frequency axis and the rest being an arbitrary number of model parameters to be fit.

Type:func(f, a, b, c, ..) -> array or None
args_empirical_model

The optimized arguments passed to the model_empirical.

Type:array or None
f_zpk

The frequency axis used during ZPK fitting.

Type:array or None
x_zpk

The frequency series data used during ZPK fitting. Note: this is a complex array.

Type:array or None
tf_zpk

The transfer function the best fit the data using ZPK model.

Type:control.xferfcn.TransferFunction or None
args_zpk_model

A 1-D list of zeros, poles, and gain. Zeros and poles are in unit of Hz.

Type:array or None
f_tf

The frequency series data used during transfer function fitting.

Type:array or None
x_tf

The frequency series data used during transfer function fitting. Note: this is a complex array.

Type:array or None
tf

The modeled transfer function.

Type:control.xferfcn.TransferFunction or None
args_tf_model

A 1-D list of numerator and denominator coefficients, from higher order to lower order.

Type:array or None
fit_empirical(model, x0=None, error_func=<function log_mse>, error_func_kwargs={}, optimizer=<function minimize>, optimizer_kwargs={})

Fit frequency series data with an empirical model. (Local optimize)

Parameters:
  • model (func(f, a, b, c, ..) -> array) – The model function whose first argument is the frequency axis and the rest being an arbitrary number of model parameters to be fit.
  • x0 (array, optional) – The initial parameters for fitting. Defaults None. If None, defaults to np.ones()
  • error_func (func(x1: array, x2: array) -> float) – The function that evaluate the error between arrays x1 and x2. Defaults to kontrol.core.math.log_mse, which evaluates the logarithmic mean square error.
  • error_func_kwargs (dict, optional) – Keyword arguments passed to the error function. Use this to specify weighting function {“weight”: weight_func}. Defaults to {}.
  • optimizer (func, optional) – The optimization function which has the signature func(func, args, bounds, **kw) -> scipy.optimize.OptimizeResult. Defaults to scipy.optimize.minimize.
  • optimizer_kwargs (dict, optional) – Keyword arguments passed to the optimization function. Defaults {} but will set to {“method”: “Powell”, “x0”=x0} if optimizer is default.
Returns:

res – The result of the fitting (optimization)
Return type:

scipy.optimize.OptimizeResult
fit_tf(x0=None, log_var=True, error_func=<function log_mse>, error_func_kwargs={}, optimizer=<function minimize>, optimizer_kwargs={})

Fit frequency series with a transfer function model

Parameters:
  • x0 (array or None, optional) – The initial parameters defining the initial transfer function. A 1-D list of numerator and denominator coefficients, from higher order to lower order. Defaults None. If None, will use the result from the ZPK fitting.
  • log_var (boolean, optional) – Optimize the log of variables instead. Useful for variables that have very large dynamic range Defaults True.
  • error_func (func(x1: array, x2: array) -> float, optional) – The function that evaluate the error between arrays x1 and x2. Defaults to kontrol.core.math.log_mse, which evaluates the logarithmic mean square error.
  • error_func_kwargs (dict, optional) – Keyword arguments passed to the error function. Use this to specify weighting function {“weight”: weight_func}. Defaults to {}.
  • optimizer (func, optional) – The optimization function which has the signature func(func, args, bounds, **kw) -> scipy.optimize.OptimizeResult. Defaults to scipy.optimize.minimize.
  • optimizer_kwargs (dict, optional) – Keyword arguments passed to the optimization function. Defaults {} but will set to {“method”: “Powell”, “x0”: x0} if optimizer is default.
Returns:

res – The result of the optimization.
Return type:

scipy.optimize.OptimizerResult

Note

Only use this after using fit_zpk().

fit_zpk(order, fit='x_empirical', padding=False, padding_order=1.0, error_func=<function log_mse>, error_func_kwargs={}, optimizer=<function differential_evolution>, optimizer_kwargs={})

Fit frequency series with a zpk model.

Parameters:
  • order (int) – The order of the ZPK model to be used.
  • fit (str, optional) – Choose the data to be fit, choose from [“x”, “x_empirical”].
  • padding (boolean, optional) – Pad the data by continuous white noise. Defaults False.
  • padding_order (float, optional) – The number of decades to be padded on both sides of the series. Defaults to 1.
  • error_func (func(x1: array, x2: array) -> float) – The function that evaluate the error between arrays x1 and x2. Defaults to kontrol.core.math.log_mse, which evaluates the logarithmic mean square error.
  • error_func_kwargs (dict, optional) – Keyword arguments passed to the error function. Use this to specify weighting function {“weight”: weight_func}. Defaults to {}.
  • optimizer (func, optional) – The optimization function which has the signature func(func, args, bounds, **kw) -> scipy.optimize.OptimizeResult. Defaults to scipy.optimize.differential_evolution.
  • optimizer_kwargs (dict, optional) – Keyword arguments passed to the optimization function. Defaults {} but will set to {“workers”: -1, “updating”:”deferred”} if optimizer is default.
Returns:

res – The result of the optimization.
Return type:

scipy.optimize.OptimizerResult

Note

For now, we support models with equal number of zeros and poles only.

Best use with log-space frequency series or after using fit_empirical().

Sensors and Actuators Utilities

Sensing Matrix Classes

Sensing matrix and diagonalization.

class kontrol.SensingMatrix(matrix, coupling_matrix=None, *args, **kwargs)

Bases: kontrol.sensact.matrix.Matrix

Base class for sensing matrices

diagonalize(coupling_matrix=None)

Diagonalize the sensing matrix, given a coupling matrix.

Parameters:coupling_matrix (array, optional.) – The coupling matrix. If None, self.coupling_matrix will be used. Default None.
Returns:The new sensing matrix.
Return type:array

Notes

The coupling matrix has (i, j) elements as coupling ratios x_i/x_j. For example, consider the 2-sensor configuration: I have a coupled sensing readout \(x_{1,\mathrm{coupled}}\) that reads \(x_{1,\mathrm{coupled}}=x_1 + 0.1x_2\). and, I have another coupled sensing readout \(x_{2,\mathrm{coupled}}\) that reads \(x_{2,\mathrm{coupled}}=-0.2x_1 + x_2\). Then, the coupling matrix is

\[\begin{split}\begin{bmatrix} 1 & 0.1\\ -0.2 & 1 \end{bmatrix}.\end{split}\]

Optical Lever Sensing Matrices

Sensing matrices for optical levers in KAGRA.

class kontrol.OpticalLeverSensingMatrix(*args, **kwargs)

Bases: kontrol.sensact.matrix.SensingMatrix

Optical lever sensing matrix base class.

This is a sensing matrix that maps tilt-sensing QPD and length-sensing QPD (placed behind a lens) readouts to the suspended optics’ longitudinal, pitch, and yaw displacements.

Parameters:
  • r_h (float) – Lever arm from the optics to the tilt-sensing QPD plane on the horizontal plane (amplifying yaw).
  • r_v (float) – Lever arm from the optics to the tilt-sensing QPD plane on the vertical plane (amplifying pitch).
  • alpha_h (float) – Angle of incidence on the horizontal plane.
  • alpha_v (float) – Angle of incidence on the vertical plane.
  • r_lens_h (float, optional) – Lever arm from optics to the lens on the horizontal plane. Defaults None. Specify if it exists.
  • r_lens_v (float, optional) – Lever arm from optics to the lens on the vertical plane. Defaults None. Specify if it exists.
  • d_h (float, optional) – Horizontal distance from the lens to the length-sensing QPD. Defaults None.
  • d_v (float, optional) – Vertical distance from the lens to the length-sensing QPD. Defaults None.
  • delta_x (float, optional) – Horizontal miscentering of the beam spot at the optics plane.
  • delta_y (float, optional) – Vertical miscentering of the beam spot at the optics plane.
  • phi_tilt (float, optional) – Angle from the tilt-sensing QPD frame to the yaw-pitch frame. Defaults None.
  • phi_len (float, optional) – Angle from the length-sensing QPD frame to the yaw-pitch frame. Defaults None.
  • f_h (float, optional,) – Focal length of the lens projected to the horizontal plane. Defaults np.inf (no lens).
  • f_v (float, optional,) – Focal length of the lens projected to the vertical plane. Defaults np.inf (no lens).
  • format (str, optional) –

    Format of the sensing matrix. Choose from

    ”OPLEV2EUL”: Default sensing matrix from KAGRA MEDM screen
    with input (TILT_PIT, TILT_YAW, LEN_PIT, LEN_YAW), and output (longitudinal, pitch and yaw).

    ”xy”: Matrix as shown in [1]_.

  • coupling_matrix (array, optional) – The coupling matrix. Default None.

Notes

We’re using equation (29) from [1]_.

\[\begin{split}\begin{pmatrix} x_L\\ \theta_P\\ \theta_Y \end{pmatrix} = \mathbf{C}_\mathrm{miscenter} \mathbf{C}_\mathrm{align} \mathbf{C}_\mathrm{rotation} \begin{pmatrix} x_\mathrm{tilt}\\ y_\mathrm{tilt}\\ x_\mathrm{len}\\ y_\mathrm{len} \end{pmatrix},\end{split}\]

where \(x_L\) is the longitudinal displacement of the optics, \(\theta_P\) is the pitch angular displacement of the optics, \(\theta_Y\) is the yaw angular displacement of the optics, \(x_\mathrm{tilt}\) is the horizontal displacement of the beam spot at the tilt-sensing QPD plane, \(y_\mathrm{tilt}\) is the vertical displacement of the beam spot at the tilt-sensing QPD plane, \(x_\mathrm{len}\) is the horizontal displacement of the beam spot at the length-sensing QPD plane, \(y_\mathrm{len}\) is the vertical displacement of the beam spot at the length-sensing QPD plane,

\[\begin{split}\mathbf{C}_\mathrm{rotation} = \begin{bmatrix} \cos\phi_\mathrm{tilt} & \sin\phi_\mathrm{tilt} & 0 & 0\\ -\sin\phi_\mathrm{tilt} & \cos\phi_\mathrm{tilt} & 0 & 0\\ 0 & 0 & \cos\phi_\mathrm{len} & \sin\phi_\mathrm{len}\\ 0 & 0 & -\sin\phi_\mathrm{len} & \cos\phi_\mathrm{len} \end{bmatrix},\end{split}\]
\[\begin{split}\mathbf{C}_\mathrm{align} = \begin{bmatrix} 2\sin\alpha_h & 0 & 2r_h\\ 2\sin\alpha_v & 2r_v & 0\\ 2\sin\alpha_h\left(1-\frac{d_h}{f_h}\right) & 0 & 2\left[\left(1-\frac{d_h}{f_h}\right)r_{\mathrm{lens},h} + d_h\right]\\ 2\sin\alpha_v\left(1-\frac{d_v}{f_v}\right) & 2\left[\left(1-\frac{d_v}{f_v}\right)r_{\mathrm{lens},v} + d_v\right] & 0\\ \end{bmatrix}^{+},\end{split}\]
\[\begin{split}\mathbf{C}_\mathrm{miscenter} = \begin{bmatrix} 1 & \delta_y & \delta_x\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{bmatrix}^{-1},\end{split}\]

\(\phi_\mathrm{tilt}\) is the angle between the tilt-sensing QPD and the yaw-pitch frame, \(\phi_\mathrm{len}\) is the angle between the length-sensing QPD and the yaw-pitch frame, \(r_h\) is the lever arm on the horzontal plane, \(r_v\) is the lever arm on the vertical plane, \(\alpha_h\) is the angle of incidence on the horizontal plane, \(\alpha_v\) is the angle of incidence on the vertical plane, \(r_{\mathrm{lens}, h}\) is the lever arm between the optics and the lens on the horizontal plane, \(r_{\mathrm{lens}, v}\) is the lever arm between the optics and the lens on the vertical plane, \(d_h\) is the distance between the lens and the length-sensing QPD on the horizontal plane, \(d_v\) is the distance between lens and the length-sensing QPD on the vertical plane, and \(f\) is the focal length of the convex lens.

References

[1]Tsang Terrence Tak Lun, Sensing Matrices for Optical Levers of the KAGRA Main Optics, https://github.com/terrencetec/kagra-optical-lever.
alpha_h

Angle of incidence on the horizontal plane.

alpha_v

Angle of incidence on the vertical plane.

d_h

Horizontal distance from the lens to the length-sensing QPD.

d_v

Vertical distance from the lens to the length-sensing QPD.

delta_x

Horizontal miscentering of the beam spot at the optics plane.

delta_y

Horizontal miscentering of the beam spot at the optics plane.

f_h

Focal length of the convex lens projected on the horizontal plane.

f_v

Focal length of the convex lens projected on the vertical plane.

format

Format of the sensing matrix.

Choose from
“OPLEV2EUL”: Default sensing matrix from KAGRA MEDM screen
with input (TILT_PIT, TILT_YAW, LEN_PIT, LEN_YAW), and output (longitudinal, pitch and yaw).

“xy”: Matrix as shown in [1]_.

References

[1]Tsang Terrence Tak Lun, Sensing Matrices for Optical Levers of the KAGRA Main Optics, https://github.com/terrencetec/kagra-optical-lever.
phi_len

Angle from the length-sensing QPD frame to the yaw-pitch frame.

phi_tilt

Angle from the tilt-sensing QPD frame to the yaw-pitch frame.

r_h

Lever arm from optics to tilt-sensing QPD on the horizontal plane.

r_lens_h

Lever arm from optics to the lens on the horizontal plane.

r_lens_v

Lever arm from optics to the lens on the vertical plane.

r_v

Lever arm from optics to tilt-sensing QPD on the vertical plane.

update_matrices_decorator()

Update matrices and self upon setting new parameters.

Horizontal Optical Lever Sensing Matrices

Sensing matrices for horizontal optical levers (Type-A, Type-Bp) in KAGRA.

class kontrol.HorizontalOpticalLeverSensingMatrix(*args, **kwargs)

Bases: kontrol.sensact.optical_lever.OpticalLeverSensingMatrix

Horizontal optical lever sensing matrix.

Parameters:
  • r (float) – Lever arm.
  • alpha_h (float) – Angle of incidence on the horizontal plane.
  • r_lens (float, optional) – Lever arm from the optics to the convex lens. Default 0.
  • f (float, optional) – Focal length of the convex lens. Default np.inf.
  • alpha_v (float, optional) – Angle of incidence on the vertical plane. Default 0.
  • phi_tilt (float, optional) – Angle from the tilt-sensing QPD frame to the yaw-pitch frame. Default 0.
  • phi_len (float, optional) – Angle from the length-sensing QPD frame to the yaw-pitch frame. Default 0.
  • delta_x (float, optional) – Horizontal miscentering of the beam spot at the optics plane. Default 0.
  • delta_y (float, optional) – Vertical miscentering of the beam spot at the optics plane. Default 0.
  • delta_d (float, optional) – Misplacement of the length-sensing QPD. Default 0.
  • format (str, optional) –

    Format of the sensing matrix. Choose from

    ”OPLEV2EUL”: Default sensing matrix from KAGRA MEDM screen
    with input (TILT_PIT, TILT_YAW, LEN_PIT, LEN_YAW), and output (longitudinal, pitch and yaw).

    ”xy”: Matrix as shown in [1]_.

  • coupling_matrix (array, optional) – The coupling matrix. Default None.
  • *args – Variable length arguments passed to OpticalLeverSensingMatrix.
  • **kwargs – Keyword arguments passed to OpticalLeverSensingMatrix.
delta_d

Misplacement of the length-sensing QPD.

f

Focal length of the convex lens.

r

Lever arm.

r_lens

Lever arm from the optics to the convex lens.

Vertical Optical Lever Sensing Matrices

Sensing matrices for vertical optical levers (Type-B) in KAGRA.

class kontrol.VerticalOpticalLeverSensingMatrix(*args, **kwargs)

Bases: kontrol.sensact.optical_lever.OpticalLeverSensingMatrix

Vertical optical lever sensing matrix.

Parameters:
  • r (float) – Lever arm.
  • alpha_v (float) – Angle of incidence on the vertical plane.
  • r_lens (float, optional) – Lever arm from the optics to the convex lens. Default 0.
  • f (float, optional) – Focal length of the convex lens. Default np.inf.
  • alpha_h (float, optional) – Angle of incidence on the horizontal plane. Default 0.
  • phi_tilt (float, optional) – Angle from the tilt-sensing QPD frame to the yaw-pitch frame. Default 0.
  • phi_len (float, optional) – Angle from the length-sensing QPD frame to the yaw-pitch frame. Default 0.
  • delta_x (float, optional) – Horizontal miscentering of the beam spot at the optics plane. Default 0.
  • delta_y (float, optional) – Vertical miscentering of the beam spot at the optics plane. Default 0.
  • delta_d (float, optional) – Misplacement of the length-sensing QPD. Default 0.
  • format (str, optional) –

    Format of the sensing matrix. Choose from

    ”OPLEV2EUL”: Default sensing matrix from KAGRA MEDM screen
    with input (TILT_PIT, TILT_YAW, LEN_PIT, LEN_YAW), and output (longitudinal, pitch and yaw).

    ”xy”: Matrix as shown in [1]_.

  • coupling_matrix (array, optional) – The coupling matrix. Default None.
  • *args – Variable length arguments passed to OpticalLeverSensingMatrix.
  • **kwargs – Keyword arguments passed to OpticalLeverSensingMatrix.
delta_d

Misplacement of the length-sensing QPD.

f

Focal length of the convex lens.

r

Lever arm.

r_lens

Lever arm from the optics to the convex lens.

Sensor Calibration Functions

Functions for calibrating sensor readouts.

Calibration library for calibrating sensors from sensors measurements

kontrol.sensact.calibration.calibrate(xdata, ydata, method='linear', **kwargs)

Fit the measurement data and returns the slope of the fit.

Parameters:
  • xdata (array) – The independent variable, e.g. the displacement of the sensor.
  • ydata (array) – The dependent variable, e.g. the sensor readout.
  • method (str, optional) –
    The method of the fit.
    Choose from [“linear”, “erf”].
    • ”linear”: use kontrol.sensact.calibration.calibrate_linear()
    • ”erf”: use kontrol.sensact.calibration.calibrate_erf()
  • **kwargs – Keyword arguments passed to the calibration methods.
Returns:

  • slope (float) – The slope of the fitted straight line.
  • intercept (float) – The y-intercept of the fitted straight line.
  • linear_range (float, optional) – The range of y where the sensor considered linear.
kontrol.sensact.calibration.calibrate_erf(xdata, ydata, nonlinearity=5, return_linear_range=False, return_model=False)

Fit an error function and return the 1st-order slope and intercept.

Parameters:
  • xdata (array) – The independent variable, e.g. the displacement of the sensor.
  • ydata (array) – The dependent variable, e.g. the sensor readout.
  • nonlinearity (float, optional) – The specification of non-linearity (%). Defined by the maximum deviation of the full_range. Default 5.
  • return_linear_range (boolean, optional) – If True, return the linear range of ydata.
  • return_model (boolean, optional) – Return a kontrol.curvefit.model.Model object with the fitted parameters.
Returns:

  • slope (float) – The slope of the fitted straight line.
  • intercept (float) – The \(y\)-intercept of the fitted straight line.
  • linear_range (float, optional) – The range of x where the sensor considered linear.

Notes

Credits to Kouseki Miyo, the inventor of this method.

We fit the following function

\[f(x; a, m, x_0, y_0) = a\,\mathrm{erf}(m(x-x_0)) + y_0\,\,\]

where \(a, m, x_0, y_0\) are some parameters to be found. The \(\mathrm{erf}(x)\) function is defined as

\[\mathrm{erf}(x) = \frac{2}{\sqrt{\pi}}\int_0^x\,e^{-x^2}\,dx\,.\]

If we taylor expand the exponential function and take the first-order approximation of the \(\mathrm{erf}(x)\) function, we get

\[\mathrm{erf}(x) \approx \frac{2am}{\sqrt{\pi}}(x-x_0) + y_0\,.\]

So the (inverse) calibration factor is \(\frac{2am}{\sqrt{\pi}}\), and the \(y\)-intercept is \(\frac{2am}{\sqrt{\pi}}x_0 + y_0\).

kontrol.sensact.calibration.calibrate_linear(xdata, ydata, nonlinearity=5, full_range=None, start_index=None, return_linear_range=False, return_model=False)

Fit a straight line to the data and returns the slope and intercept

This functions recursively fits a straight line to chosen data points and terminates when no linear data point remains. It starts from 3 points closest to the middle of the full-range. After fitting a straight line to the 3 points, other data points that are within the linearity specification will be included to the dataset and this dataset will be use for the next fit. The process repeats until no data points can be added anymore.

Parameters:
  • xdata (array) – The independent variable, e.g. the displacement of the sensor.
  • ydata (array) – The dependent variable, e.g. the sensor readout.
  • nonlinearity (float, optional) – The specification of non-linearity (%). Defined by the maximum deviation of the full_range. Default 5.
  • full_range (float, optional) – The full output range of the sensor. If not specified, it will be chosen to be max(ydata) - min(ydata)
  • start_index (int, optional) – The index of the data point to start with. If not specified, it will be chosen to be index of the ydata closest to (max(ydata) + min(ydata))/2.
  • return_linear_range (boolean, optional) – If True, return the linear range of ydata.
  • return_model (boolean, optional) – Return a kontrol.curvefit.model.Model object with the fitted parameters.
Returns:

  • slope (float) – The slope of the fitted straight line.
  • intercept (float) – The y-intercept of the fitted straight line.
  • linear_range (float, optional) – The range of x where the sensor considered linear.
  • model (kontrol.curvefit.model.Model) – The fitted model.

Spectral Analysis Functions

Spectral analysis related functions library.

Spectral analysis related function library.

kontrol.core.spectral.asd2rms(asd, f=None, df=1.0, return_series=True)

Calculate root-mean-squared value from amplitude spectral density

Parameters:
  • asd (array) – The amplitude spectral density
  • f (array, optional) – The frequency axis. Defaults None.
  • df (float, optional) – The frequency spacing. Defaults 1.
  • return_series (bool, optional) – Returns the RMS as a frequency series, where each value is the integrated RMS from highest frequency. If False, returns a single RMS value for the whole spectrum. Defaults True.
Returns:

  • array – The integrated RMS series, if return_series==True.
  • float – The integrated RMS value, if return_series==False.

Notes

The integrated RMS series is defined as

\[x_\mathrm{RMS}(f) = \int_\infty^{f}\,x(f')^2\,df'\,,\]

where \(x(f)\) is the amplitude spectral density.

When return_series is False, only \(x_\mathrm{RMS}(0)\) is returned.

kontrol.core.spectral.asd2ts(asd, f=None, fs=None, t=None, window=None, zero_mean=True)

Simulate time series from amplitude spectral density

Parameters:
  • asd (array) – The amplitude spectral density.
  • f (array, optional) – The frequency axis of the amplitude spectral density. This is used to calculate the sampling frequency. If fs is not specified, f must be specified. Defaults None.
  • fs (float, optional) – The sampling frequency in Hz. If f is specified, the last element of f will be treated as the Nyquist frequency so the double of it would be the sampling frequency. Defaults None.
  • t (array, optional) – The desired time axis for the time series. If t is specified, then the time series will be interpolated using numpy.interp() assuming that the time series is periodic. Default None.
  • window (array, optional) – The FFT window used to compute the amplitude spectral density. Defaults None.
  • zero_mean (boolean, optional) – Returns a time series with zero mean. Defaults True.
Returns:

  • t (array) – Time axis.
  • time_series (array) – The time series

Notes

Using a custom time axis is not recommended as interpolation leads to distortion of the signal. To extend the signal in time, it’s recommended to repeat the original time series instead.

kontrol.core.spectral.pad_above_maxima(series, pad_index=-1, **kwargs)

Pad series above local maxima

Parameters:
  • series (array) – The series to be padded.
  • pad_index (int, optional) – The index of the local maxima. Defaults to -1.
  • **kwargs – Keyword arguments passed to scipy.signal.argrelmin()
Returns:

series_padded – The padded series
Return type:

array
kontrol.core.spectral.pad_above_minima(series, pad_index=-1, **kwargs)

Pad series above local minima

Parameters:
  • series (array) – The series to be padded.
  • pad_index (int, optional) – The index of the local minima. Defaults to -1.
  • **kwargs – Keyword arguments passed to scipy.signal.argrelmin()
Returns:

series_padded – The padded series
Return type:

array
kontrol.core.spectral.pad_below_maxima(series, pad_index=0, **kwargs)

Pad series below local maxima

Parameters:
  • series (array) – The series to be padded.
  • pad_index (int, optional) – The index of the local maxima. Defaults to 0.
  • **kwargs – Keyword arguments passed to scipy.signal.argrelmin()
Returns:

series_padded – The padded series
Return type:

array
kontrol.core.spectral.pad_below_minima(series, pad_index=0, **kwargs)

Pad series below local minima

Parameters:
  • series (array) – The series to be padded.
  • pad_index (int, optional) – The index of the local minima. Defaults to 0.
  • **kwargs – Keyword arguments passed to scipy.signal.argrelmin()
Returns:

series_padded – The padded series
Return type:

array
kontrol.core.spectral.three_channel_correlation(psd1, psd2=None, psd3=None, csd12=None, csd13=None, csd21=None, csd23=None, csd31=None, csd32=None, returnall=True)

Noise estimation from three sensors’ readouts.

Parameters:
  • psd1 (array) – The power spectral density of the readout of the first sensor.
  • psd2 (array, optional) – The power spectral density of the readout of the second sensor. Defaults None.
  • psd3 (array, optional) – The power spectral density of the readout of the third sensor. Defaults None.
  • csd12 (array, optional) – Cross power spectral density between readout 1 and 2. If not specified, this will be estimated as psd1*psd2/csd21. Default None.
  • csd13 (array, optional) – Cross power spectral density between readout 1 and 3. If not specified, this will be estimated as psd1*psd3/csd31. Default None.
  • csd21 (array, optional) – Cross power spectral density between readout 2 and 1. If not specified, this will be estimated as psd2*psd1/csd12. Default None.
  • csd23 (array, optional) – Cross power spectral density between readout 2 and 3. If not specified, this will be estimated as psd2*psd3/csd32. Default None.
  • csd31 (array, optional) – Cross power spectral density between readout 3 and 1. If not specified, this will be estimated as psd1*psd3/csd13. Default None.
  • csd32 (array, optional) – Cross power spectral density between readout 3 and 1. If not specified, this will be estimated as psd3*psd2/csd23. Default None.
  • returnall (boolean, optional) – If True, return all three noise estimations. If False, return noise estimation of first sensor only. Defaults True.
Returns:

  • noise1 (array) – Power spectral density of the estimated noise in psd1.
  • noise2 (array, optional) – Power spectral density of the estimated noise in psd2. Returns only if returnall==True
  • noise3 (array, optional) – Power spectral density of the estimated noise in psd3. Returns only if returnall==True

Notes

The PSD of the noise in psd1 is then computed as

\[P_{n_1n_1}(f) = \left\lvert P_{x_1x_1}(f) - \frac{P_{x_1x_3}(f)}{P_{x_2x_3}(f)}P_{x_2x_1}\right\vert\,,\]

If returnall is True, at least psd1, psd2, psd3, (csd13 or csd31), (csd23 or csd32), and (csd12 and csd21) must be provided.

References

[2]R. Sleeman, A. Wettum, and J. Trampert. Three-channel correlation analysis: A new technique to measure instrumental noise of digitizers and seismic sensors. Bulletin of the Seismological Society of America, 96:258–271, 2006.
kontrol.core.spectral.two_channel_correlation(psd, coh)

Noise estimation from two identical sensors’ readouts.

Parameters:
  • psd (array) – The power spectral density of the readout of the sensor
  • coh (array) – The coherence between readout of the sensor and another identical sensor.
Returns:

noise – Power spectral density of the estimated noise.
Return type:

array

Notes

The PSD of the noise is computed as

\[P_{nn}(f) = P_{x_1x_1}(f)\left(1-C_{x_1x_2}(f)^{\frac{1}{2}}\right)\,,\]

where \(P_{x_1x_1}(f)\) is the power spectral density of the readout and \(C_{x_1x_2}(f)\) is the coherence between the two sensor readouts.

References

[1]Aaron Barzilai, Tom VanZandt, and Tom Kenny. Technique for measurement of the noise of a sensor in the presence of large background signals. Review of Scientific Instruments, 69:2767–2772, 07 1998.

Foton Utilities

KAGRA/LIGO Foton related utilities.

KAGRA/LIGO Foton related utilities.

kontrol.core.foton.foton2tf(foton_string)

Convert a Foton string to TransferFunction

Parameters:foton_string (str) – The Foton string, e.g. zpk([0], [1; 1], 1, “n”).
Returns:tf – The transfer function.
Return type:TransferFunction
kontrol.core.foton.notch(frequency, q, depth, significant_figures=6)

Returns the foton expression of a notch filter.

Parameters:
  • frequency (float) – The notch frequency (Hz).
  • q (float) – The quality factor.
  • depth (float) – The depth of the notch filter (magnitude).
  • significant_figures (int, optional) – Number of significant figures to print out. Defaults to 6.
Returns:

The foton representation of this notch filter.
Return type:

str
kontrol.core.foton.rpoly2tf(foton_string)

Converts rpoly Foton strings to TransferFunction

Parameters:foton_string (str) – The rpoly Foton string. E.g. rpoly([1; 2; 3], [2; 3; 4], 5)
kontrol.core.foton.tf2foton(tf, expression='zpk', root_location='s', significant_figures=6, itol=1e-25, epsilon=1e-25)

Convert a single transfer function to foton expression.

Parameters:
  • tf (TransferFunction) – The transfer function object.
  • expression (str, optional) – Format of the foton expression. Choose from [“zpk”, “rpoly”]. Defaults to “zpk”.
  • root_location (str, optional) – Root location of the zeros and poles for expression==”zpk”. Choose from [“s”, “f”, “n”]. “s”: roots in s-plane, i.e. zpk([…], […], …, “s”). “f”: roots in frequency plane, i.e. zpk([…], [,,,], …, “f”). “n”: roots in frequency plane but negated and gains are normalized, i.e. real parts are positive zpk([…], […], …, “n”). Defaults to “s”.
  • significant_figures (int, optional) – Number of significant figures to print out. Defaults to 6.
  • itol (float, optional) – Treating complex roots as real roots if the ratio of the imaginary part and the real part is smaller than this tolerance Defaults to 1e-25.
  • epsilon (float, optional) – Small number to add to denominator to prevent division error. Defaults to 1e-25.
Returns:

foton_expression – The foton expression in selected format.
Return type:

str
kontrol.core.foton.tf2rpoly(tf, significant_figures=6)

Convert a transfer function to foton rpoly expression.

Parameters:
  • tf (TransferFunction) – The transfer function object
  • significant_figures (int, optional) – Number of significant figures to print out. Defaults to 6.
Returns:

Foton express in foton rpoly expression.
Return type:

str

Notes

Only works for transfer functions with less than 20 orders.

kontrol.core.foton.tf2zpk(tf, root_location='s', significant_figures=6, itol=1e-25, epsilon=1e-25)

Convert a single transfer function to foton zpk expression.

Parameters:
  • tf (TransferFunction) – The transfer function object.
  • root_location (str, optional) – Root location of the zeros and poles. Choose from [“s”, “f”, “n”]. “s”: roots in s-plane, i.e. zpk([…], […], …, “s”). “f”: roots in frequency plane, i.e. zpk([…], [,,,], …, “f”). “n”: roots in frequency plane but negated and gains are normalized, i.e. real parts are positive zpk([…], […], …, “n”). Defaults to “s”.
  • significant_figures (int, optional) – Number of significant figures to print out. Defaults to 6.
  • itol (float, optional) – Treating complex roots as real roots if the ratio of the imaginary part and the real part is smaller than this tolerance Defaults to 1e-25.
  • epsilon (float, optional) – Small number to add to denominator to prevent division error. Defaults to 1e-25.
Returns:

The foton zpk expression in selected format.
Return type:

str

Notes

Only works for transfer functions with less than 20 orders.

kontrol.core.foton.zpk2tf(foton_string)

Convert a Foton ZPK string to TransferFunction

foton_string : str
The Foton ZPK string, e.g. zpk([0], [1; 1], 1, “n”).
Returns:tf – The transfer function.
Return type:TransferFunction

Curve Fitting

Curve fitting

Curve fitting class

Base class for general curve fitting.

class kontrol.curvefit.CurveFit(xdata=None, ydata=None, model=None, model_kwargs=None, cost=None, optimizer=None, optimizer_kwargs=None)

Bases: object

Base class for curve fitting

Parameters:
  • xdata (array, ydata: array) -> float) – The independent variable data / data on the x axis. Defaults to None.
  • ydata (array or None, optional) – The dependent variable data / data on the y axis. Defaults to None.
  • model (func(x: array, args: array, **kwargs) -> array, or None, optional) – The model used to fit the data. args in model is an array of parameters that define the model. Defaults to None
  • model_kwargs (dict or None, optional) – Keyword arguments passed to the model. Defaults to None.
  • cost (kontrol.curvefit.Cost or func(args, model, xdata, ydata) -> array) – Cost function.
  • xdata – The cost function to be used to fit the data. First argument is a list of parameters that will be passed to the model. This must be pickleable if multiprocessing is to be used. Defaults to None.
  • optimizer (func(func, **kwargs) -> OptimizeResult, or None, optional) – The optimization algorithm use for minimizing the cost function.
  • optimizer_kwargs (dict or None, optional) – Keyword arguments passed to the optimizer function. Defaults to None.
cost

The cost function to be used to fit the data.

fit(model_kwargs=None, cost_kwargs=None, optimizer_kwargs=None)

Fit the data

Sets self.optimized_args and self.optimize_result.

Parameters:
  • model_kwargs (dict or None, optional) – Overriding keyword arguments in self.model_kwargs. Defaults to None.
  • cost_kwargs (dict or None, optional) – Overriding keyword arguments in self.cost_kwargs. Defaults to None.
  • optimizer_kwargs (dict or None, optional) – Overriding keyword argumetns in self.optimizer_kwargs. Defaults to None.
Returns:
Return type:

scipy.optimizer.OptimizeResult
model

The model used to fit the data.

model_kwargs

Keyword arguments passed to the model.

optimize_result

Optimization Result

optimized_args

The optimized arguments

optimizer

The optimizer function to be used to fit the data.

optimizer_kwargs

Keyword arguments passed to the optimizer function.

prefit()

Something to do before fitting.

xdata

The independent variable data.

ydata

The dependent variable data

yfit

The fitted y values

Transfer function fitting class

Curve fitting class for transfer function fitting.

class kontrol.curvefit.TransferFunctionFit(xdata=None, ydata=None, model=None, model_kwargs=None, cost=None, weight=None, error_func_kwargs=None, optimizer=None, optimizer_kwargs=None, options=None, x0=None)

Bases: kontrol.curvefit.curvefit.CurveFit

Transfer function fitting class

This class is basically a CurveFit class with default cost functions and optimizer that is designed for fitting a transfer function. By default, the error function is kontrol.curvefit.error_func.tf_error(). The default optimizer is scipy.optimize.minimize( ...,method="Nelder-Mead",...), with options = {"adaptive": True, "maxiter": N*1000}, where N is the number of variables. All of these can be overridden if specified.

Parameters:
  • xdata (array or None, optional) – Independent variable data. Defaults to None.
  • ydata (array or None, optional) – Transfer function frequency response in complex numbers. Defaults to None.
  • model (func(x: ``array, args: array, **kwargs)`` -> array, or None, optional) – The model used to fit the data. args in model is an array of parameters that define the model. Defaults to None
  • model_kwargs (dict or None, optional) – Keyword arguments passed to the model. Defaults to None.
  • cost (kontrol.curvefit.Cost or func(args, model, xdata, ydata) -> array, optional) – Cost function. The cost function to be used to fit the data. First argument is a list of parameters that will be passed to the model. This must be pickleable if multiprocessing is to be used. Defaults to None.
  • weight (array or None, optional) – Weighting function. Defaults None.
  • error_func_kwargs (dict or None, optional) – Keyword arguments the will be passed to error_func, which is passed to the construct the cost function. Defaults to None.
  • optimizer (func(func, **kwargs) -> OptimizeResult, or None, optional) – The optimization algorithm use for minimizing the cost function.
  • optimizer_kwargs (dict or None, optional) – Keyword arguments passed to the optimizer function. Defaults to None.
  • options (dict or None, optional) – The option arguments passed to the optimizer
  • x0 (array, optional) – Inital guess. Defaults to None.
options

Option arguments passed to optimizer

weight

Weighting function

x0

Initial guess

Models for Curve Fitting

These classes are designed to use with kontrol.curvefit

class kontrol.curvefit.model.Model(args=None, nargs=None, log_args=False)

Bases: object

Model base class for curve fitting

args

Model parameters

nargs

Number of model parameters

class kontrol.curvefit.model.TransferFunctionModel(nzero, npole, args=None, log_args=False)

Bases: kontrol.curvefit.model.model.Model

Transfer function model class defined by numerator and denominator

Parameters:
  • nzero (int) – Number of zeros.
  • npole (int) – Number of poles.
  • args (array or None, optional.) – The model parameters. Structured as follows: [b_n, b_n-1,…, b_1, b_0, a_m, a_m-1,…, a_1, a_0], where b and a are the coefficients of the numerator and denominator respectively, ordered from higher-order to lower-order. Defaults to None.
tf

The last evaluted transfer function.

Type:kontrol.TransferFunction

Notes

The transfer function model is defined as

\[G(s, b_n, b_{n-1},..., b_1, b_0, a_m, a_{m-1},..., a_1, a_0) = \frac{\sum_{i=0}^{n} b_i s^i}{\sum_{j=0}^{m} a_j s^j}\]
den

Denominator array

npole

Number of zeros

num

Numerator array

nzero

Number of zeros

tf

The Transfer Function object.

class kontrol.curvefit.model.DampedOscillator(args=None)

Bases: kontrol.curvefit.model.transfer_function_model.TransferFunctionModel

Transfer function model for a damped oscillator.

Parameters:args (array or None, optional.) – The model parameters with three numbers. Structured as follows: args = [k, fn, q], where k is the DC gain of the transfer function, fn is the resonance frequency in Hz, and q is the Q-factor. Defaults to None.

Notes

The model is definded as

\[G(s; k, \omega_n, q) = k\frac{\omega_n^2}{s^2 + \frac{\omega_n}{q}s + \omega_n^2}\]

where \(k\) is the DC gain of the transfer function, \(\omega_n\) is the resonance frequency of the oscillator, and \(q\) is the Q-factor if the damped oscillator.

args

Model parameters

damped_oscillator_args

The model parameters with three numbers [k, fn, q]

fn

Resonance frequency

gain

DC gain

q

Q factor

class kontrol.curvefit.model.SimpleZPK(nzero, npole, args=None, log_args=False)

Bases: kontrol.curvefit.model.model.Model

ZPK model with simple poles and zeros.

Parameters:
  • nzero (int) – Number of simple zeros.
  • npole (int) – Number of simple poles.
  • args (array, optional) – The model parameters defined as args = [z1, z2,..., p1, p2,..., k] where z are the locations of the zeros in Hz, p are the locations of the pole in Hz, and k is the static gain,
  • log_args (boolean, optional) – If true, model parameters passed to the model are assumed to be passed through a log10() function. So, when the real parameters will be assumed to be 10**args instead. Defaults to False.
zero

List of zeros.

Type:array
pole

List of poles.

Type:array
gain

Static gain.

Type:float
tf

The transfer function representation of this ZPK model

Type:kontrol.TranserFunction

Notes

The simple ZPK model is defined as

\[G(s; z_1, z_2,..., p_1, p_2, ..., k) = k\frac{\prod_i(\frac{s}{2\pi z_i}+1)}{\prod_j(\frac{s}{2\pi p_j}+1)}\]
args

Model parameters in ZPK [z1, z2,…, p1, p2,…, k] format

gain

Static gain

npole

Number of complex pole pairs

nzero

Number of simple zeros

pole

List of poles

tf

Returns a TransferFunction object of this ZPK model

zero

List of zeros

class kontrol.curvefit.model.ComplexZPK(nzero_pairs, npole_pairs, args=None, log_args=False)

Bases: kontrol.curvefit.model.model.Model

ZPK model with complex poles and zeros.

Parameters:
  • nzero_pairs (int) – Number of complex zero pairs.
  • npole_pairs (int) – Number of complex pole pairs.
  • args (array, optional) – The model parameters defined as args = [f1, q1, f2, q2,..., fi, qi,..., fn, qn, k], where f are resonance frequencies of the complex zero/pole pairs, q are the quality factors, and k is the static gain, i is the number of complex zero pairs, and n-i is the number of of complex pole pairs.
  • log_args (boolean, optional) – If true, model parameters passed to the model are assumed to be passed through a log10() function. So, when the real parameters will be assumed to be 10**args instead. Defaults to False.
fn_zero

List of resonance frequencies of the complex zeros.

Type:array
fn_pole

List of resonance frequencies of the complex poles.

Type:array
q_zero

List of Q-factors of the complex zeros.

Type:array
q_pole

List of Q-factors of the complex poles.

Type:array

Notes

The complex ZPK model is defined by:

\[G(s; f_i, q_i,k) =k\frac{\prod_i(\frac{s^2}{(2\pi f_i)^2} + \frac{1}{2\pi f_i q_i}s + 1)} {\prod_j(\frac{s^2}{(2\pi f_j)^2} + \frac{1}{2\pi f_j q_j}s + 1)}\]
args

Model parameters in ZPK [f1, q1, f2, q2,…, fn, qn, k] format

fn_pole

List of resonance frequencies of complex poles.

fn_zero

List of resonance frequencies of complex zeros.

gain

Static gain.

npole_pairs

Number of complex pole pairs

nzero_pairs

Number of complex zero pairs

q_pole

List of quality factors of the complex poles

q_zero

List of quality factors of the complex zeros

tf

Returns a TransferFunction object of this ZPK model

class kontrol.curvefit.model.StraightLine(args=None)

Bases: kontrol.curvefit.model.model.Model

Simply a straight line.

Parameters:args (array, optional) – The model parameters structued as: [slope, y-intercept]. Defaults to None.
slope

The slope of the line.

Type:float
intercept

The y-intercept of the line.

Type:float

Notes

The straight line is defined as

\[y(x; m, c) = mx + c\,,\]

where \(m\) is the slope and \(c\) is the y-intercept.

intercept

The y-intercept of the line

slope

The slope of the line

class kontrol.curvefit.model.Erf(args=None)

Bases: kontrol.curvefit.model.model.Model

The error function erf(x)

Parameters:args (array, optional) – The model parameters structured as [amplitude, slope, x0, y0]. Defaults to None.

Notes

The function is defined as

\[f(x; a, m, x_0, y_0) = a\,\mathrm{erf}(m(x-x_0)) + y_0\,,\]

where \(\mathrm{erf}(x)\) is the error function [1]_, \(a\) is the peak to peak amplitude, \(m\) is the slope at the inflection point, \(x_0\) is the \(x\) offset , and \(y_0\) is the \(y\) offset.

References

[1]https://en.wikipedia.org/wiki/Error_function
amplitude

Peak to Peak amplitude of the error function.

slope

Slope at the inflection point.

x_offset

x offset

y_offset

y offset

Control Regulator Design

Functions for algorithmic design of control regulators

Feedback Control

Algorithmic designs for feedback control regulators.

kontrol.regulator.feedback.add_integral_control(plant, regulator=None, integrator_ugf=None, integrator_time_constant=None, **kwargs)

Match and returns an integral gain.

This function finds an integral gain such that the UGF of the integral control matches that of the specified regulator. If integrator_ugf or integrator_time_constant is specified instead, these will be matched instead.

Parameters:
  • plant (TransferFunction) – The transfer function representation of the system to be feedback controlled.
  • regulator (TransferFunction, optional) – The pre-regulator Use kontrol.regulator.feedback.proportional_derivative() or kontrol.regulator.feedback.critical_damping() to make one for oscillator-like systems.
  • integrator_ugf (float, optional) – The unity gain frequency (Hz) of the integral control. This is the inverse of the integration time constant. If integrator_time_constant is not None, then this value will be ignored. If set to None, it’ll be set to match the first UGF of the derivative control. Defaults to None.
  • integrator_time_constant (float, optional,) – The integration time constant (s) for integral control. Setting this will override the integrator_ugf argument. Defaults to None.
Returns:

ki – The integral control gain.
Return type:

float
kontrol.regulator.feedback.add_proportional_control(plant, regulator=None, dcgain=None, **kwargs)

Match and returns proportional gain.

This function finds a proportional gain such that the proportional control UGF matches the first UGF of the specified regulator (typically a derivative control regulator). If dcgain is specified, the DC gain is matched instead.

Parameters:
  • plant (TransferFunction) – The transfer function representation of the system to be feedback controlled. The plant must contain at least one pair of complex poles.
  • regulator (TransferFunction, optional) – The pre-regulator. If not specified, then dcgain must be specified. Defaults to None.
  • dcgain (float, optional) – The desired DC gain of the open-loop transfer function. If not specified, the portional gain is tuned such that the portional control’s first UGF matches that of the derivative control. Defaults to None.
Returns:

kp – The proportional control gain.
Return type:

float
kontrol.regulator.feedback.critical_damp_calculated(plant, nmode=1, **kwargs)

Find dominant mode and returns the approximate critical regulator.

Parameters:
  • plant (TransferFunction) – The transfer function representation of the system to be feedback controlled. The plant must contain at least one pair of complex poles.
  • nmode (int, optional) – The ``nmode``th dominant mode to be damped. This number must be less than the number of modes in the plant. Defaults 1.
Returns:

kd – The derivative gain for critically damping the dominant mode.
Return type:

float

Notes

The plant must contain at least one complex pole pair. The plant must have a non-zero DC gain.

The critical regulator is approximated by

\[K(s) \approx \frac{2}{\omega_n K_\mathrm{DC}}\,,\]

where \(\omega_n\) is the resonance frequency of the dominant mode and \(K_\mathrm{DC}\) is the sum of DC gains of the mode and that of the other higher-frequency modes..

kontrol.regulator.feedback.critical_damp_optimize(plant, gain_step=1.1, ktol=1e-06, **kwargs)

Optimize derivative damping gain and returns the critical regulator

Parameters:
  • plant (TransferFunction) – The transfer function representation of the system to be feedback controlled. The plant must contain at least one pair of complex poles.
  • gain_step (float, optional) – The multiplicative factor of the gain for finding the gain upper bound. It must be greater than 1. Defaults to 1.1.
  • ktol (float, optional) – The tolerance for the convergence condition. The convergence condition is (kd_max-kd_min)/kd_min > ktol. It must be greater than 0. Defaults to 1e-6.
Returns:

kd – The derivative gain for critically damping the dominant mode.
Return type:

float

Notes

Update on 2021-12-04: Use carefully. It only critically damps the mode that has the highest peak when multiplied by an differentiator. Note for myself: only 2 complex poles can become simple poles for plants with second-order rolloff.

Only works with plants that contain at least one complex pole pair. Works best with plants that only contain complex zeros/poles. If it returns unreasonably high gain, try lowering gain_step.

The algorithm goes as follows.

1. Find the minimum damping gain k_min such that the open-loop transfer function k_min*s*plant has maxmimum gain at unity gain frequency.

  1. Iterate i: k_i=k_min*i*gain_step for i=1,2,3…

3. Terminate when 1/(1+k_i*s*plant) has less complex pole pairs than the plant itself, i.e. one mode has been overdamped. Then, define k_max=k_i.

4. Iterate: Define k_mid as the logarithmic mean of k_max and k_min. If 1/(1+k_mid*s*plant) has less complex pole pairs than plant, i.e. overdamps, then set k_max=k_mid. Otherwise, set k_min=k_mid.

  1. Terminate when (k_max - k_min)/k_min < ktol, i.e. converges.
kontrol.regulator.feedback.critical_damping(plant, method='calculated', **kwargs)

Returns the critical damping derivative control gain.

This functions calls kontrol.regulator.feedback.critical_damp_calculated() or kontrol.regulator.feedback.cricical_damp_optimized() and returns the derivative control gain.

Parameters:
  • plant (TransferFunction) – The transfer function representation of the system to be feedback controlled. The plant must contain at least one pair of complex poles.
  • method (str, optional) –

    The method to be used for setting the gain. Choose from [“optimized”, “calculated”].

    ”optimized”: the gain is optimized until the dominant complex pole pairs become two simple poles.

    ”calculated”: the gain is set to \(2/\omega_n/K_{DC}\), where \(\omega_n\) is the resonance frequency in rad/s of the mode to be damped, and \(K_{DC}\) is the sum of DC gains of the mode and that of the other high-frequency modes.

    Both method assumes that the plant has at least one pair of complex poles.

    Defaults to “calculated”.

  • **kwargs (dict) –

    Method specific keyword arguments.

    See:

    • ”optimized”: kontrol.regulator.feedback.critical_damp_optimized
    • ”calculated”: kontrol.regulator.feedback.critical_damp_calculated
Returns:

kd – The derivative gain for critically damping the dominant mode.
Return type:

float
kontrol.regulator.feedback.mode_composition(wn, q, k)

Create a plant composed of many modes.

Parameters:
  • wn (array) –
  • (rad/s) (Frequencies) –
  • q (array) – Q factors.
  • k (array) – Dcgains of the modes.
Returns:

The composed plant.
Return type:

TransferFunction
kontrol.regulator.feedback.mode_decomposition(plant)

Returns a list of single mode transfer functions

Parameters:plant (TransferFunction) – The transfer function with at list one pair of complex poles.
Returns:
  • wn (array) – Frequencies (rad/s).
  • q (array) – Q factors.
  • k (array) – Dcgains of the modes.

Regulator for Oscillatory Systems

Control regulators designs for oscillator-like systems

kontrol.regulator.oscillator.pid(plant, regulator_type='PID', dcgain=None, integrator_ugf=None, integrator_time_constant=None, return_gain=False, **kwargs)

PID-like controller design for oscillator-like systems

Parameters:
  • plant (TransferFunction) – The transfer function of the system that needs to be controlled.
  • regulator_type (str, optional) – The type of the contorl regulator. Choose from {"PID", "PD", "PI", "I", "D"} for proportional-integral-derivative, proportional-derivative, proportional-integral, or derivative (velocity) control respectively. Defaults to “PID”.
  • dcgain (float, optional) – The DC gain of the OLTF of the proportional control. If set to None, it will be set automatically depending on the type of the controller. If regulator_type=="PI", then dcgain will be set such that the proportional control UGF matches that of the integral control. If regulator_type=="PD" or "PID", then the dcgain will be set such that it matches the UGF of the derivative control. Defaults to None.
  • integrator_ugf (float, optional) – The unity gain frequency (Hz) of the integral control. This is the inverse of the integration time constant. If integrator_time_constant is not None, then this value will be ignored. If set to None, it’ll be set to match the UGF of the derivative control. For regulator_type=="I", this must be specified. Defaults to None.
  • integrator_time_constant (float, optional,) – The integration time constant (s) for integral control. Setting this will override the integrator_ugf argument. Defaults to None.
  • return_gain (boolean, optional) – Return the PID gain instead.
Returns:

  • regulator (TransferFunction, optional) – The regulator. Return only if return_gain is False.
  • kp (float, optional) – Proportional gain. Return only if return_gain is ``True’’.
  • ki (float, optional) – Integral gain. Return only if return_gain is True.
  • kd (float, optional) – Derivative gain. Return only if return_gain is True.

Post Filtering

Functions for designing post-regulator filters

kontrol.regulator.post_filter.post_low_pass(plant, regulator, post_filter=None, ignore_ugf_above=None, decades_after_ugf=1, phase_margin=45, f_start=None, f_step=1.1, low_pass=None, mtol=1e-06, small_number=1e-06, oscillatory=True, **kwargs)

Add low-pass filter after regulator.

This function lowers/increase the the cutoff frequency of a low-pass filter until the phase margin at a dedicated ugf crosses the specified phase margin. Then, runs a bisection algorithm to poolish the cutoff frequency until the phase margin converges relative to the specified tolerance.

Parameters:
  • plant (TransferFunction) – The transfer function of the system that needs to be controlled.
  • regulator (TransferFunction) – The regulator.
  • post_filter (TransferFunction, optional) – Any post filters that will be applied on top of the regulator. Defaults None.
  • ignore_ugf_above (float, optional) – Ignore unity gain frequencies higher than ignore_ugf_above (Hz). If not specified, defaults to 1 decade higher than the last UGF. This value can be overrided by the argument decades_after_ugf Note that there’s no guarantee that the UGF will be lower than this. The priority is to match the target phase margin. Defaults to None.
  • decades_after_ugf (float, optional) – Set ignore_ugf_above some decades higher than the UGF of the OLTF ignore_ugf_above is None. Defaults to 1.
  • phase_margin (float, optional,) – The target phase margin (Degrees). Defaults to 45.
  • f_start (float, optional,) – The cutoff frequency to start iterating with. If not specified, defaults to some decades higher than the highest UGF. “Some decade” is set by decades_after_ugf. Defaults None.
  • f_step (float, optional,) – The gain that is used to multiply (or divide) the cutoff frequency during a coarse search. Defaults 1.1
  • low_pass (func(cutoff, order) -> TransferFunction, optional) – The low-pass filter. If not specified, kontrol.regulator.predefined.lowpass() with order 2 will be used. Defaults to None.
  • mtol (float, optional) – Tolerance for convergence of phase margin. Defaults to 1e-6.
  • small_number (float, optional) – A small number as a delta f to detect whether the gain is a rising or lowering edge at the unity gain frequency. Defaults to 1e-6.
  • oscillatory (boolean, optional) – Use the first mode of the oscillatory system to evaluate the phase margins to avoid having UGFs at steep phase slopes. The benefit of using this is to have a conservative phase margin estimate. The phase response of the first mode is the lower bound of the phase response. If False, use the plant itself to calculate phase margins. If the plant does not contain any complex poles, this option will be overridden to False. Defaults True.
  • **kwargs – Keyword arguments passed to low_pass.
Returns:

The low-pass filter.
Return type:

TransferFunction
kontrol.regulator.post_filter.post_notch(plant, regulator=None, post_filter=None, target_gain=None, notch_peaks_above=None, phase_margin=45, notch=None, **kwargs)

Returns a list of notch filters that suppress resonance peaks.

This functions finds the resonances peak of the plant/OLTF above certain frequencies and returns a list of notch filters that suppress these peaks to the target gains.

Parameters:
  • plant (TransferFunction) – The transfer function of the system that needs to be controlled.
  • regulator (TransferFunction, optional) – The regulator. Defaults to None.
  • post_filter (TransferFunction, optional) – Any post filters that will be applied on top of the regulator. Defaults None.
  • target_gain (float, optional) – The target open-loop gain for the suppressed peak. To ensure a stable system, a value of less than 1 is recommended. If not specified, the notch will fully suppress the peak. Default None.
  • phase_margin (float, optional) – The target phase margin. Defaults to 45.
  • notch_peaks_above (float, optional) – Notch modes that has freqeuncies above notch_peaks_above. If not specified, defaults to the highest unity gain frequency that is above phase_margin Defaults to None.
  • notch (func(frequency, q, depth) -> TransferFunction, optional) – The notch filter. If not specified, kontrol.Notch() will be used. Defaults to None.
  • **kwargs – Keyword arguments passed to notch().
Returns:

A list of notch filters.
Return type:

list of TransferFunction

Notes

This operation does not guarantee stability. It only finds resonances peaking out of the unity gain above some unity gain frequency and make notch filters to suppress them to a target gain level lower than the unity gain. The stability of the notched OLTF is not checked whatsoever.

Predefined Filters

Predefined regulator library.

kontrol.regulator.predefined.low_pass(cutoff, order=1, **kwargs)

Simple low-pass filter

Parameters:
  • cutoff (float) – Cutoff frequency (Hz)
  • order (int, optional) – The order of the filter. Defaults to be 1.
  • **kwargs – Keyword arguments holder. Not passed to anywhere.
Returns:

The low-pass filter.
Return type:

TransferFunction

Notes

The low-pass filter is defined as

\[L(s) = \left(\frac{2\pi f_c}{s+2\pi f_c}\right)^n\,,\]

where \(f_c\) is the cutoff frequency (Hz), \(n\) is the order of the filter.

kontrol.regulator.predefined.notch(frequency, q, depth=None, depth_db=None, **kwargs)

Notch filter defined in Foton.

Parameters:
  • frequency (float) – The notch frequency (Hz).
  • q (float) – The quality factor.
  • depth (float, optional) – The depth of the notch filter (magnitude). If not specified, depth_db will be used. Defaults None.
  • depth_db (float, optional) – The depth of the notch filter (decibel). If not specified, depth will be used instead. Defaults None.
Returns:

The notch filter
Return type:

TransferFunction

Notes

The notch filter is defined by Foton, as

\[N(s) = \frac{s^2 + (2\pi f_n)/(dQ/2)s + (2\pi f_n)^2} {s^2 + (2\pi f_n)/(Q/2)s + (2\pi f_n)^2}\,,\]

where \(f_n\) is the notch frequency, \(q\) is the quality factor , and \(d\) is the depth.

kontrol.regulator.predefined.pid(kp=0, ki=0, kd=0)

Alias of proportional_integral_derivative()

kontrol.regulator.predefined.proportional_integral_derivative(kp=0, ki=0, kd=0)

PID control build on proportional_derivative().

Parameters:
  • kp (float, optional) – The proportional control gain. Defaults to 0.
  • ki (float, optional) – The integral control gain. Defaults to 0.
  • kd (float, optional) – Defaults to 0.
Returns:

The PID controller.
Return type:

TransferFunction

Notes

The PID controller is defined as

\[K_\mathrm{PID}(s) = K_p + K_i/s + K_d s\,.\]

Dynamic Mode Decomposition

Dynamic mode decomposition algorithm

Dynamic mode decomposition main class

Dynamic mode decomposition class

class kontrol.dmd.dmd.DMD(snapshot_1, snapshot_2=None, truncation_value=None, dt=1, run=False)

Bases: object

Dynamic mode decomposition class

Parameters:
  • snapshot_1 (array) – The 2-D snapshot array.
  • snapshot_2 (array, optional) – The 2-D snapshot array at next time step. If not provided, snapshot_1[:, :-1] becomes snapshot_1 and snapshot_2[:, 1:] becomes snapshot_2. Defaults None.
  • truncation_value (float, optional) – The truncation value (order/rank) of the system. Defaults None.
  • dt (float, optional) – The time difference between two snapshots. Defaults 1.
  • run (bool, optional) – Run dynamic mode decomposition upon construction. Computes DMD modes, Reduced-order model, etc. truncation_value must be specified for this to work. Defaults False.
snapshot_1

The 2-D snapshot matrix

Type:array
snapshot_2

The 2-D snapshot matrix at next time step.

Type:array
truncation_value

The truncation value (order/rank) of the system.

Type:int
dt

The time difference between two snapshots.

Type:float
u

The left singular vector of the snapshot_1 matrix.

Type:array
vh

The right singular vector (conjugate-transposed) of the snapshot_1 matrix.

Type:array
sigma

The singular values of the snapshot_1 matrix

Type:array
u_truncated

The truncated left singular vector of the snapshot_1 matrix.

Type:array
vh_truncated

The truncated right singular vector (conjugate-transposed) of the snapshot_1 matrix.

Type:array
sigma_truncated

The truncated singular values of the snapshot_1 matrix

Type:array
A_reduced

The reduced-order model of the dynamics.

Type:array
w_reduced

The eigenvalues of the reduced-order model.

Type:array
v_reduced

The eigenvectors of the reduced-order model.

Type:array
dmd_modes

The DMD modes.

Type:array
complex_frequencies

The complex frequencies of the modes.

Type:array
v_constant

The constant vector of the time dynamics, determined by initial conditions.

Type:array
time_dynamics

The time dynamics of the ODE solution.

Type:array
prediction

The predicted time series.

Type:array
A_reduced

Reduced-order model

complex_frequencies

Complex frequencies

compute_complex_frequencies(dt=None)

Compute the complex frequencies

Parameters:dt (float, optional) – The time spacing between snapshots. Specified as self.dt or in constructor option. Defaults 1.
Returns:complex_frequencies – Array of complex frequencies
Return type:array
compute_dmd_modes()

Compute the DMD modes

Returns:dmd_modes – The DMD modes.
Return type:array
compute_reduced_model()

Compute the reduced-order model

Returns:A_reduced – Matrix representing the reduced-order model.
Return type:array
dt

Time difference between the snapshots

eig_reduced_model()

Eigen decomposition of the reduced-order model

Returns:
  • w_reduced (array) – The list of eigenvalues.
  • v_reduced (array) – Eigenvalues as columns.
low_rank_approximation(truncation_value=None)

Truncate the svd.

truncation_value : int, optional
The truncation value (order/rank) of the system. Specify as self.truncation_value or via the constuctor option. Defaults None.
Returns:
  • u_truncated (array) – Truncated unitary matrix having left singular vectors as columns.
  • sigma_truncated (array) – Truncated singular values.
  • vh_truncated (array) – Truncated unitary matrix having right singular vectors as rows.
predict(t, t_constant=None, i_constant=0)

Predict the future states given a time array

Parameters:
  • t (array) – The time array
  • t_constant (float, optional) – Time at which the constant vector is calculated. Defaults to be t[0] if not given.
  • i_constant (int, optional) – Index of the snapshot_1 at which the constant vector is calculated. Defaults 0.
Returns:

prediction – Predicted states.
Return type:

array

Notes

The constant vector is the vector that evolves and is evaluated using the initial/boundary values. Set t_constant and i_constant to where the prediction of the time series begins.

prediction

The ODE solution

run()

Run the DMD algorithm, compute dmd modes and complex frequencies

sigma

Singular values

sigma_truncated

Singular values truncated

snapshot_1

Snapshot 1

snapshot_2

Snapshot 1

svd()

Decompose the snapshot 1

Returns:
  • u (array) – Unitary matrix having left singular vectors as columns.
  • sigma (array) – The singular values.
  • vh (array) – Unitary matrix having right singular vectors as rows.
time_dynamics

Time dynamics of the ODE solution

truncation_value

Truncation value (order/rank) of the system

u

Left singular vectors

u_truncated

Left singular vectors (truncated)

v_constant

Constant vector of the ODE solution

vh

Right singular vectors

vh_truncated

Right singular vectors (truncated)

Dynamic mode decomposition utilities

Utility functions for Dymanic Mode Decomposition.

kontrol.dmd.utils.auto_truncate(sigma, threshold=0.99)

Automatically get truncation value

Parameters:
  • sigma (array) – Singular values. Arranged from large to small values.
  • threshold (float, optional) – Only include singular values so their sum is threshold of the original sum. Defaults 0.99.
Returns:

truncation_value – Amount of singular values that are required to reach the threshold sum.
Return type:

int
kontrol.dmd.utils.hankel(array, order)

Hankelize an array

Parameters:
  • array (array) – Array with rows as channels and columns as samples
  • order (int) – The number of rows of the hankelized matrix.
Returns:

hankel_array – The hankelized array
Return type:

array

Examples