"""Sensitivity functions of a feedback loop, built from its plant, controller and filter."""
from __future__ import annotations
import numpy as np
from minilink.analysis.frequency import transfer_function
from minilink.analysis.linearization import output_selectors
from minilink.core.feedback import free_inputs
###############################################################################
# Public API
###############################################################################
[docs]
def sensitivity(
*,
plant,
controller,
filter=None,
x_bar=None,
u_bar=None,
t=0.0,
params=None,
method: str = "auto",
eps: float = 1e-6,
):
"""Return the sensitivity function ``S`` of the loop as a ``TransferFunction``.
``S`` maps an output disturbance to ``y`` and the reference ``r`` to the
error ``e``. The loop is the plant ``H``, the controller ``C`` and the
optional filter ``F`` on the measurement, each a SISO ``System``
linearized first.
Parameters
----------
plant : System
The plant ``H``, from its one input besides ``d`` and ``n`` to its output.
controller : System
The controller ``C``, from the error ``e`` to the command ``u``. A
controller in the reference form (inputs ``r`` and ``y``) is refused:
pass its channel from ``r``, ``transfer_function(C, wrt="r")``.
filter : System, optional
The filter ``F`` on the measurement; ``F = 1`` when omitted.
x_bar, u_bar, t, params : optional
The plant's operating point, as in
:func:`~minilink.analysis.frequency.bode`; the controller and the
filter are taken at their nominal point.
method : {"auto", "fd", "jax"}, optional
Differentiation backend of the linearizations.
eps : float, optional
Central-difference step.
Returns
-------
TransferFunction
``S(s)``, ready for :func:`~minilink.analysis.frequency.bode`,
``plot_bode`` or a simulation; its ports are labelled ``r`` and ``e``
(``T``: ``r`` to ``y``, ``PS``: ``w`` to ``y``, ``CS``: ``v`` to ``u``).
"""
# The loop's pieces, each a SISO transfer function
H, C, F = loop_pieces(
plant, controller, filter, x_bar, u_bar, t, params, method, eps
)
# L(s) = C(s) H(s) F(s)
L_num = np.polymul(np.polymul(C.numerator, H.numerator), F.numerator)
L_den = np.polymul(np.polymul(C.denominator, H.denominator), F.denominator)
# S = 1 / (1 + L)
S_num = L_den
S_den = np.polyadd(L_den, L_num)
return loop_function(S_num, S_den, "Sensitivity S", "r", "e")
[docs]
def complementary_sensitivity(
*,
plant,
controller,
filter=None,
x_bar=None,
u_bar=None,
t=0.0,
params=None,
method: str = "auto",
eps: float = 1e-6,
):
"""Return the complementary sensitivity ``T``, the map from ``r`` to ``y``.
Same arguments as :func:`sensitivity`; without a filter, ``S + T = 1``.
"""
# The loop's pieces, each a SISO transfer function
H, C, F = loop_pieces(
plant, controller, filter, x_bar, u_bar, t, params, method, eps
)
# L(s) = C(s) H(s) F(s)
L_num = np.polymul(np.polymul(C.numerator, H.numerator), F.numerator)
L_den = np.polymul(np.polymul(C.denominator, H.denominator), F.denominator)
# T = C H / (1 + L)
T_num = np.polymul(np.polymul(C.numerator, H.numerator), F.denominator)
T_den = np.polyadd(L_den, L_num)
return loop_function(T_num, T_den, "Complementary sensitivity T", "r", "y")
[docs]
def load_sensitivity(
*,
plant,
controller,
filter=None,
x_bar=None,
u_bar=None,
t=0.0,
params=None,
method: str = "auto",
eps: float = 1e-6,
):
"""Return the load sensitivity ``PS``, the map from a disturbance ``d`` at the plant input to ``y``.
Same arguments as :func:`sensitivity`.
"""
# The loop's pieces, each a SISO transfer function
H, C, F = loop_pieces(
plant, controller, filter, x_bar, u_bar, t, params, method, eps
)
# L(s) = C(s) H(s) F(s)
L_num = np.polymul(np.polymul(C.numerator, H.numerator), F.numerator)
L_den = np.polymul(np.polymul(C.denominator, H.denominator), F.denominator)
# PS = H / (1 + L)
PS_num = np.polymul(np.polymul(H.numerator, C.denominator), F.denominator)
PS_den = np.polyadd(L_den, L_num)
return loop_function(PS_num, PS_den, "Load sensitivity PS", "d", "y")
[docs]
def noise_sensitivity(
*,
plant,
controller,
filter=None,
x_bar=None,
u_bar=None,
t=0.0,
params=None,
method: str = "auto",
eps: float = 1e-6,
):
"""Return the noise sensitivity ``CS``, the map from measurement noise ``n`` to the command ``u``.
Same arguments as :func:`sensitivity`. In the loop, ``n`` reaches ``u``
as ``-CS``: the sign of the negative feedback is left out.
"""
# The loop's pieces, each a SISO transfer function
H, C, F = loop_pieces(
plant, controller, filter, x_bar, u_bar, t, params, method, eps
)
# L(s) = C(s) H(s) F(s)
L_num = np.polymul(np.polymul(C.numerator, H.numerator), F.numerator)
L_den = np.polymul(np.polymul(C.denominator, H.denominator), F.denominator)
# CS = C F / (1 + L)
CS_num = np.polymul(np.polymul(C.numerator, F.numerator), H.denominator)
CS_den = np.polyadd(L_den, L_num)
return loop_function(CS_num, CS_den, "Noise sensitivity CS", "n", "u")
###############################################################################
# Internal machinery
###############################################################################
def loop_pieces(plant, controller, filter, x_bar, u_bar, t, params, method, eps):
"""``H``, ``C`` and ``F`` as SISO transfer functions; ``F = 1`` without a filter."""
# Lazy import: the block and its graphics load on first use, not with the analysis
from minilink.blocks.transfer_function import TransferFunction
# The plant at its operating point; the controller and filter at their nominal one
H = siso_piece(plant, "plant", x_bar, u_bar, t, params, method, eps)
C = siso_piece(controller, "controller", None, None, 0.0, None, method, eps)
if filter is None:
F = TransferFunction([1.0], [1.0], name="F")
else:
F = siso_piece(filter, "filter", None, None, 0.0, None, method, eps)
return H, C, F
def siso_piece(sys, role, x_bar, u_bar, t, params, method, eps):
"""The transfer function of one loop piece; a piece with several channels is refused."""
# p outputs (the state, without an output port) and m inputs (every free input:
# all but the disturbance d and the noise n)
outputs = output_selectors(sys, None)
p = sys.n if outputs is None else sum(sys.outputs[name].dim for name, _ in outputs)
m = sum(sys.inputs[port_id].dim for port_id in free_inputs(sys.inputs))
require_siso(sys, role, m, p)
return transfer_function(
sys, x_bar, u_bar, t, params, method=method, eps=eps, minimal=True
)
def require_siso(sys, role, m, p):
"""Refuse a loop piece that is not single-input, single-output.
The fix it names reads the channel from the reference when the piece has one.
"""
if m != 1 or p != 1:
channel = 'wrt="r"' if "r" in sys.inputs else "of=..., wrt=..."
raise ValueError(
f"The {role} {sys.name!r} must be SISO, got {m} input(s) and {p} "
f"output(s). Select one channel first: "
f"{role}=transfer_function({role}, {channel})."
)
def loop_function(numerator, denominator, name, wrt, of):
"""The closed-loop function as a ``TransferFunction`` block, its ports labelled by the loop signals."""
# Lazy import: the block and its graphics load on first use, not with the analysis
from minilink.blocks.transfer_function import TransferFunction
G = TransferFunction(
np.atleast_1d(numerator), np.atleast_1d(denominator), name=name
)
G.inputs["u"].labels = [wrt]
G.outputs["y"].labels = [of]
return G