visualdynamics.core.drive_prediction

drive_prediction

The drives a specification asks for, predicted from an FRF before the test, and what they ask of the amplifiers and the shakers.

A random test is controlled on spectra. The control channels' spectral density matrix is the drives' passed through the plant,

S_yy = conj(H) S_vv H^T,

at every line, where H is the FRF matrix from the drives to the control channels (one row per control channel, one column per drive, in g/V or (m/s^2)/V), S_vv the drives' spectral density matrix and H^T the transpose of H. The conjugate is on the left because a cross term between channels a and b is stored as conj(X_a) X_b (TimeHistory.compute_cpsds), so outputs Y = H V have the matrix E[conj(Y) Y^T] = conj(H) E[conj(V) V^T] H^T. Written the other way round, H S_vv H^H gives other autospectra wherever the plant and the cross terms are both complex. Turned round, the drives that would produce a target S_ref are

S_vv = conj(Z) S_ref Z^T,

with Z an inverse of H. On one drive and one control channel the matrices are single numbers and this is the drive autospectrum a controller starts from, G_vv = G_ref / |H|^2. With several, H is inverted line by line, and where it is ill-conditioned (two drives that reach the control channels in nearly the same proportions, or a line where the plant barely responds) the plain inverse is made of noise, so it is regularized, in one of two stated ways:

• Truncated (the default): the singular values of H at a line below tolerance times that line's largest are dropped, and the pseudo-inverse is built from the rest. With the default tolerance of one hundredth, a line whose condition number is over 100 loses its weakest direction, the target there is met only in the directions kept, and the reading says so. On one drive and one channel nothing is ever dropped, and the inverse is exactly 1 / H.
• Tikhonov: Z = (H^H H + lambda I)^(-1) H^H, with lambda the stated weight times the largest squared singular value of H anywhere in the band, so the weight has no units and means the same for any plant. Where |H|^2 is large against lambda this is the inverse; where the plant barely responds it stays small instead of growing without limit. On one channel it is H* / (|H|^2 + lambda).

Either way the prediction carries the control channels' response to the drives it predicts, conj(H) S_vv H^T, and the lines where that misses the target by more than UNMET_DB are the lines the regularization gave up, named as bands. Any other channel the FRF holds against the same drives (a force gauge, a current monitor, an accelerometer on the article) is predicted the same way, its autospectrum the diagonal of conj(H_u) S_vv H_u^T.

The prediction is open loop, made once, before anything is played: it is what a controller's first drive would be with this FRF as its plant estimate, and it does not iterate.

What it asks of the rig

Every quantity is a Gaussian random signal with a predicted autospectrum, so its rms is the root of the area under it (compliance.rms) and its expected largest peak over a run of duration seconds is the rms times the Gaussian peak factor sqrt(2 ln(2 nu0 T)), with nu0 its own crossing rate by Rice's formula (time_statistics.rice_rate and expected_largest_peak), counted both ways: a rating is reached by a peak of either sign, the drive's limit and the travel each way alike, so the largest of either sign is the one judged. A stated peak_factor replaces that, for a controller that clips its drive at a set multiple of the rms. Each quantity is judged on its peak:

• the drive at each drive, against drive_limit, in the drive's unit (the volts a source or amplifier input accepts, or the newtons of a force-controlled drive);
• for a drive given a ShakerRating, its shaker, read at its table channel (the acceleration measured at the drive's own DOF, or the one control channel of a single-input test): the force, the moving mass times the peak acceleration; the velocity, from the table's acceleration autospectrum over (2 pi f)^2; and the travel each way from the middle of the stroke, from the autospectrum over (2 pi f)^4, with the body's motion counted. All three go through core.shaker's one rule;
• any other channel given a peak limits entry, in that channel's unit: a force gauge in newtons, a current monitor in the volts its scale gives at the amplifier's current rating.

Each judged quantity names the line that contributes most to its mean square and the band of lines holding half of it, so the reading can say which drive or shaker runs out first and at what frequency.

References

The spectral relations between the inputs and outputs of a linear system with several inputs, and the Gaussian peak statistics, are Bendat and Piersol's; the singular value decomposition, the pseudo-inverse and its truncation are Golub and Van Loan's; the regularized inverse is Tikhonov's, in the form Hansen sets out.

1. Bendat, J. S., & Piersol, A. G. (2010). Random Data: Analysis and Measurement Procedures, 4th ed. Wiley. The output spectral density of a multiple-input linear system and the expected largest peak of a stationary Gaussian process.
2. Golub, G. H., & Van Loan, C. F. (2013). Matrix Computations, 4th ed. Johns Hopkins University Press. The singular value decomposition, the pseudo-inverse and its truncation at a relative tolerance.
3. Hansen, P. C. (1998). Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion. SIAM. Tikhonov regularization as a filter on the singular values, sigma / (sigma^2 + lambda), beside truncation.

Classes:

Name	Description
Demand	One quantity the predicted drives ask of the rig, against its
DrivePrediction	The drives a specification asks for through an FRF, the

Functions:

Name	Description
regularized_inverse	An inverse of the plant at every line, truncated or Tikhonov.
inverse_from_svd	regularized_inverse from the plant's singular value
through	The plant's conjugate times density times the plant's
predict_drives	The drives that would make the specification at its control

Classes

Demand dataclass

Demand(kind: str, channel: str, dimension: str, rms: float, peak: float, rating: float | None, share: float | None, frequency: float, band: tuple[float, float])

One quantity the predicted drives ask of the rig, against its rating.

kind is 'drive', 'force', 'velocity', 'stroke' (the travel each way from the middle of the stroke) or 'response' (another channel against its own limit). channel names the drive or the channel it is read at, dimension the quantity of rms, peak and rating, all in SI. For the force, the velocity and the travel, peak is what core.shaker judges: the force at the peak acceleration, the velocity and the travel against the body. share is the peak over the rating, by core.shaker's rule for a shaker, None where no rating was given. frequency is the line contributing most to the mean square and band the span of the lines holding half of it.

Methods:

Name	Description
describe	The demand in a sentence, in unit_system's units.

Attributes:

Name	Type	Description
exceeds	bool | None	Whether the peak is over the rating; None without one.
name	str	The quantity and where, in words: 'travel each way at 1Z+'.

Attributes

exceeds property

exceeds: bool | None

Whether the peak is over the rating; None without one.

name property

name: str

The quantity and where, in words: 'travel each way at 1Z+'.

Methods:

describe

describe(unit_system: UnitSystem | None = None) -> str

The demand in a sentence, in unit_system's units.

Source code in src/visualdynamics/core/drive_prediction.py

def describe(self, unit_system: UnitSystem | None = None) -> str:
    """The demand in a sentence, in `unit_system`'s units."""
    if self.kind == 'stroke':
        words = (f'{self.name}: {_amount(self.peak, self.dimension, unit_system)}'
                 f' expected')
    else:
        words = (f'{self.name}: {_amount(self.rms, self.dimension, unit_system)}'
                 f' rms, {_amount(self.peak, self.dimension, unit_system)}'
                 ' expected peak')
    if self.rating is not None and self.share is not None:
        words += (f' against {_amount(self.rating, self.dimension, unit_system)}'
                  f' ({100.0 * self.share:.0f} %'
                  f'{", over" if self.share > 1.0 else ""})')
    else:
        words += ', no rating given'
    if np.isfinite(self.frequency):
        low, high = self.band
        words += (f'; most at {self.frequency:.4g} Hz, half from '
                  f'{low:.4g} to {high:.4g} Hz')
    return words

DrivePrediction dataclass

DrivePrediction(frequency: ndarray, drives: list[str], drive_quantity: str, rows: list[str], row_quantities: list[str], controls: list[int], drive_csd: ndarray, asd: ndarray, target: ndarray, control: ndarray, judged: ndarray, kept: ndarray, condition: ndarray, method: str, tolerance: float, weight: float, duration: float, peak_factor: float | None, demands: list[Demand] = list(), notes: list[str] = list(), row_channels: list[tuple[str, str]] = list(), drive_channels: list[tuple[str, str]] = list(), plain_notes: int = 0)

The drives a specification asks for through an FRF, the responses they make, and what they ask of the rig.

Everything is in SI, on the FRF's lines. drive_csd is the drives' spectral density matrix (lines, drives, drives) and asd the predicted autospectrum of every FRF row (lines, rows); target is the specification's matrix on the control channels and control the control channels' predicted matrix. judged marks the lines predicted: above zero hertz, the plant finite and every control autospectrum written; the arrays are not-a-number elsewhere. kept is how many directions the inverse kept at each line and condition the plant's condition number there. demands are the quantities judged, notes what the prediction should be read with.

Methods:

Name	Description
with_ratings	This prediction judged against other ratings, the drives and
drive_asd	One drive's predicted autospectrum, by index or label.
drive_coherence	The coherence between two predicted drives at every line:
missed	The lines where a control channel's prediction misses the
describe	The prediction in sentences: the verdict, each demand, the

Attributes:

Name	Type	Description
single	bool	One drive and one control channel.
shortfall_db	ndarray	The control channels' predicted autospectra over the target,
unmet	ndarray	The lines where some control channel's prediction misses the
unmet_bands	list[tuple[float, float]]	(first, last) frequency of each run of unmet lines.
judged_demands	list[Demand]	The demands with a rating, most of their rating first.
first_out	Demand | None	The demand nearest its rating, or furthest past it: what runs
fits	bool | None	Whether every rated demand is within its rating; None when
headroom_db	float | None	How far the level could rise before the first demand reaches
title	str	What was predicted, in a few words.

Attributes

single property

single: bool

One drive and one control channel.

shortfall_db property

shortfall_db: ndarray

The control channels' predicted autospectra over the target, in dB, (lines, controls): zero where the inverse met it.

unmet property

unmet: ndarray

The lines where some control channel's prediction misses the target by more than UNMET_DB (missed).

unmet_bands property

unmet_bands: list[tuple[float, float]]

(first, last) frequency of each run of unmet lines.

judged_demands property

judged_demands: list[Demand]

The demands with a rating, most of their rating first.

first_out property

first_out: Demand | None

The demand nearest its rating, or furthest past it: what runs out first as the level comes up. None when nothing is rated.

fits property

fits: bool | None

Whether every rated demand is within its rating; None when nothing is rated.

headroom_db property

headroom_db: float | None

How far the level could rise before the first demand reaches its rating, in dB on the test level (twenty times the log of one over its share: every peak scales with the level); negative when it is already past. None when nothing is rated.

title property

title: str

What was predicted, in a few words.

Methods:

with_ratings

with_ratings(*, drive_limit: float | Mapping[str, float] | None = None, shakers: ShakerRating | Mapping[str, ShakerRating] | None = None, tables: Mapping[str, str] | None = None, limits: Mapping[str, float] | None = None, duration: float = DURATION, peak_factor: float | None = None) -> DrivePrediction

This prediction judged against other ratings, the drives and responses kept: the inverse does not depend on the ratings, so a rating moved does not predict again. The ratings are taken exactly as predict_drives takes them, and one not given here is not judged, whatever the prediction was made with.

Returns:

Type	Description
DrivePrediction	A new prediction; this one is unchanged.

Raises:

Type	Description
(ValueError, TypeError)	As predict_drives raises for a rating it cannot read.

Source code in src/visualdynamics/core/drive_prediction.py

def with_ratings(self, *, drive_limit: float | Mapping[str, float] | None = None,
                 shakers: ShakerRating | Mapping[str, ShakerRating] | None = None,
                 tables: Mapping[str, str] | None = None,
                 limits: Mapping[str, float] | None = None,
                 duration: float = DURATION,
                 peak_factor: float | None = None) -> DrivePrediction:
    """This prediction judged against other ratings, the drives and
    responses kept: the inverse does not depend on the ratings, so a
    rating moved does not predict again. The ratings are taken
    exactly as `predict_drives` takes them, and one not given here
    is not judged, whatever the prediction was made with.

    Returns
    -------
    DrivePrediction
        A new prediction; this one is unchanged.

    Raises
    ------
    ValueError, TypeError
        As `predict_drives` raises for a rating it cannot read.
    """
    from dataclasses import replace

    out = replace(self, duration=_positive(duration, 'duration'),
                  peak_factor=(None if peak_factor is None else
                               _positive(peak_factor, 'peak factor')),
                  demands=[], notes=list(self.notes[:self.plain_notes]))
    out.demands.extend(_demands(out, out.row_channels, out.drive_channels,
                                drive_limit, shakers, tables, limits))
    return out

drive_asd

drive_asd(drive: int | str = 0) -> ndarray

One drive's predicted autospectrum, by index or label.

Source code in src/visualdynamics/core/drive_prediction.py

def drive_asd(self, drive: int | str = 0) -> np.ndarray:
    """One drive's predicted autospectrum, by index or label."""
    j = self.drives.index(drive) if isinstance(drive, str) else int(drive)
    return np.real(self.drive_csd[:, j, j])

drive_coherence

drive_coherence(first: int | str, second: int | str) -> ndarray

The coherence between two predicted drives at every line: one where the specification makes them copies of one another, zero where it makes them independent.

Source code in src/visualdynamics/core/drive_prediction.py

def drive_coherence(self, first: int | str, second: int | str
                    ) -> np.ndarray:
    """The coherence between two predicted drives at every line:
    one where the specification makes them copies of one another,
    zero where it makes them independent."""
    i = self.drives.index(first) if isinstance(first, str) else int(first)
    j = (self.drives.index(second) if isinstance(second, str)
         else int(second))
    s = self.drive_csd
    with np.errstate(divide='ignore', invalid='ignore'):
        return (np.abs(s[:, i, j]) ** 2
                / np.real(s[:, i, i] * s[:, j, j])).real

missed

missed(control: int | None = None) -> ndarray

The lines where a control channel's prediction misses the target by more than UNMET_DB, by its place among the control channels; with none, where any of them does. A line where the prediction is not a number, or zero against a target, misses.

Source code in src/visualdynamics/core/drive_prediction.py

def missed(self, control: int | None = None) -> np.ndarray:
    """The lines where a control channel's prediction misses the
    target by more than `UNMET_DB`, by its place among the control
    channels; with none, where any of them does. A line where the
    prediction is not a number, or zero against a target, misses."""
    gap = np.abs(self.shortfall_db)
    off = ~(np.where(np.isfinite(gap), gap, np.inf) <= UNMET_DB)
    off = off.any(axis=1) if control is None else off[:, int(control)]
    return self.judged & off

describe

describe(unit_system: UnitSystem | None = None) -> list[str]

The prediction in sentences: the verdict, each demand, the unmet bands and the notes.

Source code in src/visualdynamics/core/drive_prediction.py

def describe(self, unit_system: UnitSystem | None = None) -> list[str]:
    """The prediction in sentences: the verdict, each demand, the
    unmet bands and the notes."""
    lines = []
    first = self.first_out
    if first is None:
        lines.append(f'{self.title}: no ratings given')
    elif first.exceeds:
        lines.append(f'{self.title}: the {first.name} runs out first, '
                     f'at {100.0 * first.share:.0f} % of its rating, '
                     f'most at {first.frequency:.4g} Hz')
    else:
        lines.append(f'{self.title}: within every rating; the '
                     f'{first.name} is nearest, at '
                     f'{100.0 * first.share:.0f} % '
                     f'({self.headroom_db:.1f} dB of headroom), most at '
                     f'{first.frequency:.4g} Hz')
    bands = self.unmet_bands
    if bands:
        # the verdict is read on the status line alone, and drives
        # that fall short of the target ask less of the rig than the
        # test will: the shortfall travels with the verdict
        lines[0] += (f'; short of the target by more than {UNMET_DB:g} dB'
                     f' in {len(bands)} band{"s" * (len(bands) != 1)}, '
                     'where these drives are less than the test asks')
    lines.extend(d.describe(unit_system) for d in self.demands)
    if bands:
        spans = ', '.join(f'{lo:.4g} Hz' if lo == hi else
                          f'{lo:.4g} to {hi:.4g} Hz'
                          for lo, hi in bands[:3])
        more = (f' and {len(bands) - 3} more band'
                f'{"s" * (len(bands) - 3 != 1)}' if len(bands) > 3
                else '')
        lines.append(f'target not met by the inverse within '
                     f'{UNMET_DB:g} dB at {spans}{more}')
    lines.extend(self.notes)
    return lines

Functions:

regularized_inverse

regularized_inverse(plant: Any, *, method: str = 'truncated', tolerance: float = TRUNCATION, weight: float = TIKHONOV, scale: float | None = None) -> tuple[ndarray, ndarray]

An inverse of the plant at every line, truncated or Tikhonov.

plant is (lines, rows, columns), one row per control channel and one column per drive. The inverse is (lines, columns, rows), built from the singular value decomposition at each line: the truncated one keeps the singular values of at least tolerance times the line's largest, and the Tikhonov one filters every singular value sigma to sigma / (sigma^2 + lambda), with lambda weight times scale, the largest squared singular value over every line when scale is not given. A line holding a value that is not a number has an inverse of not-a-number; a line that is zero has an inverse of zero. The module docstring says how the two differ.

Parameters:

Name	Type	Description	Default
plant	array_like	The FRF matrix, (lines, rows, columns), complex.	required
method	('truncated', 'tikhonov')	Which inverse.	'truncated'
tolerance	float	The truncated inverse's relative threshold, between 0 and 1.	`TRUNCATION`
weight	float	The Tikhonov weight, above zero.	`TIKHONOV`
scale	float	The squared singular value the weight is a share of.	None

Returns:

Type	Description
tuple of numpy.ndarray	The inverse, (lines, columns, rows), and the number of directions kept at each line: the singular values kept by the truncation, or those whose square is at least lambda.

Raises:

Type	Description
ValueError	If method is not one of METHODS, or the tolerance or the weight is out of its range.

Source code in src/visualdynamics/core/drive_prediction.py

def regularized_inverse(plant: Any, *, method: str = 'truncated',
                        tolerance: float = TRUNCATION,
                        weight: float = TIKHONOV,
                        scale: float | None = None
                        ) -> tuple[np.ndarray, np.ndarray]:
    """An inverse of the plant at every line, truncated or Tikhonov.

    `plant` is (lines, rows, columns), one row per control channel and
    one column per drive. The inverse is (lines, columns, rows), built
    from the singular value decomposition at each line: the truncated
    one keeps the singular values of at least `tolerance` times the
    line's largest, and the Tikhonov one filters every singular value
    sigma to sigma / (sigma^2 + lambda), with lambda `weight` times
    `scale`, the largest squared singular value over every line when
    `scale` is not given. A line holding a value that is not a number
    has an inverse of not-a-number; a line that is zero has an inverse
    of zero. The module docstring says how the two differ.

    Parameters
    ----------
    plant : array_like
        The FRF matrix, (lines, rows, columns), complex.
    method : {'truncated', 'tikhonov'}, default 'truncated'
        Which inverse.
    tolerance : float, default `TRUNCATION`
        The truncated inverse's relative threshold, between 0 and 1.
    weight : float, default `TIKHONOV`
        The Tikhonov weight, above zero.
    scale : float, optional
        The squared singular value the weight is a share of.

    Returns
    -------
    tuple of numpy.ndarray
        The inverse, (lines, columns, rows), and the number of
        directions kept at each line: the singular values kept by the
        truncation, or those whose square is at least lambda.

    Raises
    ------
    ValueError
        If `method` is not one of `METHODS`, or the tolerance or the
        weight is out of its range.
    """
    _check_inverse(method, tolerance, weight)
    plant = np.asarray(plant, dtype=np.complex128)
    lines, rows, columns = plant.shape
    inverse = np.full((lines, columns, rows), np.nan, dtype=np.complex128)
    kept = np.zeros(lines, dtype=int)
    finite = np.isfinite(plant).all(axis=(1, 2))
    if not finite.any():
        return inverse, kept
    inverse[finite], kept[finite] = inverse_from_svd(
        *np.linalg.svd(plant[finite], full_matrices=False), method=method,
        tolerance=tolerance, weight=weight, scale=scale)
    return inverse, kept

inverse_from_svd

inverse_from_svd(u: Any, sigma: Any, vh: Any, *, method: str = 'truncated', tolerance: float = TRUNCATION, weight: float = TIKHONOV, scale: float | None = None) -> tuple[ndarray, ndarray]

regularized_inverse from the plant's singular value decomposition, (lines, rows, k), (lines, k) and (lines, k, columns) as numpy.linalg.svd returns them with full_matrices=False: the prediction decomposes the plant once and reads its condition number off the same singular values.

Parameters:

Name	Type	Description	Default
u	array_like	The decomposition, every line finite.	required
sigma	array_like	The decomposition, every line finite.	required
vh	array_like	The decomposition, every line finite.	required
method	str	As regularized_inverse takes them.	'truncated'
tolerance	str	As regularized_inverse takes them.	'truncated'
weight	str	As regularized_inverse takes them.	'truncated'
scale	str	As regularized_inverse takes them.	'truncated'

Returns:

Type	Description
tuple of numpy.ndarray	The inverse, (lines, columns, rows), and the directions kept at each line.

Source code in src/visualdynamics/core/drive_prediction.py

def inverse_from_svd(u: Any, sigma: Any, vh: Any, *,
                     method: str = 'truncated', tolerance: float = TRUNCATION,
                     weight: float = TIKHONOV, scale: float | None = None
                     ) -> tuple[np.ndarray, np.ndarray]:
    """`regularized_inverse` from the plant's singular value
    decomposition, (lines, rows, k), (lines, k) and (lines, k, columns)
    as `numpy.linalg.svd` returns them with `full_matrices=False`:
    the prediction decomposes the plant once and reads its condition
    number off the same singular values.

    Parameters
    ----------
    u, sigma, vh : array_like
        The decomposition, every line finite.
    method, tolerance, weight, scale
        As `regularized_inverse` takes them.

    Returns
    -------
    tuple of numpy.ndarray
        The inverse, (lines, columns, rows), and the directions kept
        at each line.
    """
    # The prediction decomposed the plant twice, once for the condition
    # number and once here: two thirds of a 20 x 20 prediction's time
    # (2026-10-28).
    _check_inverse(method, tolerance, weight)
    method, tolerance, weight = str(method), float(tolerance), float(weight)
    u = np.asarray(u, dtype=np.complex128)
    sigma = np.asarray(sigma, dtype=float)
    vh = np.asarray(vh, dtype=np.complex128)
    largest = sigma[:, :1]
    if method == 'truncated':
        keep = (sigma > 0.0) & (sigma >= tolerance * largest)
        with np.errstate(divide='ignore'):
            gain = np.where(keep, 1.0 / np.where(keep, sigma, 1.0), 0.0)
    else:
        if scale is None:
            scale = float(np.max(largest) ** 2)
        lam = weight * float(scale)
        keep = sigma ** 2 >= lam if lam > 0.0 else sigma > 0.0
        denominator = sigma ** 2 + lam
        with np.errstate(divide='ignore', invalid='ignore'):
            gain = np.where(denominator > 0.0, sigma / denominator, 0.0)
    # V diag(gain) U^H, line by line
    inverse = np.einsum('lji,lj,lkj->lik', vh.conj(), gain, u.conj())
    return inverse, keep.sum(axis=1)

through

through(plant: Any, density: Any) -> ndarray

The plant's conjugate times density times the plant's transpose, conj(H) S H^T, at every line: the spectral density matrix of the outputs of a linear system whose inputs have the spectral density matrix density, each cross term stored as compute_cpsds stores it, conj(X_a) X_b. With an inverse for plant, the inputs that produce a target. Both (lines, ...) arrays; the result is (lines, rows, rows).

Source code in src/visualdynamics/core/drive_prediction.py

def through(plant: Any, density: Any) -> np.ndarray:
    """The plant's conjugate times `density` times the plant's
    transpose, `conj(H) S H^T`, at every line: the spectral density
    matrix of the outputs of a linear system whose inputs have the
    spectral density matrix `density`, each cross term stored as
    `compute_cpsds` stores it, conj(X_a) X_b. With an inverse for
    `plant`, the inputs that produce a target. Both (lines, ...)
    arrays; the result is (lines, rows, rows)."""
    plant = np.asarray(plant, dtype=np.complex128)
    density = np.asarray(density, dtype=np.complex128)
    return np.conj(plant) @ density @ np.swapaxes(plant, 1, 2)

predict_drives

predict_drives(frf: Frf, specification: Psd, *, method: str = 'truncated', tolerance: float = TRUNCATION, weight: float = TIKHONOV, drive_limit: float | Mapping[str, float] | None = None, shakers: ShakerRating | Mapping[str, ShakerRating] | None = None, tables: Mapping[str, str] | None = None, limits: Mapping[str, float] | None = None, duration: float = DURATION, peak_factor: float | None = None, progress: Callable[[int, int], None] | None = None) -> DrivePrediction

The drives that would make the specification at its control channels through this FRF, the responses they would make at every channel the FRF holds, and what they ask of the amplifiers and the shakers.

prediction = predict_drives(
    frf, specification, drive_limit=10.0,
    shakers=ShakerRating(force=26700.0, stroke=0.0508,
                         velocity=1.8, moving_mass=60.0))

The FRF is from the drives (volts at the amplifier inputs, or force) to the control channels and to any other channel measured against the same drives; the specification names the control channels, a full cross-spectral matrix when there are several and its autospectrum when there is one. The drives are conj(Z) S_ref Z^T at every line, with Z the regularized inverse method names, and the module docstring says how each is regularized and why the conjugate is on the left, how the peaks are read and what each rating is judged against. The prediction is open loop: what a controller's first drive would be with this FRF as its plant, made once.

Parameters:

Name	Type	Description	Default
frf	Frf	Drives to responses, units declared, every control channel measured against every drive.	required
specification	Psd	The target at the control channels (a Specification, or any spectral density matrix), units declared.	required
method	('truncated', 'tikhonov')	The inverse (regularized_inverse).	'truncated'
tolerance	float	The truncated inverse's share of each line's largest singular value.	`TRUNCATION`
weight	float	The Tikhonov weight, a share of the largest squared singular value in the band.	`TIKHONOV`
drive_limit	float or mapping	The peak drive each drive may reach, in the drive's SI unit (V, or N for a force drive); one number for every drive, or by drive.	None
shakers	ShakerRating or mapping	The shaker behind each drive (core.shaker.ShakerRating); one for every drive, or by drive.	None
tables	mapping	The channel each shaker's table is measured at, by drive. By default the acceleration at the drive's own DOF, or the control channel of a single-input test.	None
limits	mapping	A peak limit for any channel of the FRF, by channel, in its SI unit: a force gauge in N, a current monitor in its volts.	None
duration	float	The run's length in seconds, for the expected largest peaks.	`DURATION`
peak_factor	float	A stated peak over rms for every quantity, in place of the expected largest peak.	None
progress	callable	Told (done, total) as the stages finish (the plant read, its inverse, the drives and responses, the demands); what it raises stops the prediction. The app's progress bar uses it.	None

Returns:

Type	Description
DrivePrediction	The drives' matrix, every channel's predicted autospectrum, the control channels' prediction beside the target, the demands against their ratings, and the notes.

Raises:

Type	Description
TypeError	If frf is not an FRF, specification not a spectral density, or a shaker's rating not a ShakerRating.
ValueError	If either object holds no records or undeclared units, the drives are of two quantities, a control channel is missing from the FRF or held there as another quantity, a control channel is not measured against every drive, the specification lacks an autospectrum or (with several control channels) a cross term, the two share no line, the plant is zero or not a number at every line shared, a setting is out of range, or a rating names a drive or channel the FRF does not have or a shaker has no acceleration measured at its table.

Source code in src/visualdynamics/core/drive_prediction.py

def predict_drives(frf: Frf, specification: Psd, *,
                   method: str = 'truncated', tolerance: float = TRUNCATION,
                   weight: float = TIKHONOV,
                   drive_limit: float | Mapping[str, float] | None = None,
                   shakers: ShakerRating | Mapping[str, ShakerRating]
                   | None = None,
                   tables: Mapping[str, str] | None = None,
                   limits: Mapping[str, float] | None = None,
                   duration: float = DURATION,
                   peak_factor: float | None = None,
                   progress: Callable[[int, int], None] | None = None
                   ) -> DrivePrediction:
    """The drives that would make the specification at its control
    channels through this FRF, the responses they would make at every
    channel the FRF holds, and what they ask of the amplifiers and the
    shakers.

        prediction = predict_drives(
            frf, specification, drive_limit=10.0,
            shakers=ShakerRating(force=26700.0, stroke=0.0508,
                                 velocity=1.8, moving_mass=60.0))

    The FRF is from the drives (volts at the amplifier inputs, or force)
    to the control channels and to any other channel measured against
    the same drives; the specification names the control channels, a
    full cross-spectral matrix when there are several and its
    autospectrum when there is one. The drives are `conj(Z) S_ref Z^T`
    at every line, with `Z` the regularized inverse `method` names, and
    the module docstring says how each is regularized and why the
    conjugate is on the left, how the peaks are read and what each
    rating is judged against. The prediction is open
    loop: what a controller's first drive would be with this FRF as its
    plant, made once.

    Parameters
    ----------
    frf : Frf
        Drives to responses, units declared, every control channel
        measured against every drive.
    specification : Psd
        The target at the control channels (a `Specification`, or any
        spectral density matrix), units declared.
    method : {'truncated', 'tikhonov'}, default 'truncated'
        The inverse (`regularized_inverse`).
    tolerance : float, default `TRUNCATION`
        The truncated inverse's share of each line's largest singular
        value.
    weight : float, default `TIKHONOV`
        The Tikhonov weight, a share of the largest squared singular
        value in the band.
    drive_limit : float or mapping, optional
        The peak drive each drive may reach, in the drive's SI unit
        (V, or N for a force drive); one number for every drive, or by
        drive.
    shakers : ShakerRating or mapping, optional
        The shaker behind each drive (`core.shaker.ShakerRating`); one
        for every drive, or by drive.
    tables : mapping, optional
        The channel each shaker's table is measured at, by drive. By
        default the acceleration at the drive's own DOF, or the control
        channel of a single-input test.
    limits : mapping, optional
        A peak limit for any channel of the FRF, by channel, in its SI
        unit: a force gauge in N, a current monitor in its volts.
    duration : float, default `DURATION`
        The run's length in seconds, for the expected largest peaks.
    peak_factor : float, optional
        A stated peak over rms for every quantity, in place of the
        expected largest peak.
    progress : callable, optional
        Told ``(done, total)`` as the stages finish (the plant read, its
        inverse, the drives and responses, the demands); what it raises
        stops the prediction. The app's progress bar uses it.

    Returns
    -------
    DrivePrediction
        The drives' matrix, every channel's predicted autospectrum, the
        control channels' prediction beside the target, the demands
        against their ratings, and the notes.

    Raises
    ------
    TypeError
        If `frf` is not an FRF, `specification` not a spectral density,
        or a shaker's rating not a `ShakerRating`.
    ValueError
        If either object holds no records or undeclared units, the
        drives are of two quantities, a control channel is missing from
        the FRF or held there as another quantity, a control channel is
        not measured against every drive, the specification lacks an
        autospectrum or (with several control channels) a cross term,
        the two share no line, the plant is zero or not a number at
        every line shared, a setting is out of range, or a rating names
        a drive or channel the FRF does not have or a shaker has no
        acceleration measured at its table.
    """
    from .data import Frf, Psd

    if not isinstance(frf, Frf):
        raise TypeError(f'a prediction reads an FRF, not a '
                         f'{type(frf).__name__}')
    if not isinstance(specification, Psd):
        raise TypeError(f'a prediction reads a specification, not a '
                         f'{type(specification).__name__}')
    if frf.num_records == 0 or specification.num_records == 0:
        raise ValueError('nothing to predict: '
                         + ('the FRF' if frf.num_records == 0
                            else 'the specification') + ' holds no records')
    duration = _positive(duration, 'duration')
    if peak_factor is not None:
        peak_factor = _positive(peak_factor, 'peak factor')
    # the settings are checked before any arithmetic, so a refusal never
    # waits on the inverse
    regularized_inverse(np.zeros((0, 1, 1)), method=method,
                        tolerance=tolerance, weight=weight)

    ticker = Ticker(progress)
    ticker.add(4)
    rows, columns, matrix, filled = _plant(frf)
    row_labels = channel_labels(rows)
    drive_labels = channel_labels(columns)
    controls = _control_channels(specification, rows)
    control_index = [rows.index(key) for key in controls]
    for r in control_index:
        missing = [drive_labels[c] for c in range(len(columns))
                   if not filled[r, c]]
        if missing:
            raise ValueError(
                f'the control channel {row_labels[r]} is not measured '
                f'against {", ".join(missing)}; every control channel '
                'needs its FRF from every drive')
    notes: list[str] = []
    complete = filled.all(axis=1)
    if not complete.all():
        dropped = [row_labels[r] for r in np.flatnonzero(~complete)]
        notes.append(f'{", ".join(dropped)} not predicted: not measured '
                     'against every drive')

    f = np.asarray(frf.abscissa, dtype=float)
    target = _read_target(specification, controls, f)
    plant = matrix[:, control_index, :]
    finite = np.isfinite(plant).all(axis=(1, 2))
    written = np.isfinite(np.real(np.einsum('lii->li', target))).all(axis=1)
    judged = (f > 0.0) & written & finite
    if not (written & (f > 0.0)).any():
        raise ValueError('the specification and the FRF share no line '
                         'where both are written')
    nan_lines = int((written & (f > 0.0) & ~finite).sum())
    where = _not_a_number(frf)
    if not judged.any():
        raise ValueError('the FRF is not a number at every line the '
                         f'specification is written over, first at {where}')
    if nan_lines:
        notes.append(f'{nan_lines} line{"s" * (nan_lines != 1)} where the '
                     f'FRF is not a number not predicted, first at {where}')
    ticker.tick()
    ticker.busy()
    # one decomposition: the condition number, the dead lines and the
    # inverse all read from it
    u, gains, vh = np.linalg.svd(plant[judged], full_matrices=False)
    if not gains.max() > 0.0:
        raise ValueError('the FRF is zero at every line the specification '
                         'is written over: the drives reach nothing')
    # a drive that reaches no control channel, or a control channel no
    # drive reaches, is a dead channel: named, since the inverse sets
    # its direction aside without a word otherwise
    reached = np.abs(plant[judged]).max(axis=0)
    for c in range(len(columns)):
        if not reached[:, c].any():
            notes.append(f'the drive {drive_labels[c]} reaches no control '
                         'channel at any line: it is predicted to do '
                         'nothing')
    for k, r in enumerate(control_index):
        if not reached[k].any():
            notes.append(f'no drive reaches the control channel '
                         f'{row_labels[r]} at any line: its target cannot '
                         'be met')
    with np.errstate(divide='ignore', invalid='ignore'):
        condition = np.full(len(f), np.nan)
        condition[judged] = gains[:, 0] / gains[:, -1]
    dead = int((gains[:, 0] == 0.0).sum())
    if dead:
        notes.append(f'the plant is zero at {dead} line{"s" * (dead != 1)}: '
                     'nothing can be driven there')

    inverse = np.full((len(f), len(columns), len(controls)), np.nan,
                      dtype=np.complex128)
    kept = np.zeros(len(f), dtype=int)
    inverse[judged], kept[judged] = inverse_from_svd(
        u, gains, vh, method=method, tolerance=tolerance, weight=weight)
    ticker.tick()
    drive_csd = np.full((len(f), len(columns), len(columns)), np.nan,
                        dtype=np.complex128)
    drive_csd[judged] = through(inverse[judged], target[judged])
    responses = np.full((len(f), len(rows), len(rows)), np.nan,
                        dtype=np.complex128)
    responses[judged] = through(matrix[judged], drive_csd[judged])
    asd = np.real(np.einsum('lii->li', responses))
    asd[:, ~complete] = np.nan
    control = responses[:, control_index][:, :, control_index]
    target = np.where(judged[:, None, None], target, np.nan)
    ticker.tick()

    if len(columns) > len(controls):
        notes.append('more drives than control channels: of the drives '
                     'that meet the target, the inverse takes the one of '
                     'least power')
    elif len(columns) < len(controls):
        notes.append('fewer drives than control channels: the target is '
                     'met in the least-squares sense, and the cross terms '
                     'the drives cannot make are not')
    short = judged & (kept < min(len(columns), len(controls)))
    if short.any():
        count = int(short.sum())
        notes.append(f'the inverse set aside a direction at {count} '
                     f'line{"s" * (count != 1)} where the plant is '
                     'ill-conditioned')

    prediction = DrivePrediction(
        frequency=f, drives=drive_labels, drive_quantity=columns[0][1],
        rows=row_labels, row_quantities=[q for _d, q in rows],
        controls=control_index, drive_csd=drive_csd, asd=asd,
        target=target, control=control, judged=judged, kept=kept,
        condition=condition, method=str(method), tolerance=float(tolerance),
        weight=float(weight), duration=duration, peak_factor=peak_factor,
        notes=notes, row_channels=list(rows), drive_channels=list(columns),
        plain_notes=len(notes))
    prediction = prediction.with_ratings(
        drive_limit=drive_limit, shakers=shakers, tables=tables,
        limits=limits, duration=duration, peak_factor=peak_factor)
    ticker.tick()
    return prediction