diff --git a/CITATION.cff b/CITATION.cff index 71a5d66..9d40ccc 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -14,5 +14,5 @@ authors: orcid: "https://orcid.org/0000-0002-1987-9268" affiliation: "McGill University" title: "decargroup/dkpy" -version: v0.1.11 +version: v0.1.12 url: "https://github.com/decargroup/dkpy" diff --git a/README.rst b/README.rst index 2ec16e4..677ff7b 100644 --- a/README.rst +++ b/README.rst @@ -141,6 +141,6 @@ If you use this software in your research, please cite it as below or see url={https://github.com/decargroup/dkpy}, publisher={Zenodo}, author={Steven Dahdah and Timothy Everett Adams and James Richard Forbes}, - version = {{v0.1.11}}, + version = {{v0.1.12}}, year={2025}, } diff --git a/conftest.py b/conftest.py index 5db9dbe..29cb351 100644 --- a/conftest.py +++ b/conftest.py @@ -32,3 +32,14 @@ def add_example_skogestad2006_p325(doctest_namespace): eg["n_u"], eg["K"], ) + + +@pytest.fixture(autouse=True) +def add_example_multimodel_uncertainty(doctest_namespace): + """Generate uncertain models using parameter variation.""" + eg = dkpy.example_multimodel_uncertainty() + doctest_namespace["example_multimodel_uncertainty"] = ( + eg["complex_response_nominal"], + eg["complex_response_offnominal_list"], + eg["omega"], + ) diff --git a/doc/_static/example_7/7_magnitude_nom_offnom.png b/doc/_static/example_7/7_magnitude_nom_offnom.png new file mode 100644 index 0000000..7ac542f Binary files /dev/null and b/doc/_static/example_7/7_magnitude_nom_offnom.png differ diff --git a/doc/_static/example_7/7_phase_nom_offnom.png b/doc/_static/example_7/7_phase_nom_offnom.png new file mode 100644 index 0000000..c0eb40b Binary files /dev/null and b/doc/_static/example_7/7_phase_nom_offnom.png differ diff --git a/doc/_static/example_7/7_sval_max_residual.png b/doc/_static/example_7/7_sval_max_residual.png new file mode 100644 index 0000000..d030c9a Binary files /dev/null and b/doc/_static/example_7/7_sval_max_residual.png differ diff --git a/doc/_static/example_7/7_sval_nom_offnom.png b/doc/_static/example_7/7_sval_nom_offnom.png new file mode 100644 index 0000000..62c76cc Binary files /dev/null and b/doc/_static/example_7/7_sval_nom_offnom.png differ diff --git a/doc/_static/example_7/7_sval_residual_A.png b/doc/_static/example_7/7_sval_residual_A.png new file mode 100644 index 0000000..c7c152f Binary files /dev/null and b/doc/_static/example_7/7_sval_residual_A.png differ diff --git a/doc/_static/example_7/7_sval_residual_I.png b/doc/_static/example_7/7_sval_residual_I.png new file mode 100644 index 0000000..7a1fd1a Binary files /dev/null and b/doc/_static/example_7/7_sval_residual_I.png differ diff --git a/doc/_static/example_7/7_sval_residual_O.png b/doc/_static/example_7/7_sval_residual_O.png new file mode 100644 index 0000000..b565735 Binary files /dev/null and b/doc/_static/example_7/7_sval_residual_O.png differ diff --git a/doc/_static/example_7/7_sval_residual_iA.png b/doc/_static/example_7/7_sval_residual_iA.png new file mode 100644 index 0000000..5e49540 Binary files /dev/null and b/doc/_static/example_7/7_sval_residual_iA.png differ diff --git a/doc/_static/example_7/7_sval_residual_iI.png b/doc/_static/example_7/7_sval_residual_iI.png new file mode 100644 index 0000000..906962f Binary files /dev/null and b/doc/_static/example_7/7_sval_residual_iI.png differ diff --git a/doc/_static/example_7/7_sval_residual_iO.png b/doc/_static/example_7/7_sval_residual_iO.png new file mode 100644 index 0000000..6cc0d1e Binary files /dev/null and b/doc/_static/example_7/7_sval_residual_iO.png differ diff --git a/doc/_static/example_7/7_uncertainty_weight.png b/doc/_static/example_7/7_uncertainty_weight.png new file mode 100644 index 0000000..dccd187 Binary files /dev/null and b/doc/_static/example_7/7_uncertainty_weight.png differ diff --git a/doc/conf.py b/doc/conf.py index bfc9073..df5b6ff 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -9,8 +9,8 @@ project = "dkpy" copyright = "2024, Steven Dahdah, Timothy Everett Adams and James Richard Forbes" author = "Steven Dahdah, Timothy Everett Adams and James Richard Forbes" -version = "0.1.11" -release = "0.1.11" +version = "0.1.12" +release = "0.1.12" # -- General configuration --------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration diff --git a/doc/examples.rst b/doc/examples.rst index 7ecd069..e3ef9bc 100644 --- a/doc/examples.rst +++ b/doc/examples.rst @@ -177,3 +177,74 @@ particular, the response of the system states and actuator inputs are shown. .. image:: _static/example_6/6_plot_states.png .. image:: _static/example_6/6_plot_inputs.png + +Multi-Model Uncertainty Characterization +---------------------------------------- + +In this example, an unstructured uncertainty set is characterized from a set +of nominal and off-nominal frequency responses. The off-nominal models are +generated from a nominal model by randomly perturbing the parameters within a +known bound. + +.. literalinclude:: ../examples/7_uncertainty_characterization.py + :language: python + +`dkpy` provides plotting functionality for the nominal and off-nominal systems. +The frequency response of the magnitude, phase, and singular values are +as follows. These plots show the variation in the frequency response caused +by the variation in parameters. + +.. image:: _static/example_7/7_magnitude_nom_offnom.png + +.. image:: _static/example_7/7_phase_nom_offnom.png + +.. image:: _static/example_7/7_sval_nom_offnom.png + +Next, the uncertainty residual can be computed for each off-nominal system and +uncertainty model. `dkpy` implements six unstructured uncertainty models: + +* Additive uncertainty ("A"); + .. image:: _static/example_7/7_sval_residual_A.png +* Multiplicative input uncertainty ("I"); + .. image:: _static/example_7/7_sval_residual_I.png +* Multiplicative output uncertainty ("O"); + .. image:: _static/example_7/7_sval_residual_O.png +* Inverse additive uncertainty ("iA"); + .. image:: _static/example_7/7_sval_residual_iA.png +* Inverse multiplicative input uncertainty ("iI"); + .. image:: _static/example_7/7_sval_residual_iI.png +* Inverse multiplicative output uncertainty ("iO"). + .. image:: _static/example_7/7_sval_residual_iO.png + +See Chapter 8.2.3 of [SP06]_ for more information on these unstructured +uncertainty models. The maximum singular value responses of all uncertainty +models are shown here. + +.. image:: _static/example_7/7_sval_max_residual.png + +As a rule of thumb, the uncertainty model that is selected should have the +smallest maximum singular value over the control bandwidth and high uncertainty +at large frequencies to account for unmodeled dynamics. In this case, the +inverse additive uncertainty ("iA") model seems to be the best option, which will be +used moving forward. + +The uncertainty set is parametrized as `E = WL Δ WR`. `E` is the uncertainty +residual, `Δ` is the normalized perturbation, and `WL`, `WR` are the dynamic +weights used to parametrize the uncertainty set. The frequency response of the +optimal weights with minimal magnitude at each frequency can be solved for from +the uncertainty residuals for different assumptions on their structure. It is +assumed that both weights are diagonal moving forward. Then, an overbounding +stable and minimum phase linear time-invariant (LTI) system can be fit to the +frequency response of the weights to obtain a LTI description of the uncertainty +set. The optimal weight frequency responses and fitted weights are shown below. + +.. image:: _static/example_7/7_uncertainty_weight.png + + + + + + + + + diff --git a/examples/7_uncertainty_characterization.py b/examples/7_uncertainty_characterization.py new file mode 100644 index 0000000..0154fa1 --- /dev/null +++ b/examples/7_uncertainty_characterization.py @@ -0,0 +1,92 @@ +"""Multi-model uncertainty characterization from frequency response data.""" + +from matplotlib import pyplot as plt + +import dkpy + + +def example_uncertainty_characterization(): + """Multi-model uncertainty characterization from frequency response data.""" + + # Generate example data + eg = dkpy.example_multimodel_uncertainty() + response_nom = eg["complex_response_nominal"] + response_offnom_list = eg["complex_response_offnominal_list"] + omega = eg["omega"] + + # Plot: Magnitude response of nominal and off-nominal systems + fig, _ = dkpy.plot_magnitude_response_uncertain_model_set( + response_nom, + response_offnom_list, + omega, + ) + + # Plot: Phase response of nominal and off-nominal systems + fig, _ = dkpy.plot_phase_response_uncertain_model_set( + response_nom, + response_offnom_list, + omega, + ) + + # Plot: Singular value response of nominal and off-nominal systems + fig, ax = dkpy.plot_singular_value_response_uncertain_model_set( + response_nom, + response_offnom_list, + omega, + ) + + # Uncertainty models + uncertainty_models = { + "additive", + "multiplicative_input", + "multiplicative_output", + "inverse_additive", + "inverse_multiplicative_input", + "inverse_multiplicative_output", + } + response_residuals_dict = dkpy.compute_uncertainty_residual_response( + response_nom, + response_offnom_list, + uncertainty_models, + ) + + # Plot: Singular value response of uncerainty residuals + figure_dict = dkpy.plot_singular_value_response_residual( + response_residuals_dict, omega + ) + + # Plot: Comparison of singular value response of uncerainty residuals for each + # uncertainty model + fig, _ = dkpy.plot_singular_value_response_residual_comparison( + response_residuals_dict, omega + ) + + # Compute uncertainty weight frequency response + response_weight_left, response_weight_right = ( + dkpy.compute_optimal_uncertainty_weight_response( + response_residuals_dict["inverse_additive"], "diagonal", "diagonal" + ) + ) + + # Fit overbounding stable and minimum-phase uncertainty weight system + weight_left = dkpy.fit_overbounding_uncertainty_weight( + response_weight_left, omega, [4, 5] + ) + weight_right = dkpy.fit_overbounding_uncertainty_weight( + response_weight_right, omega, [3, 5] + ) + + # Plot: Magnitude response of uncertainty weight frequency response and overbounding + # fit + fig, _ = dkpy.plot_magnitude_response_uncertainty_weight( + response_weight_left, + response_weight_right, + omega, + weight_left, + weight_right, + ) + plt.show() + + +if __name__ == "__main__": + example_uncertainty_characterization() diff --git a/pyproject.toml b/pyproject.toml index e107e49..44b8b2b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "hatchling.build" [project] name = "dkpy" -version = "0.1.11" +version = "0.1.12" dependencies = [ "numpy>=1.21.0", "scipy>=1.7.0", diff --git a/src/dkpy/__init__.py b/src/dkpy/__init__.py index de6b21e..4a48a14 100644 --- a/src/dkpy/__init__.py +++ b/src/dkpy/__init__.py @@ -5,4 +5,5 @@ from .d_scale_fit import * from .structured_singular_value import * from .uncertainty_structure import * +from .uncertainty_characterization import * from .utilities import * diff --git a/src/dkpy/uncertainty_characterization.py b/src/dkpy/uncertainty_characterization.py new file mode 100644 index 0000000..b5b2a9d --- /dev/null +++ b/src/dkpy/uncertainty_characterization.py @@ -0,0 +1,1764 @@ +"""Uncertainty characterization utilties.""" + +__all__ = [ + "compute_uncertainty_residual_response", + "compute_optimal_uncertainty_weight_response", + "fit_overbounding_uncertainty_weight", + "plot_magnitude_response_uncertain_model_set", + "plot_phase_response_uncertain_model_set", + "plot_singular_value_response_uncertain_model_set", + "plot_singular_value_response_residual", + "plot_singular_value_response_residual_comparison", + "plot_magnitude_response_uncertainty_weight", +] + +import warnings + +import control +import numpy as np +import cvxpy +import scipy +from matplotlib import pyplot as plt + +from typing import List, Optional, Union, Tuple, Dict, Callable, Set, Any +from matplotlib.figure import Figure +from matplotlib.axes import Axes + +from . import utilities + + +def compute_uncertainty_residual_response( + complex_response_nom: Union[np.ndarray, control.FrequencyResponseData], + complex_response_offnom_list: Union[np.ndarray, control.FrequencyResponseList], + uncertainty_model: Union[str, List[str], Set[str]] = { + "additive", + "multiplicative_input", + "multiplicative_output", + "inverse_additive", + "inverse_multiplicative_input", + "inverse_multiplicative_output", + }, + tol_residual_existence: float = 1e-12, +) -> Dict[str, np.ndarray]: + """Compute the residual response of unstructured uncertainty models. + + Parameters + ---------- + complex_response_nom : Union[np.ndarray, control.FrequencyResponseData] + Frequency response of the nominal system. + complex_response_offnom_list : Union[np.ndarray, control.FrequencyResponseList] + Frequency response of the off-nominal system. + uncertainty_model : Union[str, List[str], Set[str]] + Uncertainty model identifiers to compute the residual response. + tol_residual_existence : float + Tolerance for the existence of an uncertainty residual. + + Returns + ------- + Dict[str, np.ndarray] + Frequency response of the uncertainty residuals. + + Raises + ------ + ValueError + An invalid uncertainty model identifier was found in `uncertainty_model`. + + Examples + -------- + Compute the residuals for the six unstructured uncertainty models from the nominal + and off-nominal frequency responses. + + >>> complex_response_nom, complex_response_offnom_list, omega = ( + ... example_multimodel_uncertainty + ... ) + >>> uncertainty_models = { + ... "additive", + ... "multiplicative_input", + ... "multiplicative_output", + ... "inverse_additive", + ... "inverse_multiplicative_input", + ... "inverse_multiplicative_output", + ... } + >>> complex_response_residuals_dict = dkpy.compute_uncertainty_residual_response( + ... complex_response_nom, + ... complex_response_offnom_list, + ... uncertainty_models, + ... ) + """ + + # Convert frequency response data to expected type + complex_response_nom = _convert_frequency_response_data_to_array( + complex_response_nom + ) + complex_response_offnom_list = _convert_frequency_response_list_to_array( + complex_response_offnom_list + ) + + # Uncertainty residual response dictionary + complex_response_residual_dict = {} + + # Check uncertainty model identifiers + uncertainty_model = set(uncertainty_model) + valid_uncertainty_model = { + "additive", + "multiplicative_input", + "multiplicative_output", + "inverse_additive", + "inverse_multiplicative_input", + "inverse_multiplicative_output", + } + if not uncertainty_model.issubset(valid_uncertainty_model): + raise ValueError( + "The uncertainty model identifiers provided in `uncertainty_model` do not " + "all correspond to valid uncertainty models. In particular, " + f"{uncertainty_model.difference(valid_uncertainty_model)} are not valid " + "uncertainty model identifiers. The identifiers are: " + '"additive": Additive uncertainty, ' + '"multiplicative_input": Multiplicative input uncertainty, ' + '"multiplicative_output": Multiplicative output uncertainty, ' + '"inverse_additive": Inverse additive uncertainty, ' + '"inverse_multiplicative_input": Inverse multiplicative input uncertainty, ' + '"inverse_multiplicative_output": Inverse multiplicative output uncertainty.' + ) + + # Additive uncertainty residual response + if "additive" in uncertainty_model: + complex_response_residual_list = _compute_uncertainty_residual_response( + complex_response_nom, + complex_response_offnom_list, + _compute_uncertainty_residual_additive_freq, + tol_residual_existence, + ) + complex_response_residual_dict["additive"] = complex_response_residual_list + + # Multiplicative input uncertainty residual response + if "multiplicative_input" in uncertainty_model: + complex_response_residual_list = _compute_uncertainty_residual_response( + complex_response_nom, + complex_response_offnom_list, + _compute_uncertainty_residual_multiplicative_input_freq, + tol_residual_existence, + ) + complex_response_residual_dict["multiplicative_input"] = ( + complex_response_residual_list + ) + + # Multiplicative output uncertainty residual response + if "multiplicative_output" in uncertainty_model: + complex_response_residual_list = _compute_uncertainty_residual_response( + complex_response_nom, + complex_response_offnom_list, + _compute_uncertainty_residual_multiplicative_output_freq, + tol_residual_existence, + ) + complex_response_residual_dict["multiplicative_output"] = ( + complex_response_residual_list + ) + + # Inverse additive uncertainty residual response + if "inverse_additive" in uncertainty_model: + complex_response_residual_list = _compute_uncertainty_residual_response( + complex_response_nom, + complex_response_offnom_list, + _compute_uncertainty_residual_inverse_additive_freq, + tol_residual_existence, + ) + complex_response_residual_dict["inverse_additive"] = ( + complex_response_residual_list + ) + + # Inverse multiplicative input uncertainty residual response + if "inverse_multiplicative_input" in uncertainty_model: + complex_response_residual_list = _compute_uncertainty_residual_response( + complex_response_nom, + complex_response_offnom_list, + _compute_uncertainty_residual_inverse_multiplicative_input_freq, + tol_residual_existence, + ) + complex_response_residual_dict["inverse_multiplicative_input"] = ( + complex_response_residual_list + ) + + # Inverse multiplicative output uncertainty residual response + if "inverse_multiplicative_output" in uncertainty_model: + complex_response_residual_list = _compute_uncertainty_residual_response( + complex_response_nom, + complex_response_offnom_list, + _compute_uncertainty_residual_inverse_multiplicative_output_freq, + tol_residual_existence, + ) + complex_response_residual_dict["inverse_multiplicative_output"] = ( + complex_response_residual_list + ) + + return complex_response_residual_dict + + +def _compute_uncertainty_residual_response( + complex_response_nom: np.ndarray, + complex_response_offnom_list: np.ndarray, + compute_uncertainty_residual_freq: Callable, + tol_residual_existence: float, +) -> np.ndarray: + """Compute the uncertainty residual response for a given uncertainty model. + + Parameters + ---------- + complex_response_nom : np.ndarray + Frequency response of the nominal system. + complex_response_offnom_list : np.ndarray + Frequency response of the off-nominal system. + compute_uncertainty_residual_freq : Callable, + Uncertainty residual computation function at a given frequency. + tol_residual_existence : float + Tolerance for the existence of an uncertainty residual. + + Returns + ------- + np.ndarray + Uncertainty residual response of the given uncertainty model. + """ + + # Frequency response parameters + num_offnom = complex_response_offnom_list.shape[0] + num_frequency = complex_response_offnom_list.shape[1] + + complex_response_residual_list = [] + for idx_offnom in range(num_offnom): + complex_response_residual = [] + complex_response_offnom = complex_response_offnom_list[idx_offnom, :, :, :] + for idx_freq in range(num_frequency): + complex_response_nom_freq = complex_response_nom[idx_freq, :, :] + complex_response_offnom_freq = complex_response_offnom[idx_freq, :, :] + complex_response_residual_freq = compute_uncertainty_residual_freq( + complex_response_nom_freq, + complex_response_offnom_freq, + tol_residual_existence, + ) + complex_response_residual.append(complex_response_residual_freq) + complex_response_residual = np.array(complex_response_residual, dtype=complex) + complex_response_residual_list.append(complex_response_residual) + complex_response_residual_list = np.array(complex_response_residual_list) + + return complex_response_residual_list + + +def _compute_uncertainty_residual_additive_freq( + complex_response_nom_freq: np.ndarray, + complex_response_offnom_freq: np.ndarray, + tol_residual_existence: Optional[float] = None, +) -> np.ndarray: + """Compute the additive uncertainty residual at a frequency. + + Parameters + ---------- + complex_response_nom_freq : np.ndarray + Nominal frequency response matrix evaluated at a given frequency. + complex_response_offnom_freq : np.ndarray + Frequency response matrix of a single off-nominal system evaluated at a given + frequency. + tol_residual_existence : float + Tolerance for the existence of a multiplicative input uncertainty residual. + + Returns + ------- + np.ndarray + Additive uncertainty residual at a given frequency. + """ + complex_response_residual_freq = ( + complex_response_offnom_freq - complex_response_nom_freq + ) + + return complex_response_residual_freq + + +def _compute_uncertainty_residual_multiplicative_input_freq( + complex_response_nom_freq: np.ndarray, + complex_response_offnom_freq: np.ndarray, + tol_residual_existence: float = 1e-8, +) -> np.ndarray: + """Compute the multiplicative input uncertainty residual at a frequency. + + Parameters + ---------- + complex_response_nom_freq : np.ndarray + Nominal frequency response matrix evaluated at a given frequency. + complex_response_offnom_freq : np.ndarray + Frequency response matrix of a single off-nominal system evaluated at a given + frequency. + tol_residual_existence : float + Tolerance for the existence of a multiplicative input uncertainty residual. + + Returns + ------- + np.ndarray + Multiplicative input uncertainty residual at a given frequency. + + Raises + ------ + ValueError + An input multiplicative uncertainty residual does not exist at the given + frequency as the linear system used to compute the residual does not have a + solution. + """ + + num_inputs = complex_response_nom_freq.shape[1] + num_outputs = complex_response_nom_freq.shape[0] + if num_inputs < num_outputs: + warnings.warn( + "Multipliative input uncertainty models cannot include all possible " + "off-nominal systems when the number of inputs is less than the number of " + "outputs. This may result in an error if an uncertainty residual cannot " + "be found for a given nominal and off-nominal frequency response matrix." + ) + + A = complex_response_nom_freq + B = complex_response_offnom_freq - complex_response_nom_freq + X, residues_lstsq, _, _ = scipy.linalg.lstsq(A, B) + complex_response_residual_freq = X + + if num_inputs >= num_outputs: + return complex_response_residual_freq + else: + if np.all(residues_lstsq <= tol_residual_existence): + return complex_response_residual_freq + else: + raise ValueError( + "A multiplicative input uncertainty residual does not exist for the " + "given nominal and off-nominal frequency response matrix. This is due " + "to the fact that the multiplicative input uncertainty model cannot " + "account for all off-nominal models when the number of inputs is less " + "than the number of outputs. Please consider using a different " + "uncertainty model." + ) + + +def _compute_uncertainty_residual_multiplicative_output_freq( + complex_response_nom_freq: np.ndarray, + complex_response_offnom_freq: np.ndarray, + tol_residual_existence: float = 1e-8, +) -> np.ndarray: + """Compute the multiplicative output uncertainty residual at a frequency. + + Parameters + ---------- + complex_response_nom_freq : np.ndarray + Nominal frequency response matrix evaluated at a given frequency. + complex_response_offnom_freq : np.ndarray + Frequency response matrix of a single off-nominal system evaluated at a given + frequency. + tol_residual_existence : float + Tolerance for the existence of a multiplicative input uncertainty residual. + + Returns + ------- + np.ndarray + Multiplicative output uncertainty residual at a given frequency. + + Raises + ------ + ValueError + An output multiplicative uncertainty residual does not exist at the given + frequency as the linear system used to compute the residual does not have a + solution. + """ + + num_inputs = complex_response_nom_freq.shape[1] + num_outputs = complex_response_nom_freq.shape[0] + if num_inputs > num_outputs: + warnings.warn( + "Multipliative output uncertainty models cannot include all possible " + "off-nominal systems when the number of outputs is less than the number of " + "inputs. This may result in an error if an uncertainty residual cannot " + "be found for a given nominal and off-nominal frequency response matrix." + ) + A = complex_response_nom_freq.T + B = complex_response_offnom_freq.T - complex_response_nom_freq.T + X, residues_lstsq, _, _ = scipy.linalg.lstsq(A, B) + complex_response_residual_freq = X.T + + if num_inputs <= num_outputs: + return complex_response_residual_freq + else: + if np.all(residues_lstsq <= tol_residual_existence): + return complex_response_residual_freq + else: + raise ValueError( + "A multiplicative output uncertainty residual does not exist for the " + "given nominal and off-nominal frequency response matrix. This is " + "due to the fact that the multiplicative output uncertainty model " + "cannot account for all off-nominal models when the number of outputs " + "is less than the number of inputs. Please consider using a different " + "uncertainty model." + ) + + +def _compute_uncertainty_residual_inverse_additive_freq( + complex_response_nom_freq: np.ndarray, + complex_response_offnom_freq: np.ndarray, + tol_residual_existence: float = 1e-8, +) -> np.ndarray: + """Compute the inverse additive uncertainty residual at a frequency. + + Parameters + ---------- + complex_response_nom_freq : np.ndarray + Nominal frequency response matrix evaluated at a given frequency. + complex_response_offnom_freq : np.ndarray + Frequency response matrix of a single off-nominal system evaluated at a given + frequency. + tol_residual_existence : float + Tolerance for the existence of a inverse additive uncertainty residual. + + Returns + ------- + np.ndarray + Inverse additive uncertainty residual at a given frequency. + + Raises + ------ + ValueError + An inverse additive uncertainty residual does not exist at the given + frequency as the linear system used to compute the residual does not have a + solution. + """ + + num_inputs = complex_response_nom_freq.shape[1] + num_outputs = complex_response_nom_freq.shape[0] + if num_inputs != num_outputs: + warnings.warn( + "Inverse additive uncertainty models cannot include all possible " + "off-nominal systems when the number of inputs is not equal to the number " + "of outputs. This may result in an error if an uncertainty residual cannot " + "be found for a given nominal and off-nominal frequency response matrix." + ) + + A1 = complex_response_offnom_freq + B1 = complex_response_offnom_freq - complex_response_nom_freq + Y, residues_lstsq_1, _, _ = scipy.linalg.lstsq(A1, B1) + A2 = complex_response_nom_freq.T + B2 = Y.T + X, residues_lstsq_2, _, _ = scipy.linalg.lstsq(A2, B2) + complex_response_residual_freq = X.T + + if num_inputs == num_outputs: + return complex_response_residual_freq + else: + if np.all(residues_lstsq_1 <= tol_residual_existence) and np.all( + residues_lstsq_2 <= tol_residual_existence + ): + return complex_response_residual_freq + else: + raise ValueError( + "An inverse additive uncertainty residual does not exist for the " + "given nominal and off-nominal frequency response matrix. This is " + "due to the fact that the inverse additive uncertainty model " + "cannot account for all off-nominal models when the number of inputs " + "not equal to the number of outputs. Please consider using a different " + "uncertainty model." + ) + + +def _compute_uncertainty_residual_inverse_multiplicative_input_freq( + complex_response_nom_freq: np.ndarray, + complex_response_offnom_freq: np.ndarray, + tol_residual_existence: float = 1e-8, +) -> np.ndarray: + """Compute the inverse multiplicative input uncertainty residual at a frequency. + + Parameters + ---------- + complex_response_nom_freq: np.ndarray + Nominal frequency response matrix evaluated at a given frequency. + complex_response_offnom_freq : np.ndarray + Frequency response matrix of a single off-nominal system evaluated at a given + frequency. + tol_residual_existence : float + Tolerance for the existence of an inverse multiplicative input uncertainty + residual. + + Returns + ------- + np.ndarray + Inverse multiplicative input uncertainty residual at a given frequency. + + Raises + ------ + ValueError + An inverse multiplicative input uncertainty residual does not exist at the + given frequency as the linear system used to compute the residual does not have + a solution. + """ + + num_inputs = complex_response_nom_freq.shape[1] + num_outputs = complex_response_nom_freq.shape[0] + if num_inputs < num_outputs: + warnings.warn( + "Inverse multipliative input uncertainty models cannot include all " + "possible off-nominal systems when the number of inputs is less than the " + "number of outputs. This may result in an error if an uncertainty residual " + "cannot be found for a given nominal and off-nominal frequency response " + "matrix." + ) + + A = complex_response_offnom_freq + B = complex_response_offnom_freq - complex_response_nom_freq + X, residues_lstsq, _, _ = scipy.linalg.lstsq(A, B) + complex_response_residual_freq = X + + if num_inputs >= num_outputs: + return complex_response_residual_freq + else: + if np.all(residues_lstsq <= tol_residual_existence): + return complex_response_residual_freq + else: + raise ValueError( + "An inverse multiplicative input uncertainty residual does not exist " + "for the given nominal and off-nominal frequency response matrix. This " + "is due to the fact that the inverse multiplicative input uncertainty " + "model cannot account for all off-nominal models when the number of " + "inputs is less than the number of outputs. Please consider using a " + "different uncertainty model." + ) + + +def _compute_uncertainty_residual_inverse_multiplicative_output_freq( + complex_response_nom_freq: np.ndarray, + complex_response_offnom_freq: np.ndarray, + tol_residual_existence: float = 1e-8, +) -> np.ndarray: + """Compute the inverse multiplicative output uncertainty residual at a frequency. + + Parameters + ---------- + complex_response_nom_freq : np.ndarray + Nominal frequency response matrix evaluated at a given frequency. + complex_response_offnom_freq : np.ndarray + Frequency response matrix of a single off-nominal system evaluated at a given + frequency. + tol_residual_existence : float + Tolerance for the existence of a multiplicative input uncertainty residual. + + Returns + ------- + np.ndarray + Inverse multiplicative output uncertainty residual at a given frequency. + + Raises + ------ + ValueError + An inverse multiplicative output uncertainty residual does not exist at the + given frequency as the linear system used to compute the residual does not have + a solution. + """ + + num_inputs = complex_response_nom_freq.shape[1] + num_outputs = complex_response_nom_freq.shape[0] + if num_inputs > num_outputs: + warnings.warn( + "Inverse multipliative output uncertainty models cannot include all " + "possible off-nominal systems when the number of outputs is less than the " + "number of inputs. This may result in an error if an uncertainty residual " + "cannot be found for a given nominal and off-nominal frequency response " + "matrix." + ) + A = complex_response_offnom_freq.T + B = complex_response_offnom_freq.T - complex_response_nom_freq.T + X, residues_lstsq, _, _ = scipy.linalg.lstsq(A, B) + complex_response_residual_freq = X.T + + if num_inputs <= num_outputs: + return complex_response_residual_freq + else: + if np.all(residues_lstsq <= tol_residual_existence): + return complex_response_residual_freq + else: + raise ValueError( + "An inverse multiplicative output uncertainty residual does not exist " + "for the given nominal and off-nominal frequency response matrix. This " + "is due to the fact that the inverse multiplicative output uncertainty " + "model cannot account for all off-nominal models when the number of " + "outputs is less than the number of inputs. Please consider using a " + "different uncertainty model." + ) + + +def compute_optimal_uncertainty_weight_response( + complex_response_residual_list: Union[np.ndarray, control.FrequencyResponseList], + weight_left_structure: str, + weight_right_structure: str, + solver_params: Optional[Dict[str, Any]] = None, +) -> Tuple[np.ndarray, np.ndarray]: + """Compute the optimal uncertainty weight frequency response. + + Parameters + ---------- + complex_response_residual_list : Union[np.ndarray, control.FrequencyResponseList] + Frequency response of the residuals for which to compute the optimal uncertainty + weights. + weight_left_structure : str + Structure of the left uncertainty weight. Valid structures include: "diagonal", + "scalar", and "identity". + weight_right_structure : str + Structure of the right uncertainty weight. Valid structures include: "diagonal", + "scalar", and "identity". + solver_param : Dict[str, Any] + Keyword arguments for the convex optimization solver. See [#cvxpy_solver]_ for + more information. + + Returns + ------- + Tuple[np.ndarray, np.ndarray] + Frequency response of the diagonal elements of the left and right uncertainty + weights. + + Examples + -------- + Compute the optimal left and right uncertainty weights for the multiplicative input + residual frequency response where the left and right uncertainty weights are assumed + to be diagonal. + + >>> complex_response_nom, complex_response_offnom_list, omega = ( + ... example_multimodel_uncertainty + ... ) + >>> uncertainty_models = { + ... "additive", + ... "multiplicative_input", + ... "multiplicative_output", + ... "inverse_additive", + ... "inverse_multiplicative_input", + ... "inverse_multiplicative_output", + ... } + >>> complex_response_residual_dict = compute_uncertainty_residual_response( + ... complex_response_nom, + ... complex_response_offnom_list, + ... uncertainty_models, + ... ) + >>> complex_response_weight_left, complex_response_weight_right = ( + ... dkpy.compute_optimal_uncertainty_weight_response( + ... complex_response_residual_dict["multiplicative_input"], + ... "diagonal", + ... "diagonal", + ... ) + ... ) + """ + + # Convert frequency response data to expected type + complex_response_residual_list = _convert_frequency_response_list_to_array( + complex_response_residual_list + ) + + # Frequency response parameters + num_frequency = complex_response_residual_list.shape[1] + + # Compute optimal uncertainty weights + complex_response_weight_left = [] + complex_response_weight_right = [] + for idx_freq in range(num_frequency): + complex_residual_freq = complex_response_residual_list[:, idx_freq, :, :] + weight_left_freq, weight_right_freq = _compute_optimal_weight_freq( + complex_residual_freq, + weight_left_structure, + weight_right_structure, + solver_params, + ) + complex_response_weight_left.append(weight_left_freq) + complex_response_weight_right.append(weight_right_freq) + + # Generate uncertainty weight complex frequency reponses + complex_response_weight_left = np.array(complex_response_weight_left) + complex_response_weight_right = np.array(complex_response_weight_right) + + return complex_response_weight_left, complex_response_weight_right + + +def _compute_optimal_weight_freq( + complex_residual_offnom_set_freq: np.ndarray, + weight_left_structure: str, + weight_right_structure: str, + solver_params: Optional[Dict[str, Any]] = None, +) -> Tuple[np.ndarray, np.ndarray]: + """Compute the optimal uncertainty weight at a given frequency. + + Parameters + ---------- + complex_residual_offnom_set_freq : np.ndarray + Frequency response matrix of the off-nominal models at a given frequency. + weight_left_structure : str + Structure of the left uncertainty weight. Valid structures include: "diagonal", + "scalar", and "identity". + weight_right_structure : str + Structure of the right uncertainty weight. Valid structures include: "diagonal", + "scalar", and "identity". + solver_param : Dict[str, Any] + Keyword arguments for the convex optimization solver. See [#cvxpy_solver]_ for + more information. + + Returns + ------- + Tuple[np.ndarray, np.ndarray] + The left and right uncertainty weight frequency response matrices at a given + frequency. + + References + ---------- + .. [#cxvpy_solver] https://www.cvxpy.org/tutorial/solvers/index.html + """ + + # Solver settings + solver_params = ( + { + "solver": cvxpy.CLARABEL, + "tol_gap_abs": 1e-6, + "tol_gap_rel": 1e-6, + "tol_feas": 1e-6, + "tol_infeas_abs": 1e-6, + "tol_infeas_rel": 1e-6, + } + if solver_params is None + else solver_params + ) + + # System parameters + num_left = complex_residual_offnom_set_freq.shape[1] + num_right = complex_residual_offnom_set_freq.shape[2] + num_offnom = complex_residual_offnom_set_freq.shape[0] + + # Generate left weight variable + if weight_left_structure == "diagonal": + L = cvxpy.Variable((num_left, num_left), diag=True) + elif weight_left_structure == "scalar": + L_scalar = cvxpy.Variable() + L = L_scalar * scipy.sparse.eye_array(num_left) + elif weight_left_structure == "identity": + L = cvxpy.Parameter( + shape=(num_left, num_left), value=np.eye(num_left), diag=True + ) + else: + raise ValueError( + f'"{weight_left_structure}" is not a valid value for ' + '`weight_right_structure`. It must take a value of either "diagonal" ' + '"scalar", or "identity".' + ) + + # Generate right weight variable + if weight_right_structure == "diagonal": + R = cvxpy.Variable((num_right, num_right), diag=True) + elif weight_right_structure == "scalar": + R_scalar = cvxpy.Variable() + R = R_scalar * scipy.sparse.eye_array(num_right) + elif weight_right_structure == "identity": + R = cvxpy.Parameter( + shape=(num_right, num_right), value=np.eye(num_right), diag=True + ) + else: + raise ValueError( + f'"{weight_right_structure}" is not a valid value for ' + '`weight_right_structure`. It must take a value of either "diagonal" ' + '"scalar", or "identity".' + ) + + # Generate optimal uncertainty weight constraints + constraint_freq_list = [] + for idx_offnom in range(num_offnom): + E_k = cvxpy.Parameter( + shape=(num_left, num_right), + value=complex_residual_offnom_set_freq[idx_offnom, :, :], + complex=True, + ) + constraint_matrix_freq = cvxpy.bmat( + [ + [L, E_k], + [E_k.H, R], + ] + ) + constraint_freq = constraint_matrix_freq >> 0 + constraint_freq_list.append(constraint_freq) + constraint_freq_list.append(L.H == L) + constraint_freq_list.append(L >> 0) + constraint_freq_list.append(R.H == R) + constraint_freq_list.append(R >> 0) + + # Semidefinite program + objective = cvxpy.Minimize(cvxpy.trace(L) + cvxpy.trace(R)) + problem = cvxpy.Problem(objective, constraint_freq_list) + problem.solve(**solver_params) + + # Extract left weight + if weight_left_structure == "identity": + L_value = np.array(L.value) + else: + L_value = np.array(L.value.toarray()) + complex_response_weight_left_freq = np.sqrt(L_value) + # Extract right weight + if weight_right_structure == "identity": + R_value = np.array(R.value) + else: + R_value = np.array(R.value.toarray()) + complex_response_weight_right_freq = np.sqrt(R_value) + + return complex_response_weight_left_freq, complex_response_weight_right_freq + + +def fit_overbounding_uncertainty_weight( + complex_response_uncertainty_weight: Union[ + np.ndarray, control.FrequencyResponseData + ], + omega: np.ndarray, + order: Union[int, List[int], np.ndarray], + weight: Optional[np.ndarray] = None, + linear_solver_params: Optional[Dict[str, Any]] = None, + tol_bisection: float = 1e-3, + max_iter_bisection: int = 500, + num_spec_constr: int = 500, +) -> control.StateSpace: + """ + Fit an overbounding stable and minimum-phase state-space uncertainty weight to + frequency response data. + + Parameters + ---------- + complex_response_uncertainty_weight : Union[np.ndarray, control.FrequencyResponseData] + Uncertainty weight frequency response used for the overbounding fit. + order : Union[int, List[int], np.ndarray] + Order of the LTI system fit. If `order` is an `int`, the order will be + used for all elements of the weight. If `order` is a `List` or `np.ndarray`, + the order can be specified for each element of the weight. + weight : Optional[np.ndarray] = None + Frequency-dependent weight used to improve the fit over certain bandwidths. The + weight is a 2D array with the first dimension representing the number of + elements in the weight and the second dimension representing the number of + frequency points. + linear_solver_params : Dict[str, Any] + Keyword arguments for the linear feasibility problem solver. See + [#cvxpy_solver]_ for more information. + tol_bisection : float + Numerical tolerance for the bisection algorithm. + max_iter_bisection : int + Maximum allowable number of iterations in the bisection algorithm. + num_spec_constr : int + Number of constraints used to enforce the spectral factorizability of the + fitted autocorrelation. + + Returns + ------- + control.StateSpace + Fitted overbounding uncertainty weight state-space systems. + + Examples + -------- + Fit an overbounding stable and minimum phase system to the left and right + uncertainty weights for the multiplicative input uncertainty model. + + >>> complex_response_nom, complex_response_offnom_list, omega = ( + ... example_multimodel_uncertainty + ... ) + >>> uncertainty_models = { + ... "additive", + ... "multiplicative_input", + ... "multiplicative_output", + ... "inverse_additive", + ... "inverse_multiplicative_input", + ... "inverse_multiplicative_output", + ... } + >>> complex_response_residual_dict = compute_uncertainty_residual_response( + ... complex_response_nom, + ... complex_response_offnom_list, + ... uncertainty_models, + ... ) + >>> complex_response_weight_left, complex_response_weight_right = ( + ... dkpy.compute_optimal_uncertainty_weight_response( + ... complex_response_residual_dict["multiplicative_input"], + ... "diagonal", + ... "diagonal", + ... ) + ... ) + >>> weight_left = dkpy.fit_overbounding_uncertainty_weight( + ... complex_response_weight_left, omega, [4, 5] + ... ) + >>> weight_right = dkpy.fit_overbounding_uncertainty_weight( + ... complex_response_weight_right, omega, [3, 5] + ... ) + + References + ---------- + .. [#cxvpy_solver] https://www.cvxpy.org/tutorial/solvers/index.html + """ + + # Convert frequency response data to expected type + complex_response_uncertainty_weight = _convert_frequency_response_data_to_array( + complex_response_uncertainty_weight + ) + + # Solver settings + linear_solver_params = ( + { + "solver": cvxpy.CLARABEL, + "tol_gap_abs": 1e-9, + "tol_gap_rel": 1e-9, + "tol_feas": 1e-9, + "tol_infeas_abs": 1e-9, + "tol_infeas_rel": 1e-9, + } + if linear_solver_params is None + else linear_solver_params + ) + + # Parse arguments + num_elements = complex_response_uncertainty_weight.shape[1] + order_list = ( + order * np.ones(num_elements, dtype=int) + if isinstance(order, int) + else np.array(order, dtype=int) + ) + + if weight is None: + # Take the default frequency-dependent weight as the normalized magnitude the + # uncertainty weight magnitude in order to place greater importance on tightly + # overbounding at the largest uncertainties + weight = np.diagonal( + np.abs(complex_response_uncertainty_weight), axis1=1, axis2=2 + ) + weight = weight / np.max(weight, axis=0) + + uncertainty_weight_list = [] + for idx_element in range(num_elements): + # Extract the parameters relevant to each SISO uncertainty weight element + magnitude_response_weight_element = np.abs( + complex_response_uncertainty_weight[:, idx_element, idx_element] + ) + order_element = order_list[idx_element] + weight_element = weight[:, idx_element] + + # Fit the uncertainty weight to each SISO element + uncertainty_weight_element = utilities._fit_magnitude_log_chebyshev_siso( + omega=omega, + magnitude_fit=magnitude_response_weight_element, + order=order_element, + magnitude_lower_bound=magnitude_response_weight_element, + weight=weight_element, + linear_solver_params=linear_solver_params, + tol_bisection=tol_bisection, + max_iter_bisection=max_iter_bisection, + num_spec_constr=num_spec_constr, + ) + uncertainty_weight_list.append(uncertainty_weight_element) + + # Construct uncertainty weight fit from SISO elements + uncertainty_weight = control.append(*uncertainty_weight_list) + + return uncertainty_weight + + +def _convert_frequency_response_data_to_array( + frequency_response: Union[ + np.ndarray, np.typing.ArrayLike, control.FrequencyResponseData + ], +) -> np.ndarray: + if isinstance(frequency_response, control.FrequencyResponseData): + complex_response = np.array(frequency_response.complex, dtype=complex) + + # SISO `FrequencyResponseData` objects return 1D arrays for frequency response + # data and 3D arrays for MIMO systems. `dkpy` uses a 3D array in all cases, so + # the 1D array is converted into a 3D array. + if complex_response.ndim == 1: + complex_response = complex_response[None, None, :] + + # The `control.FrequencyResponseData.complex` uses an array with shape + # (output, input, frequency) whereas the format used in `dkpy` is + # (frequency, output, input). Therefore, a transpose to switch the axes is + # required. + complex_response = complex_response.transpose(2, 0, 1) + else: + complex_response = np.array(frequency_response, dtype=complex) + if complex_response.ndim != 3: + raise ValueError( + f"The dimension of `complex_response` is {complex_response.ndim}, " + "whereas it should be 3 (frequency, output, input)." + ) + return complex_response + + +def _convert_frequency_response_list_to_array( + frequency_response_list: Union[ + np.ndarray, np.typing.ArrayLike, control.FrequencyResponseList + ], +) -> np.ndarray: + if isinstance(frequency_response_list, control.FrequencyResponseList): + complex_response_list = [] + for frequency_response in frequency_response_list: + complex_response = np.array(frequency_response.complex, dtype=complex) + + # SISO `FrequencyResponseData` objects return 1D arrays for frequency response + # data and 3D arrays for MIMO systems. `dkpy` uses a 3D array in all cases, so + # the 1D array is converted into a 3D array. + if complex_response.ndim == 1: + complex_response = complex_response[None, None, :] + + # The `control.FrequencyResponseData.complex` uses an array with shape + # (output, input, frequency) whereas the format used in `dkpy` is + # (frequency, output, input). Therefore, a transpose to switch the axes is + # required. + complex_response = complex_response.transpose(2, 0, 1) + complex_response_list.append(complex_response) + complex_response_list = np.array(complex_response_list, dtype=complex) + else: + complex_response_list = np.array(frequency_response_list, dtype=complex) + if complex_response_list.ndim != 4: + raise ValueError( + "The dimension of `complex_response_list` is " + f"{complex_response_list.ndim}, whereas it should be 4 (off-nominals, " + "frequency, output, input)." + ) + return complex_response_list + + +def plot_magnitude_response_uncertain_model_set( + complex_response_nom: np.ndarray, + complex_response_offnom_list: np.ndarray, + omega: np.ndarray, + db: bool = True, + hz: bool = False, + frequency_log_scale: bool = True, + plot_nom_kw: Dict[str, Any] = {}, + plot_offnom_kw: Dict[str, Any] = {}, + subplot_kw: Dict[str, Any] = {}, +) -> Tuple[Figure, Union[Axes, np.ndarray]]: + """Plot magnitude response of nominal model and set of off-nominal models. + + Parameters + ---------- + complex_response_nom : np.ndarray + Frequency response matrices over a grid of frequencies of the nominal system. + complex_response_offnom_list : np.ndarray + Frequency response matrices over a grid of frequencies of the set of off-nominal + systems. + omega : np.ndarray + Angular frequency grid. + db : bool + If True, plot the magnitude in units of dB. Otherwise, plot the magnitude in + absolute units. + hz : bool + If True, plot the frequency in units of Hz. Otherwise, plot the frequency in + units of rad/s. + frequency_log_scale : bool + If True, plot the frequency using a logarithmic axis. Otherwise, plot the + the frequency using a linear axis. + plot_nom_kw : Dict[str, Any] + Keyword arguments for the nominal model plot. See [#plot_kw]_ for more + information on plotting keywords. + plot_offnom_kw : Dict[str, Any] + Keyword arguments for the off-nominal model plot. See [#plot_kw]_ for more + information on plotting keywords. + subplot_kw : Dict[str, Any] + Keyword arguments for the subplot. See [#subplot_kw]_ for more information on + the subplot keywords. + + References + ---------- + .. [#plot_kw] https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html + .. [#subplot_kw] https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.subplots.html + """ + + # System paramters + num_offnom = complex_response_offnom_list.shape[0] + num_outputs = complex_response_offnom_list.shape[2] + num_inputs = complex_response_offnom_list.shape[3] + + # Nominal model plot keyword arguments + plot_nom_kwargs = {"color": "C0", "label": "Nominal"} + plot_nom_kwargs.update(plot_nom_kw) + + # Off-nominal model plot keyword arguments + plot_offnom_kwargs = {"color": "C1", "alpha": 0.25, "label": "Off-Nominal"} + plot_offnom_kwargs.update(plot_offnom_kw) + + # Subplot keyword arguments + subplot_kwargs = {"sharex": True, "layout": "constrained"} + subplot_kwargs.update(subplot_kw) + + # Initialize figure + fig, ax = plt.subplots(num_outputs, num_inputs, squeeze=False, **subplot_kwargs) + + # Off-nominal frequency response + for idx_offnom in range(num_offnom): + magnitude_offnom = np.abs(complex_response_offnom_list[idx_offnom, :, :, :]) + for idx_input in range(num_inputs): + for idx_output in range(num_outputs): + ax[idx_output, idx_input].plot( + omega / (2 * np.pi) if hz else omega, + control.mag2db(magnitude_offnom[:, idx_output, idx_input]) + if db + else magnitude_offnom[:, idx_output, idx_input], + **plot_offnom_kwargs, + ) + + # Nominal frequency response + magnitude_nom = np.abs(complex_response_nom) + for idx_input in range(num_inputs): + for idx_output in range(num_outputs): + ax[idx_output, idx_input].plot( + omega / (2 * np.pi) if hz else omega, + control.mag2db(magnitude_nom[:, idx_output, idx_input]) + if db + else magnitude_nom[:, idx_output, idx_input], + **plot_nom_kwargs, + ) + + # Plot settings + for ax_output in ax: + for ax_output_input in ax_output: + ax_output_input.set_ylabel("Magnitude (dB)" if db else "Magnitude (-)") + ax_output_input.grid() + if frequency_log_scale: + ax_output_input.set_xscale("log") + for idx_input in range(num_inputs): + ax[-1, idx_input].set_xlabel("$f$ (Hz)" if hz else r"$\omega$ (rad/s)") + handles, labels = ax[0, 0].get_legend_handles_labels() + legend_dict = dict(zip(labels, handles)) + fig.legend( + labels=legend_dict.keys(), + handles=legend_dict.values(), + loc="outside lower center", + ncol=2, + ) + + return fig, ax + + +def plot_phase_response_uncertain_model_set( + complex_response_nom: np.ndarray, + complex_response_offnom_list: np.ndarray, + omega: np.ndarray, + deg: bool = True, + hz: bool = False, + frequency_log_scale: bool = True, + plot_nom_kw: Dict[str, Any] = {}, + plot_offnom_kw: Dict[str, Any] = {}, + subplot_kw: Dict[str, Any] = {}, +) -> Tuple[Figure, Union[Axes, np.ndarray]]: + """Plot phase response of nominal model and set of off-nominal models. + + Parameters + ---------- + complex_response_nom : np.ndarray + Frequency response matrices over a grid of frequencies of the nominal system. + complex_response_offnom_list : np.ndarray + Frequency response matrices over a grid of frequencies of the set of off-nominal + systems. + omega : np.ndarray + Angular frequency grid. + db : bool + If True, plot the magnitude in units of dB. Otherwise, plot the magnitude in + absolute units. + hz : bool + If True, plot the frequency in units of Hz. Otherwise, plot the frequency in + units of rad/s. + frequency_log_scale : bool + If True, plot the frequency using a logarithmic axis. Otherwise, plot the + the frequency using a linear axis. + plot_nom_kw : Dict[str, Any] + Keyword arguments for the nominal model plot. See [#plot_kw]_ for more + information on plotting keywords. + plot_offnom_kw : Dict[str, Any] + Keyword arguments for the off-nominal model plot. See [#plot_kw]_ for more + information on plotting keywords. + subplot_kw : Dict[str, Any] + Keyword arguments for the subplot. See [#subplot_kw]_ for more information on + the subplot keywords. + + References + ---------- + .. [#plot_kw] https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html + .. [#subplot_kw] https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.subplots.html + """ + + # System paramters + num_offnom = complex_response_offnom_list.shape[0] + num_outputs = complex_response_offnom_list.shape[2] + num_inputs = complex_response_offnom_list.shape[3] + + # Nominal model plot keyword arguments + plot_nom_kwargs = {"color": "C0", "label": "Nominal"} + plot_nom_kwargs.update(plot_nom_kw) + + # Off-nominal model plot keyword arguments + plot_offnom_kwargs = {"color": "C1", "alpha": 0.25, "label": "Off-Nominal"} + plot_offnom_kwargs.update(plot_offnom_kw) + + # Subplot keyword arguments + subplot_kwargs = {"sharex": True, "layout": "constrained"} + subplot_kwargs.update(subplot_kw) + + # Initialize figure + fig, ax = plt.subplots(num_outputs, num_inputs, squeeze=False, **subplot_kwargs) + + # Off-nominal frequency response + for idx_offnom in range(num_offnom): + phase_offnom = np.angle(complex_response_offnom_list[idx_offnom, :, :, :]) + for idx_input in range(num_inputs): + for idx_output in range(num_outputs): + ax[idx_output, idx_input].plot( + omega / (2 * np.pi) if hz else omega, + 180 / np.pi * phase_offnom[:, idx_output, idx_input] + if deg + else phase_offnom[:, idx_output, idx_input], + **plot_offnom_kwargs, + ) + + # Nominal frequency response + phase_nom = np.angle(complex_response_nom) + for idx_input in range(num_inputs): + for idx_output in range(num_outputs): + ax[idx_output, idx_input].plot( + omega / (2 * np.pi) if hz else omega, + 180 / np.pi * phase_nom[:, idx_output, idx_input] + if deg + else phase_nom[:, idx_output, idx_input], + **plot_nom_kwargs, + ) + + # Plot settings + for ax_output in ax: + for ax_output_input in ax_output: + ax_output_input.set_ylabel("Phase (deg)" if deg else "Phase (rad)") + ax_output_input.grid() + if frequency_log_scale: + ax_output_input.set_xscale("log") + for idx_input in range(num_inputs): + ax[-1, idx_input].set_xlabel("$f$ (Hz)" if hz else r"$\omega$ (rad/s)") + handles, labels = ax[0, 0].get_legend_handles_labels() + legend_dict = dict(zip(labels, handles)) + fig.legend( + labels=legend_dict.keys(), + handles=legend_dict.values(), + loc="outside lower center", + ncol=2, + ) + + return fig, ax + + +def plot_singular_value_response_uncertain_model_set( + complex_response_nom: np.ndarray, + complex_response_offnom_list: np.ndarray, + omega: np.ndarray, + db: bool = True, + hz: bool = False, + frequency_log_scale: bool = True, + plot_nom_kw: Dict[str, Any] = {}, + plot_offnom_kw: Dict[str, Any] = {}, + subplot_kw: Dict[str, Any] = {}, +) -> Tuple[Figure, Union[Axes, np.ndarray]]: + """Plot singular value response of nominal model and set of off-nominal models. + + Parameters + ---------- + complex_response_nom : np.ndarray + Frequency response matrices over a grid of frequencies of the nominal system. + complex_response_offnom_list : np.ndarray + Frequency response matrices over a grid of frequencies of the set of off-nominal + systems. + omega : np.ndarray + Angular frequency grid. + db : bool + If True, plot the magnitude in units of dB. Otherwise, plot the magnitude in + absolute units. + hz : bool + If True, plot the frequency in units of Hz. Otherwise, plot the frequency in + units of rad/s. + frequency_log_scale : bool + If True, plot the frequency using a logarithmic axis. Otherwise, plot the + the frequency using a linear axis. + plot_nom_kw : Dict[str, Any] + Keyword arguments for the nominal model plot. See [#plot_kw]_ for more + information on plotting keywords. + plot_offnom_kw : Dict[str, Any] + Keyword arguments for the off-nominal model plot. See [#plot_kw]_ for more + information on plotting keywords. + subplot_kw : Dict[str, Any] + Keyword arguments for the subplot. See [#subplot_kw]_ for more information on + the subplot keywords. + + References + ---------- + .. [#plot_kw] https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html + .. [#subplot_kw] https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.subplots.html + """ + + # System paramters + num_offnom = complex_response_offnom_list.shape[0] + + # Nominal model plot keyword arguments + plot_nom_kwargs = {"color": "C0", "label": "Nominal"} + plot_nom_kwargs.update(plot_nom_kw) + + # Off-nominal model plot keyword arguments + plot_offnom_kwargs = {"color": "C1", "alpha": 0.25, "label": "Off-Nominal"} + plot_offnom_kwargs.update(plot_offnom_kw) + + # Subplot keyword arguments + subplot_kwargs = {"ncols": 1, "nrows": 1, "layout": "constrained"} + subplot_kwargs.update(subplot_kw) + + # Initialize figure + fig, ax = plt.subplots(**subplot_kwargs) + + # Off-nominal frequency response + for idx_offnom in range(num_offnom): + sval_offnom = np.linalg.svdvals( + complex_response_offnom_list[idx_offnom, :, :, :] + ) + for idx_sval in range(sval_offnom.shape[1]): + ax.semilogx( + omega / (2 * np.pi) if hz else omega, + control.mag2db(sval_offnom[:, idx_sval]) + if db + else sval_offnom[:, idx_sval], + **plot_offnom_kwargs, + ) + + # Nominal frequency response + sval_nom = np.linalg.svdvals(complex_response_nom) + for idx_sval in range(sval_nom.shape[1]): + ax.plot( + omega / (2 * np.pi) if hz else omega, + control.mag2db(sval_nom[:, idx_sval]) if db else sval_nom[:, idx_sval], + **plot_nom_kwargs, + ) + + # Plot settings + if frequency_log_scale: + ax.set_xscale("log") + ax.set_ylabel( + "Singular Value Magnitude (dB)" if db else " Singular Value Magnitude (-)" + ) + ax.grid() + ax.set_xlabel("$f$ (Hz)" if hz else r"$\omega$ (rad/s)") + handles, labels = ax.get_legend_handles_labels() + legend_dict = dict(zip(labels, handles)) + fig.legend( + labels=legend_dict.keys(), + handles=legend_dict.values(), + loc="outside lower center", + ncol=2, + ) + + return fig, ax + + +def plot_singular_value_response_residual( + complex_response_residual_dict: Dict[str, np.ndarray], + omega: np.ndarray, + db: bool = True, + hz: bool = False, + frequency_log_scale: bool = True, + plot_sval_max_kw: Dict[str, Any] = {}, + plot_sval_kw: Dict[str, Any] = {}, + subplot_kw: Dict[str, Any] = {}, +) -> Dict[str, Tuple[Figure, Union[Axes, np.ndarray]]]: + """Plot the singular value response of the uncertainty residuals for different + unstructured uncertainty models on separate figures. + + Parameters + ---------- + complex_response_residual_dict : Dict[str, np.ndarray] + Dictionary of the uncertainty residual frequency response matrices over a grid + of frequencies for different uncertainty models. + omega : np.narray + Angular frequency grid. + db : bool + If True, plot the magnitude in units of dB. Otherwise, plot the magnitude in + absolute units. + hz : bool + If True, plot the frequency in units of Hz. Otherwise, plot the frequency in + units of rad/s. + frequency_log_scale : bool + If True, plot the frequency using a logarithmic axis. Otherwise, plot the + the frequency using a linear axis. + plot_sval_max_kw : Dict[str, Any] + Keyword arguments for the maximum singular value plot. See [#plot_kw]_ for more + information on plotting keywords. + plot_sval_kw : Dict[str, Any] + Keyword arguments for the singular value plots. See [#plot_kw]_ for more + information on plotting keywords. + subplot_kw : Dict[str, Any] + Keyword arguments for the subplot. See [#subplot_kw]_ for more information on + the subplot keywords. + + References + ---------- + .. [#plot_kw] https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html + .. [#subplot_kw] https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.subplots.html + """ + + # Residual singular value plot keyword arguments + plot_sval_max_kwargs = {"color": "black"} + plot_sval_max_kwargs.update(plot_sval_max_kw) + + # Residual maximum singular value plot keyword arguments + plot_sval_kwargs = {"color": "grey", "alpha": 0.25} + plot_sval_kwargs.update(plot_sval_kw) + + # Subplot keyword arguments + subplot_kwargs = {"ncols": 1, "nrows": 1, "layout": "constrained"} + subplot_kwargs.update(subplot_kw) + + # Initialize dictionary for storing figure and axes + figure_dict = {} + + # Uncertainty model residual suffixes + residual_suffix_dict = { + "additive": "A", + "multiplicative_input": "I", + "multiplicative_output": "O", + "inverse_additive": "iA", + "inverse_multiplicative_input": "iI", + "inverse_multiplicative_output": "iO", + } + + # Iterate over each uncertainty model + for ( + uncertainty_model_id, + complex_response_residual, + ) in complex_response_residual_dict.items(): + # Off-nominal frequency response parameters + num_offnom = complex_response_residual.shape[0] + + # Uncertainty model suffix + residual_suffix = residual_suffix_dict[uncertainty_model_id] + + # Compute the singular value and maximum singlar value response of the residuals + sval_response_residual = np.linalg.svdvals(complex_response_residual) + sval_max_response_residual = np.max(sval_response_residual, axis=(0, 2)) + + # Intialize the plot + fig, ax = plt.subplots(**subplot_kwargs) + + # Singular value response of the residuals + for idx_offnom in range(num_offnom): + for idx_sval in range(sval_response_residual.shape[2]): + ax.plot( + omega / (2 * np.pi) if hz else omega, + control.mag2db(sval_response_residual[idx_offnom, :, idx_sval]) + if db + else sval_response_residual[idx_offnom, :, idx_sval], + label=rf"$\sigma(E_{{{residual_suffix}}})$", + **plot_sval_kwargs, + ) + # Maximum singular value response of the residuals + ax.plot( + omega / (2 * np.pi) if hz else omega, + control.mag2db(sval_max_response_residual) + if db + else sval_max_response_residual, + label=rf"$\max \; \sigma(E_{{{residual_suffix}}})$", + **plot_sval_max_kwargs, + ) + + # Plot settings + if frequency_log_scale: + ax.set_xscale("log") + ax.set_ylabel( + "Singular Value Magnitude (dB)" if db else "Singular Value Magnitude (-)" + ) + ax.grid() + ax.set_xlabel("$f$ (Hz)" if hz else r"$\omega$ (rad/s)") + handles, labels = ax.get_legend_handles_labels() + legend_dict = dict(zip(labels, handles)) + fig.legend( + labels=legend_dict.keys(), + handles=legend_dict.values(), + loc="outside lower center", + ncol=2, + ) + + figure_dict[uncertainty_model_id] = (fig, ax) + + return figure_dict + + +def plot_singular_value_response_residual_comparison( + complex_response_residual_dict: Dict[str, np.ndarray], + omega: np.ndarray, + db: bool = True, + hz: bool = False, + frequency_log_scale: bool = True, + plot_sval_max_kw: Dict[str, Any] = {}, + subplot_kw: Dict[str, Any] = {}, +) -> Tuple[Figure, Union[Axes, np.ndarray]]: + """Plot the maximum singular value response of the uncertainty residuals for + different unstructured uncertainty models on the same figure for comparision of the + various uncertainty models. + + Parameters + ---------- + complex_response_residual_dict : Dict[str, np.ndarray] + Dictionary of the uncertainty residual frequency response matrices over a grid + of frequencies for different uncertainty models. + omega : np.narray + Angular frequency grid. + db : bool + If True, plot the magnitude in units of dB. Otherwise, plot the magnitude in + absolute units. + hz : bool + If True, plot the frequency in units of Hz. Otherwise, plot the frequency in + units of rad/s. + frequency_log_scale : bool + If True, plot the frequency using a logarithmic axis. Otherwise, plot the + the frequency using a linear axis. + plot_sval_max_kw : Dict[str, Any] + Keyword arguments for the maximum singular value plot. See [#plot_kw]_ for more + information on plotting keywords. + subplot_kw : Dict[str, Any] + Keyword arguments for the subplot. See [#subplot_kw]_ for more information on + the subplot keywords. + + References + ---------- + .. [#plot_kw] https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html + .. [#subplot_kw] https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.subplots.html + """ + + # Residual singular value plot keyword arguments + plot_sval_max_kwargs = {} + plot_sval_max_kwargs.update(plot_sval_max_kw) + + # Subplot keyword arguments + subplot_kwargs = {"ncols": 1, "nrows": 1, "layout": "constrained"} + subplot_kwargs.update(subplot_kw) + + # Uncertainty model suffixes + residual_suffix_dict = { + "additive": "A", + "multiplicative_input": "I", + "multiplicative_output": "O", + "inverse_additive": "iA", + "inverse_multiplicative_input": "iI", + "inverse_multiplicative_output": "iO", + } + + # Maximum singular value response of uncertainty residuals + sval_max_response_residual_dict = {} + for ( + uncertainty_model_id, + complex_response_residual_list, + ) in complex_response_residual_dict.items(): + sval_response_residual_list = np.linalg.svdvals(complex_response_residual_list) + sval_max_response_residual = np.max(sval_response_residual_list, axis=(0, 2)) + sval_max_response_residual_dict[uncertainty_model_id] = ( + sval_max_response_residual + ) + + # Initialize figure + fig, ax = plt.subplots(**subplot_kwargs) + + # Maximum singular value reponse of uncertainty residuals + for ( + uncertainty_model_id, + sval_max_response_residual, + ) in sval_max_response_residual_dict.items(): + residual_suffix = residual_suffix_dict[uncertainty_model_id] + ax.semilogx( + omega / (2 * np.pi) if hz else omega, + control.mag2db(sval_max_response_residual) + if db + else sval_max_response_residual, + label=rf"$\max \; {{\sigma}}(E_{{{residual_suffix}}})$", + **plot_sval_max_kwargs, + ) + + # Plot settings + if frequency_log_scale: + ax.set_xscale("log") + ax.set_ylabel( + "Singular Value Magnitude (dB)" if db else "Singular Value Magnitude (-)" + ) + ax.grid() + ax.set_xlabel("$f$ (Hz)" if hz else r"$\omega$ (rad/s)") + handles, labels = ax.get_legend_handles_labels() + legend_dict = dict(zip(labels, handles)) + fig.legend( + labels=legend_dict.keys(), + handles=legend_dict.values(), + loc="outside lower center", + ncol=3, + ) + + return fig, ax + + +def plot_magnitude_response_uncertainty_weight( + complex_response_weight_left: np.ndarray, + complex_response_weight_right: np.ndarray, + omega: np.ndarray, + weight_left: Optional[control.StateSpace] = None, + weight_right: Optional[control.StateSpace] = None, + db: bool = True, + hz: bool = False, + frequency_log_scale: bool = True, + plot_response_kw: Dict[str, Any] = {}, + plot_response_fit_kw: Dict[str, Any] = {}, + subplot_kw: Dict[str, Any] = {}, +) -> Tuple[Figure, Union[Axes, np.ndarray]]: + """Plot the diagonal elements of the optimal left and right uncertainty weight + frequency responses. Optionally, the fitted overbounding left and right uncertainty + weights can also be displayed. + + Parameters + ---------- + complex_response_weight_left : np.ndarray, + Frequency response matrices of the left uncertainty weight over a grid of + frequencies. + complex_response_weight_right : np.ndarray, + Frequency response matrices of the right uncertainty weight over a grid of + frequencies. + omega : np.ndarray + Angular frequency grid. + weight_left: Optional[control.StateSpace] = None, + State-space model if the fitted overbounding left uncertainty weight. + weight_right: Optional[control.StateSpace] = None, + State-space model if the fitted overbounding right uncertainty weight. + db : bool + If True, plot the magnitude in units of dB. Otherwise, plot the magnitude in + absolute units. + hz : bool + If True, plot the frequency in units of Hz. Otherwise, plot the frequency in + units of rad/s. + frequency_log_scale : bool + If True, plot the frequency using a logarithmic axis. Otherwise, plot the + the frequency using a linear axis. + plot_response_kw : Dict[str, Any] + Keyword arguments for the frequency response of the uncertainty weight plot. + See [#plot_kw]_ for more information on plotting keywords. + plot_response_fit_kw : Dict[str, Any] + Keyword arguments for the frequency response of the fitted uncertainty weight + plot. See [#plot_kw]_ for more information on plotting keywords. + subplot_kw : Dict[str, Any] + Keyword arguments for the subplot. See [#subplot_kw]_ for more information on + the subplot keywords. + + References + ---------- + .. [#plot_kw] https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html + .. [#subplot_kw] https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.subplots.html + """ + + # Uncertainty weight response plot keyword arguments + plot_response_kwargs = { + "color": "C0", + "marker": "*", + "linestyle": " ", + "label": "Response", + } + plot_response_kwargs.update(plot_response_kw) + + # Uncertainty weight fit response plot keyword arguments + plot_response_fit_kwargs = { + "color": "C1", + "linestyle": "-", + "label": "Fit", + } + plot_response_fit_kwargs.update(plot_response_fit_kw) + + # Subplot keyword arguments + subplot_kwargs = {"sharex": True, "layout": "constrained"} + subplot_kw = {} if subplot_kw is None else subplot_kw + subplot_kwargs.update(subplot_kw) + + # Uncertainty weight parameters + num_left = complex_response_weight_left.shape[1] + num_right = complex_response_weight_right.shape[1] + + # Magnitude response of the uncertainty weights + magnitude_response_weight_left = np.abs(complex_response_weight_left) + magnitude_response_weight_right = np.abs(complex_response_weight_right) + + # Initialize figure + fig, ax = plt.subplots( + max(num_left, num_right), 2, sharex=True, layout="constrained" + ) + + # Plot left uncertainty weight frequency response + for idx_left in range(num_left): + ax[idx_left, 0].plot( + omega / (2 * np.pi) if hz else omega, + control.mag2db(magnitude_response_weight_left[:, idx_left, idx_left]) + if db + else magnitude_response_weight_left[:, idx_left, idx_left], + **plot_response_kwargs, + ) + ax[idx_left, 0].set_ylabel( + f"$|W_{{L, ({idx_left + 1}, {idx_left + 1})}}|$ (dB)" + if db + else f"$|W_{{L, ({idx_left + 1}, {idx_left + 1})}}|$ (-)" + ) + ax[idx_left, 0].grid() + # Plot right uncertainty weight frequency response + for idx_right in range(num_left): + ax[idx_right, 1].plot( + omega / (2 * np.pi) if hz else omega, + control.mag2db(magnitude_response_weight_right[:, idx_right, idx_right]) + if db + else magnitude_response_weight_right[:, idx_right, idx_right], + **plot_response_kwargs, + ) + ax[idx_right, 1].set_ylabel( + f"$|W_{{R, ({idx_right + 1}, {idx_right + 1})}}|$ (dB)" + if db + else f"$|W_{{R, ({idx_right + 1}, {idx_right + 1})}}|$ (-)" + ) + ax[idx_right, 1].grid() + + # Plot left uncertainty weight fit frequency response + if weight_left is not None: + response_fit_weight_left = control.frequency_response(weight_left, omega) + magnitude_response_fit_weight_left = np.array( + response_fit_weight_left.magnitude + ) + for idx_left in range(num_left): + ax[idx_left, 0].plot( + omega / (2 * np.pi) if hz else omega, + control.mag2db( + magnitude_response_fit_weight_left[idx_left, idx_left, :] + ) + if db + else magnitude_response_fit_weight_left[idx_left, idx_left, :], + **plot_response_fit_kwargs, + ) + # Plot right uncertainty weight fit frequency response + if weight_right is not None: + response_fit_weight_right = control.frequency_response(weight_right, omega) + magnitude_response_fit_weight_right = np.array( + response_fit_weight_right.magnitude + ) + for idx_right in range(num_right): + ax[idx_right, 1].plot( + omega / (2 * np.pi) if hz else omega, + control.mag2db( + magnitude_response_fit_weight_right[idx_right, idx_right, :] + ) + if db + else magnitude_response_fit_weight_right[idx_right, idx_right, :], + **plot_response_fit_kwargs, + ) + + # Plot settings + for idx_col in range(2): + ax[-1, idx_col].set_xlabel("$f$ (Hz)" if hz else r"$\omega$ (rad/s)") + for ax_row in ax: + for ax_row_col in ax_row: + if frequency_log_scale: + ax_row_col.set_xscale("log") + if not ax_row_col.has_data(): + fig.delaxes(ax_row_col) + handles, labels = ax[0, 0].get_legend_handles_labels() + legend_dict = dict(zip(labels, handles)) + fig.legend( + labels=legend_dict.keys(), + handles=legend_dict.values(), + loc="outside lower center", + ncol=2, + ) + + return fig, ax diff --git a/src/dkpy/utilities.py b/src/dkpy/utilities.py index 2d48e31..4495f97 100644 --- a/src/dkpy/utilities.py +++ b/src/dkpy/utilities.py @@ -3,18 +3,22 @@ __all__ = [ "example_scherer1997_p907", "example_skogestad2006_p325", + "example_mackenroth2004_p430", + "example_multimodel_uncertainty", "_ensure_tf", "_tf_close_coeff", "_auto_lmi_strictness", + "_fit_magnitude_log_chebyshev_siso", ] -from typing import Any, Dict, Union +from typing import Any, Dict, Union, Optional, Tuple +from numpy.typing import ArrayLike import control import cvxpy import numpy as np import scipy.linalg -from numpy.typing import ArrayLike +import warnings def example_scherer1997_p907() -> Dict[str, Any]: @@ -318,6 +322,124 @@ def example_mackenroth2004_p430(): return out +def example_multimodel_uncertainty(): + """Generate uncertain models using parameter variation.""" + + def generate_sys( + omega_n_11, omega_n_22, gamma_11, gamma_22, gain_12, gain_21, tau_12, tau_21 + ): + """Generate example transfer function system from parameters.""" + tf_11 = control.TransferFunction( + [omega_n_11**2], + [1, 2 * gamma_11 * omega_n_11, omega_n_11**2], + ) + tf_12 = control.TransferFunction([gain_12], [tau_12, 1]) + tf_21 = control.TransferFunction([gain_21], [tau_21, 1]) + tf_22 = control.TransferFunction( + [omega_n_22**2], + [1, 2 * gamma_22 * omega_n_22, omega_n_22**2], + ) + tf = control.combine_tf([[tf_11, tf_12], [tf_21, tf_22]]) + + return tf + + # Nominal model paramters + omega_n_11_nom = 5 + omega_n_22_nom = 1 + gamma_11_nom = 0.4 + gamma_22_nom = 0.6 + gain_12_nom = 3 + gain_21_nom = 10 + tau_12_nom = 0.1 + tau_21_nom = 0.5 + # Generate nominal model + sys_nom = generate_sys( + omega_n_11_nom, + omega_n_22_nom, + gamma_11_nom, + gamma_22_nom, + gain_12_nom, + gain_21_nom, + tau_12_nom, + tau_21_nom, + ) + + # Off-nominal system relative parameter variation + rel_variation_omega_n = 0.20 + rel_variation_gamma = 0.20 + rel_variation_gain = 0.05 + rel_variation_tau = 0.25 + num_offnom = 50 + + # Generate off-nominal models + np.random.seed(0) + sys_offnom_list = [] + for _ in range(num_offnom): + # Off-nominal system paramters + omega_n_11_offnom = omega_n_11_nom * ( + 1 + rel_variation_omega_n * (2 * np.random.rand() - 1) + ) + omega_n_22_offnom = omega_n_22_nom * ( + 1 + rel_variation_omega_n * (2 * np.random.rand() - 1) + ) + gamma_11_offnom = gamma_11_nom * ( + 1 + rel_variation_gamma * (2 * np.random.rand() - 1) + ) + gamma_22_offnom = gamma_22_nom * ( + 1 + rel_variation_gamma * (2 * np.random.rand() - 1) + ) + gain_12_offnom = gain_12_nom * ( + 1 + rel_variation_gain * (2 * np.random.rand() - 1) + ) + gain_21_offnom = gain_21_nom * ( + 1 + rel_variation_gain * (2 * np.random.rand() - 1) + ) + tau_12_offnom = tau_12_nom * ( + 1 + rel_variation_tau * (2 * np.random.rand() - 1) + ) + tau_21_offnom = tau_21_nom * ( + 1 + rel_variation_tau * (2 * np.random.rand() - 1) + ) + # Generate off-nominal system + sys_offnom = generate_sys( + omega_n_11_offnom, + omega_n_22_offnom, + gamma_11_offnom, + gamma_22_offnom, + gain_12_offnom, + gain_21_offnom, + tau_12_offnom, + tau_21_offnom, + ) + sys_offnom_list.append(sys_offnom) + + # Frequency grid + omega_min = 0.05 + omega_max = 100 + num_omega = 100 + omega = np.logspace(np.log10(omega_min), np.log10(omega_max), num_omega) + + # Frequency response data of nominal and off-nominal systems + frequency_response_nom = control.frequency_response(sys_nom, omega) + frequency_response_offnom_list = control.frequency_response(sys_offnom_list, omega) + + # Complex response of nominal and off-nominal systems + complex_response_nom = frequency_response_nom.complex.transpose(2, 0, 1) + complex_response_offnom_list = [] + for frequency_response_offnom in frequency_response_offnom_list: + complex_response_offnom = frequency_response_offnom.complex.transpose(2, 0, 1) + complex_response_offnom_list.append(complex_response_offnom) + complex_response_offnom_list = np.array(complex_response_offnom_list) + + return { + "system_nominal": sys_nom, + "system_offnominal_list": sys_offnom_list, + "omega": omega, + "complex_response_nominal": complex_response_nom, + "complex_response_offnominal_list": complex_response_offnom_list, + } + + def _tf_close_coeff( tf_a: control.TransferFunction, tf_b: control.TransferFunction, @@ -495,3 +617,412 @@ def _auto_lmi_strictness( ) strictness = scale * tol return strictness + + +# TODO: Change variable names to be more descriptive and add more comments +def _fit_magnitude_log_chebyshev_siso( + omega: np.ndarray, + magnitude_fit: np.ndarray, + order: int, + magnitude_upper_bound: Optional[np.ndarray] = None, + magnitude_lower_bound: Optional[np.ndarray] = None, + weight: Optional[np.ndarray] = None, + linear_solver_params: Dict[str, Any] = {}, + tol_bisection: float = 1e-3, + max_iter_bisection: int = 500, + num_spec_constr: int = 500, +) -> control.StateSpace: + """ + Fit a stable and minimum-phase SISO LTI system to magnitude data using a + log-Chebyshev approach. + + Parameters + ---------- + omega : np.ndarray + Angular frequencies (rad/s). + magnitude_fit : np.ndarray + Magnitude response to fit the LTI system. + order : int + Order of the LTI system. + magnitude_upper_bound : Optional[np.ndarray] + Magnitude response for the upper bound constraint on the fitted LTI system + magnitude response. + magnitude_lower_bound : Optional[np.ndarray] + Magnitude response for the lower bound constraint on the fitted LTI system + magnitude response. + weight : np.ndarray + Frequency-dependent weight to encode bandwidths over which to prioritize the + accuracy of the LTI system fit. + linear_solver_param : Dict[str, Any] + Keyword arguments for the linear feasibility problem solver. + tol_bisection : float + Numerical tolerance for the bisection algorithm. + max_iter_bisection : int + Maximum allowable number of iterations in the bisection algorithm. + num_spec_constr : int + Number of constraints used to enforce the spectral factorizability of the + fitted autocorrelation. + + Returns + ------- + control.StateSpace + State-space system fitted to the magnitude data. + """ + + # Discrete-time frequency + theta, alpha = _convert_freq_halfplane_to_disk(omega) + + # Parse frequency-dependent weight + if weight is None: + weight = np.ones_like(omega) + + # Discrete-time autocorrelation + num_auto, den_auto = _fit_autocorrelation( + theta, + magnitude_fit, + order, + magnitude_upper_bound, + magnitude_lower_bound, + weight, + linear_solver_params, + tol_bisection, + max_iter_bisection, + num_spec_constr, + ) + + # Discrete-time transfer function + tf_dt = _compute_spectral_factorization(num_auto, den_auto) + + # Continous-time state-space + ss_ct = _convert_discrete_to_continuous_bilinear(tf_dt, alpha) + + return ss_ct + + +def _convert_freq_halfplane_to_disk(omega: np.ndarray) -> Tuple[np.ndarray, float]: + """ + Convert continous-time angular frequency data on the jw-axis in the s-domain into + discrete-time frequency data on the unit-circle of the z-domain via a bilinear + transformation. + + Parameters + ---------- + omega : np.ndarray + Continuous-time angular frequency data [rad/s]. + + Returns + ------- + Tuple[np.ndarray, float] + - Discrete-time frequency data, + - Equivalent sampling period [s]. + """ + + # Angular frequency bounds + omega_min = omega[0] + omega_max = omega[-1] + + # Bilinear transformation constant + alpha = np.sqrt(omega_min * omega_max) + + # Discrete-time frequency + theta = np.angle((-1j * omega - alpha) / (1j * omega - alpha)) + + return theta, alpha + + +# TODO: Update docstring +def _fit_autocorrelation( + theta: np.ndarray, + magnitude_fit: np.ndarray, + order: int, + magnitude_upper_bound: Optional[np.ndarray], + magnitude_lower_bound: Optional[np.ndarray], + weight: np.ndarray, + linear_solver_params: Dict[str, Any], + tol_bisection: float, + max_iter_bisection: int, + num_spec_constr: int, +) -> Tuple[np.ndarray, np.ndarray]: + """ + Fit a discrete-time biproper autocorrelation overbounding transfer function to + magnitude data. + + Parameters + ---------- + theta : np.ndarray + Discrete-time frequency. + magnitude_fit : np.ndarray + Magnitude response to fit the LTI system. + order : int + Order of the LTI system. + magnitude_upper_bound : Optional[np.ndarray] + Magnitude response for the upper bound constraint on the fitted LTI system + magnitude response. + magnitude_lower_bound : Optional[np.ndarray] + Magnitude response for the lower bound constraint on the fitted LTI system + magnitude response. + weight : np.ndarray + Frequency-dependent weight to encode bandwidths over which to prioritize the + accuracy of the LTI system fit. + linear_solver_params : Dict[str, Any] + Keyword arguments for the linear feasibility problem solver. + tol_bisection : float + Numerical tolerance for the bisection algorithm. + max_iter_bisection : int + Maximum allowable number of iterations in the bisection algorithm. + num_spec_constr : int + Number of constraints used to enforce the spectral factorizability of the + fitted autocorrelation. + + Returns + ------- + Tuple[np.ndarray, np.ndarray] + - Discrete-time autocorrelation transfer function denominator polynomial + coefficients. + - Discrete-time autocorrelation transfer function numerator polynomial + coefficients. + """ + + # Optimization problem parameters + psi_num = np.array( + [[np.cos(i * theta_idx) for i in range(order, -1, -1)] for theta_idx in theta] + ) + psi_den = np.array( + [[np.cos(i * theta_idx) for i in range(order, 0, -1)] for theta_idx in theta] + ) + magnitude_fit_sq_inv = np.diag(1 / magnitude_fit**2) + num_theta = theta.size + + # Spectral factorization constraint + A_spec, B_spec = _compute_spec_fac_constr_mat(num_spec_constr, order) + + # Lower bound constraint + if magnitude_lower_bound is not None: + magnitude_lower_bound_sq_inv = np.diag(1 / magnitude_lower_bound**2) + A_lower_bound = np.block([[-magnitude_lower_bound_sq_inv @ psi_num, psi_den]]) + B_lower_bound = -np.ones((num_theta, 1)) + + # Upper bound constraint + if magnitude_upper_bound is not None: + magnitude_upper_bound_sq_inv = np.diag(1 / magnitude_upper_bound**2) + A_upper_bound = np.block([[magnitude_upper_bound_sq_inv @ psi_num, -psi_den]]) + B_upper_bound = np.ones((num_theta, 1)) + + # Autocorrelation transfer function coefficient variables + num_auto_coef = cvxpy.Variable(shape=(order + 1, 1)) + den_auto_coef = cvxpy.Variable(shape=(order, 1)) + tf_auto_coef = cvxpy.bmat([[num_auto_coef], [den_auto_coef]]) + + # Error bounds for bisection algorithm + t_max = np.max(magnitude_fit) ** 2 / np.min(magnitude_fit) ** 2 - 1 + t_upper = t_max + t_lower = 0 + + # Initialize bisection algorithm + iter_bisection = 0 + feasibility_status = "infeasible" + + t = cvxpy.Parameter() + + # Bisection minimization of error upper bound + while np.abs(t_upper - t_lower) >= tol_bisection or feasibility_status != "optimal": + # Increment bisection iteration and stop if max number of iterations is reached. + iter_bisection += 1 + if iter_bisection >= max_iter_bisection: + break + + # Update weighted error bounds + t.value = 0.5 * (t_upper + t_lower) + gamma = cvxpy.diag(np.ones_like(weight) + t * 1 / weight) + gamma_inv = cvxpy.diag(1 / (np.ones_like(weight) + t * 1 / weight)) + + # Log-Chebyshev constraint + A_log_cheby_lower = cvxpy.bmat( + [[-magnitude_fit_sq_inv @ psi_num, gamma_inv @ psi_den]] + ) + B_log_cheby_lower = -gamma_inv @ np.ones((num_theta, 1)) + A_log_cheby_upper = cvxpy.bmat( + [[magnitude_fit_sq_inv @ psi_num, -gamma @ psi_den]] + ) + B_log_cheby_upper = gamma @ np.ones((num_theta, 1)) + + # Total constraint matrix + A_tot = cvxpy.bmat([[A_log_cheby_upper], [A_log_cheby_lower], [A_spec]]) + B_tot = cvxpy.bmat([[B_log_cheby_upper], [B_log_cheby_lower], [B_spec]]) + if magnitude_lower_bound is not None: + A_tot = cvxpy.bmat([[A_tot], [A_lower_bound]]) + B_tot = cvxpy.bmat([[B_tot], [B_lower_bound]]) + if magnitude_upper_bound is not None: + A_tot = cvxpy.bmat([[A_tot], [A_upper_bound]]) + B_tot = cvxpy.bmat([[B_tot], [B_upper_bound]]) + + # Linear feasibility problem + objective = cvxpy.Minimize(0) + constraint = [A_tot @ tf_auto_coef <= B_tot] + problem = cvxpy.Problem(objective, constraint) + try: + # Ignore warnings as some solvers issue warnings that are not applicable + # as we are expecting infeasible solutions in some cases due to the nature + # of the bisection algorithm + with warnings.catch_warnings(action="ignore"): + problem.solve(**linear_solver_params, ignore_dpp=True) + feasibility_status = problem.status + except cvxpy.SolverError: + t_lower = t.value + if feasibility_status == "optimal": + t_upper = t.value + else: + t_lower = t.value + + if num_auto_coef.value is None or den_auto_coef.value is None: + raise cvxpy.SolverError("The linear feasibility problem is infeasible.") + + # Autocorrelation polynomial coefficients + num_auto_coef_opt = num_auto_coef.value.reshape(-1) + den_auto_coef_opt = np.append(den_auto_coef.value.reshape(-1), 1) + + # Autocorrelation numerator and denominator polynominals + num_auto = np.concatenate( + ( + 0.5 * num_auto_coef_opt[:-1], + num_auto_coef_opt[[-1]], + np.flip(0.5 * num_auto_coef_opt[:-1]), + ) + ) + den_auto = np.concatenate( + ( + 0.5 * den_auto_coef_opt[:-1], + den_auto_coef_opt[[-1]], + np.flip(0.5 * den_auto_coef_opt[:-1]), + ) + ) + + return num_auto, den_auto + + +def _compute_spec_fac_constr_mat( + num_spec_constr: int, order: int +) -> Tuple[np.ndarray, np.ndarray]: + """ + Compute the matrices for the approximate linear spectral factorization constraint of + the form Ax <= B. This ensures that the autocorrelation transfer function can be + factored into the fitted transfer function. + + Parameters + ---------- + num_spec_constr : int + Number of discretization points for spectral factorization constraint + approximation. + order : int + Biproper transfer function order. + + Returns + ------- + Tuple[np.ndarray, np.ndarray] + - Spectral factorization constraint A matrix + - Spectral factorization constraint B matrix + """ + + # Discrete-time frequency grid approximation + theta_spec = np.linspace(0, np.pi, num_spec_constr) + + # Numerator and denominator regression matrices + psi_num_spec = np.array( + [ + [np.cos(i * theta_idx) for i in range(order, -1, -1)] + for theta_idx in theta_spec + ] + ) + psi_den_spec = np.array( + [ + [np.cos(i * theta_idx) for i in range(order, 0, -1)] + for theta_idx in theta_spec + ] + ) + + # Linear constraint matrices + A_spec = np.block( + [ + [-psi_num_spec, np.zeros((num_spec_constr, order))], + [np.zeros((num_spec_constr, order + 1)), -psi_den_spec], + ] + ) + B_spec = np.block( + [[np.zeros((num_spec_constr, 1))], [np.ones((num_spec_constr, 1))]] + ) + + return A_spec, B_spec + + +def _compute_spectral_factorization( + num_auto: np.ndarray, den_auto: np.ndarray +) -> control.TransferFunction: + # Autocorrelation numerator and denominator polynominal roots + roots_auto_num = np.roots(num_auto) + roots_auto_den = np.roots(den_auto) + + # Stable autocorrelation numerator and denominator polynominal roots + roots_num_stable = roots_auto_num[np.abs(roots_auto_num) < 1] + roots_den_stable = roots_auto_den[np.abs(roots_auto_den) < 1] + + # Spectral factorization + spectral_const_num = np.sqrt(num_auto[0] / np.prod(-roots_num_stable)) + spectral_const_den = np.sqrt(den_auto[0] / np.prod(-roots_den_stable)) + num = np.real( + spectral_const_num + * np.flip(np.polynomial.polynomial.polyfromroots(roots_num_stable)) + ) + den = np.real( + spectral_const_den + * np.flip(np.polynomial.polynomial.polyfromroots(roots_den_stable)) + ) + + # Transfer function + tf = control.TransferFunction(num, den) + + return tf + + +def _convert_discrete_to_continuous_bilinear( + sys_dt: Union[control.TransferFunction, control.StateSpace], + alpha: float, +) -> control.StateSpace: + """ + Convert a discrete-time LTI system (transfer function or state-space) into a + continous-time state-space system using the bilinear transformation. + + Parameters + ---------- + sys_dt : Union[control.TransferFunction, control.StateSpace] + Discrete-time system to be converted to continuous time. + alpha : float + Bilinear transformation constant. + + Returns + ------- + control.StateSpace + Continuous-time system. + """ + # Convert transfer function to state-space system + sys_dt = control.StateSpace(control.ss(sys_dt)) + + # Discrete-time state-space matrices + Ad = sys_dt.A + Bd = sys_dt.B + Cd = sys_dt.C + Dd = sys_dt.D + + # Additional matrices + In = np.eye(Ad.shape[0]) + iden_Ad_inv = np.linalg.solve(In + Ad, In) + + # Continous-time matrices + Ac = alpha * (Ad - In) @ iden_Ad_inv + Bc = alpha * (In - (Ad - In) @ iden_Ad_inv) @ Bd + Cc = Cd @ iden_Ad_inv + Dc = Dd - Cd @ iden_Ad_inv @ Bd + + # Continuous-time system + sys_dt = control.StateSpace(Ac, Bc, Cc, Dc) + + return sys_dt diff --git a/tests/test_uncertainty_characterization.py b/tests/test_uncertainty_characterization.py new file mode 100644 index 0000000..963eeb2 --- /dev/null +++ b/tests/test_uncertainty_characterization.py @@ -0,0 +1,1086 @@ +import control +import numpy as np +import pytest + +import dkpy + + +class TestComputeUncertaintyResidualResponse: + """Test :func:`compute_uncertainty_residual_response`.""" + + @pytest.mark.parametrize( + "sys_nom, sys_offnom_list, omega", + [ + ( + control.TransferFunction([1], [0.5, 1]), + [ + control.TransferFunction([1], [0.3, 1]), + control.TransferFunction([1], [0.4, 1]), + control.TransferFunction([1], [0.6, 1]), + control.TransferFunction([1], [0.7, 1]), + ], + np.logspace(-2, 2, 100), + ), + ( + control.TransferFunction([1], [1, 2 * 0.5 * 1, 1**2]), + [ + control.TransferFunction([1], [1, 2 * 0.3 * 1.5, 1.5**2]), + control.TransferFunction([1], [1, 2 * 0.7 * 1.3, 1.3**2]), + control.TransferFunction([1], [1, 2 * 0.2 * 0.9, 0.9**2]), + control.TransferFunction([1], [1, 2 * 0.9 * 0.7, 0.7**2]), + ], + np.logspace(-1, 1, 100), + ), + ( + control.TransferFunction( + [ + [[1], [3]], + [[2], [1]], + ], + [ + [[1, 2 * 0.5 * 1, 1**2], [0.5, 1]], + [[1, 1], [1, 2 * 0.3 * 5, 5**2]], + ], + ), + [ + control.TransferFunction( + [ + [[1], [3.25]], + [[1.9], [1]], + ], + [ + [[1, 2 * 0.4 * 1.1, 1.1**2], [0.6, 1]], + [[0.8, 1], [1, 2 * 0.4 * 4.2, 4.2**2]], + ], + ), + control.TransferFunction( + [ + [[1.1], [2.85]], + [[1.7], [0.9]], + ], + [ + [[1, 2 * 0.45 * 0.9, 0.9**2], [0.55, 1]], + [[0.95, 1], [1, 2 * 0.8 * 5.5, 5.5**2]], + ], + ), + control.TransferFunction( + [ + [[1], [3.05]], + [[1.7], [1.0]], + ], + [ + [[1, 2 * 0.5 * 1.3, 1.3**2], [0.42, 1]], + [[1, 1], [1, 2 * 0.9 * 5.25, 5.25**2]], + ], + ), + ], + np.logspace(-1, 1.5, 100), + ), + ], + ) + def test_compute_uncertainty_residual_response( + self, ndarrays_regression, sys_nom, sys_offnom_list, omega + ): + """Regression test :func:`compute_uncertainty_residual_response`.""" + + # Frequency response of systems + frequency_response_nom = control.frequency_response( + sys_nom, omega, squeeze=False + ) + frequency_response_offnom_list = control.frequency_response( + sys_offnom_list, omega, squeeze=False + ) + # Complex frequency response of systems + complex_response_nom = frequency_response_nom.complex.transpose(2, 0, 1) + complex_response_offnom_list = [] + for frequency_response_offnom in frequency_response_offnom_list: + complex_response_offnom_list.append( + frequency_response_offnom.complex.transpose(2, 0, 1) + ) + complex_response_offnom_list = np.array( + complex_response_offnom_list, dtype=complex + ) + + # Uncertainty residual computation + uncertainty_model_list = [ + "additive", + "multiplicative_input", + "multiplicative_output", + "inverse_additive", + "inverse_multiplicative_input", + "inverse_multiplicative_output", + ] + frequency_response_residual_dict = dkpy.compute_uncertainty_residual_response( + complex_response_nom, + complex_response_offnom_list, + uncertainty_model_list, + ) + ndarrays_regression.check( + frequency_response_residual_dict, + default_tolerance=dict(atol=1e-5, rtol=0), + ) + + @pytest.mark.parametrize( + "complex_response_nom_freq, complex_response_offnom_freq", + [ + ( + np.array([[1, 0], [0, 1], [1, 1]]), + np.array([[7, 8], [9, 10], [11, 12]]), + ), + ], + ) + def test_compute_uncertainty_residual_multiplicative_input_freq( + self, complex_response_nom_freq, complex_response_offnom_freq + ): + """Test ValueError of + :func:`_test_compute_uncertainty_residual_multiplicative_input_freq`. + """ + + with pytest.raises(ValueError): + residual_freq = dkpy.uncertainty_characterization._compute_uncertainty_residual_multiplicative_input_freq( + complex_response_nom_freq, complex_response_offnom_freq + ) + + @pytest.mark.parametrize( + "complex_response_nom_freq, complex_response_offnom_freq", + [ + ( + np.array([[1, 0], [0, 1], [1, 1]]).T, + np.array([[7, 8], [9, 10], [11, 12]]).T, + ), + ], + ) + def test_compute_uncertainty_residual_multiplicative_output_freq( + self, complex_response_nom_freq, complex_response_offnom_freq + ): + """Test ValueError of + :func:`_test_compute_uncertainty_residual_multiplicative_output_freq`. + """ + with pytest.raises(ValueError): + residual_freq = dkpy.uncertainty_characterization._compute_uncertainty_residual_multiplicative_output_freq( + complex_response_nom_freq, complex_response_offnom_freq + ) + + @pytest.mark.parametrize( + "complex_response_nom_freq, complex_response_offnom_freq", + [ + ( + np.array([[1, 0], [0, 1], [1, 1]]), + np.array([[7, 8], [9, 10], [11, 12]]), + ), + ( + np.array([[1, 0], [0, 1], [1, 1]]).T, + np.array([[7, 8], [9, 10], [11, 12]]).T, + ), + ], + ) + def test_compute_uncertainty_residual_inverse_additive_freq( + self, complex_response_nom_freq, complex_response_offnom_freq + ): + """Test ValueError of + :func:`_test_compute_uncertainty_residual_inverse_additive_freq`. + """ + with pytest.raises(ValueError): + residual_freq = dkpy.uncertainty_characterization._compute_uncertainty_residual_inverse_additive_freq( + complex_response_nom_freq, complex_response_offnom_freq + ) + + @pytest.mark.parametrize( + "complex_response_nom_freq, complex_response_offnom_freq", + [ + ( + np.array([[1, 0], [0, 1], [1, 1]]), + np.array([[7, 8], [9, 10], [11, 12]]), + ), + ], + ) + def test_compute_uncertainty_residual_inverse_multiplicative_input_freq( + self, complex_response_nom_freq, complex_response_offnom_freq + ): + """Test ValueError of + :func:`_test_compute_uncertainty_residual_inverse_multiplicative_input_freq`. + """ + with pytest.raises(ValueError): + residual_freq = dkpy.uncertainty_characterization._compute_uncertainty_residual_inverse_multiplicative_input_freq( + complex_response_nom_freq, complex_response_offnom_freq + ) + + @pytest.mark.parametrize( + "complex_response_nom_freq, complex_response_offnom_freq", + [ + ( + np.array([[1, 0], [0, 1], [1, 1]]).T, + np.array([[7, 8], [9, 10], [11, 12]]).T, + ), + ], + ) + def test_compute_uncertainty_residual_inverse_multiplicative_output_freq( + self, complex_response_nom_freq, complex_response_offnom_freq + ): + """Test ValueError of + :func:`_test_compute_uncertainty_residual_inverse_multiplicative_output_freq`. + """ + with pytest.raises(ValueError): + residual_freq = dkpy.uncertainty_characterization._compute_uncertainty_residual_inverse_multiplicative_output_freq( + complex_response_nom_freq, complex_response_offnom_freq + ) + + +class TestComputeOptimalUncertaintyWeightResponse: + """Test :func:`compute_optimal_uncertainty_weight_response`.""" + + @pytest.mark.parametrize( + "sys_nom, sys_offnom_list, omega, weight_left_structure, weight_right_structure", + [ + ( + control.TransferFunction([1], [0.5, 1]), + [ + control.TransferFunction([1], [0.3, 1]), + control.TransferFunction([1], [0.4, 1]), + control.TransferFunction([1], [0.6, 1]), + control.TransferFunction([1], [0.7, 1]), + ], + np.logspace(-2, 2, 100), + "diagonal", + "diagonal", + ), + ( + control.TransferFunction([1], [0.5, 1]), + [ + control.TransferFunction([1], [0.3, 1]), + control.TransferFunction([1], [0.4, 1]), + control.TransferFunction([1], [0.6, 1]), + control.TransferFunction([1], [0.7, 1]), + ], + np.logspace(-2, 2, 100), + "diagonal", + "identity", + ), + ( + control.TransferFunction([1], [0.5, 1]), + [ + control.TransferFunction([1], [0.3, 1]), + control.TransferFunction([1], [0.4, 1]), + control.TransferFunction([1], [0.6, 1]), + control.TransferFunction([1], [0.7, 1]), + ], + np.logspace(-2, 2, 100), + "identity", + "diagonal", + ), + ( + control.TransferFunction([1], [0.5, 1]), + [ + control.TransferFunction([1], [0.3, 1]), + control.TransferFunction([1], [0.4, 1]), + control.TransferFunction([1], [0.6, 1]), + control.TransferFunction([1], [0.7, 1]), + ], + np.logspace(-2, 2, 100), + "scalar", + "identity", + ), + ( + control.TransferFunction([1], [0.5, 1]), + [ + control.TransferFunction([1], [0.3, 1]), + control.TransferFunction([1], [0.4, 1]), + control.TransferFunction([1], [0.6, 1]), + control.TransferFunction([1], [0.7, 1]), + ], + np.logspace(-2, 2, 100), + "identity", + "scalar", + ), + ( + control.TransferFunction([1], [1, 2 * 0.5 * 1, 1**2]), + [ + control.TransferFunction([1], [1, 2 * 0.3 * 1.5, 1.5**2]), + control.TransferFunction([1], [1, 2 * 0.7 * 1.3, 1.3**2]), + control.TransferFunction([1], [1, 2 * 0.2 * 0.9, 0.9**2]), + control.TransferFunction([1], [1, 2 * 0.9 * 0.7, 0.7**2]), + ], + np.logspace(-1, 1, 100), + "diagonal", + "diagonal", + ), + ( + control.TransferFunction([1], [1, 2 * 0.5 * 1, 1**2]), + [ + control.TransferFunction([1], [1, 2 * 0.3 * 1.5, 1.5**2]), + control.TransferFunction([1], [1, 2 * 0.7 * 1.3, 1.3**2]), + control.TransferFunction([1], [1, 2 * 0.2 * 0.9, 0.9**2]), + control.TransferFunction([1], [1, 2 * 0.9 * 0.7, 0.7**2]), + ], + np.logspace(-1, 1, 100), + "diagonal", + "identity", + ), + ( + control.TransferFunction([1], [1, 2 * 0.5 * 1, 1**2]), + [ + control.TransferFunction([1], [1, 2 * 0.3 * 1.5, 1.5**2]), + control.TransferFunction([1], [1, 2 * 0.7 * 1.3, 1.3**2]), + control.TransferFunction([1], [1, 2 * 0.2 * 0.9, 0.9**2]), + control.TransferFunction([1], [1, 2 * 0.9 * 0.7, 0.7**2]), + ], + np.logspace(-1, 1, 100), + "identity", + "diagonal", + ), + ( + control.TransferFunction([1], [1, 2 * 0.5 * 1, 1**2]), + [ + control.TransferFunction([1], [1, 2 * 0.3 * 1.5, 1.5**2]), + control.TransferFunction([1], [1, 2 * 0.7 * 1.3, 1.3**2]), + control.TransferFunction([1], [1, 2 * 0.2 * 0.9, 0.9**2]), + control.TransferFunction([1], [1, 2 * 0.9 * 0.7, 0.7**2]), + ], + np.logspace(-1, 1, 100), + "scalar", + "identity", + ), + ( + control.TransferFunction([1], [1, 2 * 0.5 * 1, 1**2]), + [ + control.TransferFunction([1], [1, 2 * 0.3 * 1.5, 1.5**2]), + control.TransferFunction([1], [1, 2 * 0.7 * 1.3, 1.3**2]), + control.TransferFunction([1], [1, 2 * 0.2 * 0.9, 0.9**2]), + control.TransferFunction([1], [1, 2 * 0.9 * 0.7, 0.7**2]), + ], + np.logspace(-1, 1, 100), + "identity", + "scalar", + ), + ( + control.TransferFunction( + [ + [[1], [3]], + [[2], [1]], + ], + [ + [[1, 2 * 0.5 * 1, 1**2], [0.5, 1]], + [[1, 1], [1, 2 * 0.3 * 5, 5**2]], + ], + ), + [ + control.TransferFunction( + [ + [[1], [3.25]], + [[1.9], [1]], + ], + [ + [[1, 2 * 0.4 * 1.1, 1.1**2], [0.6, 1]], + [[0.8, 1], [1, 2 * 0.4 * 4.2, 4.2**2]], + ], + ), + control.TransferFunction( + [ + [[1.1], [2.85]], + [[1.7], [0.9]], + ], + [ + [[1, 2 * 0.45 * 0.9, 0.9**2], [0.55, 1]], + [[0.95, 1], [1, 2 * 0.8 * 5.5, 5.5**2]], + ], + ), + control.TransferFunction( + [ + [[1], [3.05]], + [[1.7], [1.0]], + ], + [ + [[1, 2 * 0.5 * 1.3, 1.3**2], [0.42, 1]], + [[1, 1], [1, 2 * 0.9 * 5.25, 5.25**2]], + ], + ), + ], + np.logspace(-1, 1.5, 100), + "diagonal", + "diagonal", + ), + ( + control.TransferFunction( + [ + [[1], [3]], + [[2], [1]], + ], + [ + [[1, 2 * 0.5 * 1, 1**2], [0.5, 1]], + [[1, 1], [1, 2 * 0.3 * 5, 5**2]], + ], + ), + [ + control.TransferFunction( + [ + [[1], [3.25]], + [[1.9], [1]], + ], + [ + [[1, 2 * 0.4 * 1.1, 1.1**2], [0.6, 1]], + [[0.8, 1], [1, 2 * 0.4 * 4.2, 4.2**2]], + ], + ), + control.TransferFunction( + [ + [[1.1], [2.85]], + [[1.7], [0.9]], + ], + [ + [[1, 2 * 0.45 * 0.9, 0.9**2], [0.55, 1]], + [[0.95, 1], [1, 2 * 0.8 * 5.5, 5.5**2]], + ], + ), + control.TransferFunction( + [ + [[1], [3.05]], + [[1.7], [1.0]], + ], + [ + [[1, 2 * 0.5 * 1.3, 1.3**2], [0.42, 1]], + [[1, 1], [1, 2 * 0.9 * 5.25, 5.25**2]], + ], + ), + ], + np.logspace(-1, 1.5, 100), + "diagonal", + "identity", + ), + ( + control.TransferFunction( + [ + [[1], [3]], + [[2], [1]], + ], + [ + [[1, 2 * 0.5 * 1, 1**2], [0.5, 1]], + [[1, 1], [1, 2 * 0.3 * 5, 5**2]], + ], + ), + [ + control.TransferFunction( + [ + [[1], [3.25]], + [[1.9], [1]], + ], + [ + [[1, 2 * 0.4 * 1.1, 1.1**2], [0.6, 1]], + [[0.8, 1], [1, 2 * 0.4 * 4.2, 4.2**2]], + ], + ), + control.TransferFunction( + [ + [[1.1], [2.85]], + [[1.7], [0.9]], + ], + [ + [[1, 2 * 0.45 * 0.9, 0.9**2], [0.55, 1]], + [[0.95, 1], [1, 2 * 0.8 * 5.5, 5.5**2]], + ], + ), + control.TransferFunction( + [ + [[1], [3.05]], + [[1.7], [1.0]], + ], + [ + [[1, 2 * 0.5 * 1.3, 1.3**2], [0.42, 1]], + [[1, 1], [1, 2 * 0.9 * 5.25, 5.25**2]], + ], + ), + ], + np.logspace(-1, 1.5, 100), + "identity", + "diagonal", + ), + ( + control.TransferFunction( + [ + [[1], [3]], + [[2], [1]], + ], + [ + [[1, 2 * 0.5 * 1, 1**2], [0.5, 1]], + [[1, 1], [1, 2 * 0.3 * 5, 5**2]], + ], + ), + [ + control.TransferFunction( + [ + [[1], [3.25]], + [[1.9], [1]], + ], + [ + [[1, 2 * 0.4 * 1.1, 1.1**2], [0.6, 1]], + [[0.8, 1], [1, 2 * 0.4 * 4.2, 4.2**2]], + ], + ), + control.TransferFunction( + [ + [[1.1], [2.85]], + [[1.7], [0.9]], + ], + [ + [[1, 2 * 0.45 * 0.9, 0.9**2], [0.55, 1]], + [[0.95, 1], [1, 2 * 0.8 * 5.5, 5.5**2]], + ], + ), + control.TransferFunction( + [ + [[1], [3.05]], + [[1.7], [1.0]], + ], + [ + [[1, 2 * 0.5 * 1.3, 1.3**2], [0.42, 1]], + [[1, 1], [1, 2 * 0.9 * 5.25, 5.25**2]], + ], + ), + ], + np.logspace(-1, 1.5, 100), + "scalar", + "identity", + ), + ( + control.TransferFunction( + [ + [[1], [3]], + [[2], [1]], + ], + [ + [[1, 2 * 0.5 * 1, 1**2], [0.5, 1]], + [[1, 1], [1, 2 * 0.3 * 5, 5**2]], + ], + ), + [ + control.TransferFunction( + [ + [[1], [3.25]], + [[1.9], [1]], + ], + [ + [[1, 2 * 0.4 * 1.1, 1.1**2], [0.6, 1]], + [[0.8, 1], [1, 2 * 0.4 * 4.2, 4.2**2]], + ], + ), + control.TransferFunction( + [ + [[1.1], [2.85]], + [[1.7], [0.9]], + ], + [ + [[1, 2 * 0.45 * 0.9, 0.9**2], [0.55, 1]], + [[0.95, 1], [1, 2 * 0.8 * 5.5, 5.5**2]], + ], + ), + control.TransferFunction( + [ + [[1], [3.05]], + [[1.7], [1.0]], + ], + [ + [[1, 2 * 0.5 * 1.3, 1.3**2], [0.42, 1]], + [[1, 1], [1, 2 * 0.9 * 5.25, 5.25**2]], + ], + ), + ], + np.logspace(-1, 1.5, 100), + "identity", + "scalar", + ), + ], + ) + def test_compute_optimal_uncertainty_weight_response( + self, + ndarrays_regression, + sys_nom, + sys_offnom_list, + omega, + weight_left_structure, + weight_right_structure, + ): + """Regression test :func:`compute_optimal_uncertainty_weight_response`.""" + # Frequency response of systems + frequency_response_nom = control.frequency_response( + sys_nom, omega, squeeze=False + ) + frequency_response_offnom_list = control.frequency_response( + sys_offnom_list, omega, squeeze=False + ) + # Complex frequency response of systems + complex_response_nom = frequency_response_nom.complex.transpose(2, 0, 1) + complex_response_offnom_list = [] + for frequency_response_offnom in frequency_response_offnom_list: + complex_response_offnom_list.append( + frequency_response_offnom.complex.transpose(2, 0, 1) + ) + complex_response_offnom_list = np.array( + complex_response_offnom_list, dtype=complex + ) + + # Uncertainty residual computation + uncertainty_model_list = ["multiplicative_input"] + complex_response_residual_dict = dkpy.compute_uncertainty_residual_response( + complex_response_nom, + complex_response_offnom_list, + uncertainty_model_list, + ) + + # Optimal uncertainty weight + complex_response_weight_left, complex_response_weight_right = ( + dkpy.compute_optimal_uncertainty_weight_response( + complex_response_residual_dict["multiplicative_input"], + weight_left_structure, + weight_right_structure, + ) + ) + complex_response_weight_dict = { + "left": complex_response_weight_left, + "right": complex_response_weight_right, + } + + # Regression testing + ndarrays_regression.check( + complex_response_weight_dict, + default_tolerance=dict(atol=1e-5, rtol=0), + ) + + +class TestFitOverboundingUncertaintyWeight: + """Test :func:`fit_overbounding_uncertainty_weight`.""" + + @pytest.mark.parametrize( + "sys_nom, sys_offnom_list, omega, weight_left_structure, weight_right_structure, fit_order", + [ + ( + control.TransferFunction([1], [0.5, 1]), + [ + control.TransferFunction([1], [0.3, 1]), + control.TransferFunction([1], [0.4, 1]), + control.TransferFunction([1], [0.6, 1]), + control.TransferFunction([1], [0.7, 1]), + ], + np.logspace(-2, 2, 100), + "diagonal", + "diagonal", + 4, + ), + ( + control.TransferFunction([1], [1, 2 * 0.5 * 1, 1**2]), + [ + control.TransferFunction([1], [1, 2 * 0.3 * 1.5, 1.5**2]), + control.TransferFunction([1], [1, 2 * 0.7 * 1.3, 1.3**2]), + control.TransferFunction([1], [1, 2 * 0.2 * 0.9, 0.9**2]), + control.TransferFunction([1], [1, 2 * 0.9 * 0.7, 0.7**2]), + ], + np.logspace(-1, 1, 100), + "diagonal", + "diagonal", + 4, + ), + ( + control.TransferFunction( + [ + [[1], [3]], + [[2], [1]], + ], + [ + [[1, 2 * 0.5 * 1, 1**2], [0.5, 1]], + [[1, 1], [1, 2 * 0.3 * 5, 5**2]], + ], + ), + [ + control.TransferFunction( + [ + [[1], [3.25]], + [[1.9], [1]], + ], + [ + [[1, 2 * 0.4 * 1.1, 1.1**2], [0.6, 1]], + [[0.8, 1], [1, 2 * 0.4 * 4.2, 4.2**2]], + ], + ), + control.TransferFunction( + [ + [[1.1], [2.85]], + [[1.7], [0.9]], + ], + [ + [[1, 2 * 0.45 * 0.9, 0.9**2], [0.55, 1]], + [[0.95, 1], [1, 2 * 0.8 * 5.5, 5.5**2]], + ], + ), + control.TransferFunction( + [ + [[1], [3.05]], + [[1.7], [1.0]], + ], + [ + [[1, 2 * 0.5 * 1.3, 1.3**2], [0.42, 1]], + [[1, 1], [1, 2 * 0.9 * 5.25, 5.25**2]], + ], + ), + ], + np.logspace(-1, 1.5, 100), + "diagonal", + "diagonal", + 4, + ), + ( + control.TransferFunction( + [ + [[1], [3]], + [[2], [1]], + ], + [ + [[1, 2 * 0.5 * 1, 1**2], [0.5, 1]], + [[1, 1], [1, 2 * 0.3 * 5, 5**2]], + ], + ), + [ + control.TransferFunction( + [ + [[1], [3.25]], + [[1.9], [1]], + ], + [ + [[1, 2 * 0.4 * 1.1, 1.1**2], [0.6, 1]], + [[0.8, 1], [1, 2 * 0.4 * 4.2, 4.2**2]], + ], + ), + control.TransferFunction( + [ + [[1.1], [2.85]], + [[1.7], [0.9]], + ], + [ + [[1, 2 * 0.45 * 0.9, 0.9**2], [0.55, 1]], + [[0.95, 1], [1, 2 * 0.8 * 5.5, 5.5**2]], + ], + ), + control.TransferFunction( + [ + [[1], [3.05]], + [[1.7], [1.0]], + ], + [ + [[1, 2 * 0.5 * 1.3, 1.3**2], [0.42, 1]], + [[1, 1], [1, 2 * 0.9 * 5.25, 5.25**2]], + ], + ), + ], + np.logspace(-1, 1.5, 100), + "diagonal", + "diagonal", + [4, 4], + ), + ( + control.TransferFunction( + [ + [[1], [3]], + [[2], [1]], + ], + [ + [[1, 2 * 0.5 * 1, 1**2], [0.5, 1]], + [[1, 1], [1, 2 * 0.3 * 5, 5**2]], + ], + ), + [ + control.TransferFunction( + [ + [[1], [3.25]], + [[1.9], [1]], + ], + [ + [[1, 2 * 0.4 * 1.1, 1.1**2], [0.6, 1]], + [[0.8, 1], [1, 2 * 0.4 * 4.2, 4.2**2]], + ], + ), + control.TransferFunction( + [ + [[1.1], [2.85]], + [[1.7], [0.9]], + ], + [ + [[1, 2 * 0.45 * 0.9, 0.9**2], [0.55, 1]], + [[0.95, 1], [1, 2 * 0.8 * 5.5, 5.5**2]], + ], + ), + control.TransferFunction( + [ + [[1], [3.05]], + [[1.7], [1.0]], + ], + [ + [[1, 2 * 0.5 * 1.3, 1.3**2], [0.42, 1]], + [[1, 1], [1, 2 * 0.9 * 5.25, 5.25**2]], + ], + ), + ], + np.logspace(-1, 1.5, 100), + "diagonal", + "diagonal", + np.array([4.0, 4.0]), + ), + ], + ) + def test_fit_overbounding_uncertainty_weight( + self, + ndarrays_regression, + sys_nom, + sys_offnom_list, + omega, + weight_left_structure, + weight_right_structure, + fit_order, + ): + """Regression test :func:`test_fit_overbounding_uncertainty_weight`.""" + # Frequency response of systems + frequency_response_nom = control.frequency_response( + sys_nom, omega, squeeze=False + ) + frequency_response_offnom_list = control.frequency_response( + sys_offnom_list, omega, squeeze=False + ) + # Complex frequency response of systems + complex_response_nom = frequency_response_nom.complex.transpose(2, 0, 1) + complex_response_offnom_list = [] + for frequency_response_offnom in frequency_response_offnom_list: + complex_response_offnom_list.append( + frequency_response_offnom.complex.transpose(2, 0, 1) + ) + complex_response_offnom_list = np.array( + complex_response_offnom_list, dtype=complex + ) + + # Uncertainty residual computation + uncertainty_model_list = ["multiplicative_input"] + complex_response_residual_dict = dkpy.compute_uncertainty_residual_response( + complex_response_nom, + complex_response_offnom_list, + uncertainty_model_list, + ) + + # Optimal uncertainty weight + complex_response_weight_left, complex_response_weight_right = ( + dkpy.compute_optimal_uncertainty_weight_response( + complex_response_residual_dict["multiplicative_input"], + weight_left_structure, + weight_right_structure, + ) + ) + + # Overbounding transfer function fit + weight_left_fit = dkpy.fit_overbounding_uncertainty_weight( + complex_response_weight_left, omega, fit_order + ) + weight_right_fit = dkpy.fit_overbounding_uncertainty_weight( + complex_response_weight_right, omega, fit_order + ) + + # Uncertainty weight fit frequency response + frequency_response_weight_left_fit = control.frequency_response( + weight_left_fit, omega, squeeze=False + ) + frequency_response_weight_right_fit = control.frequency_response( + weight_right_fit, omega, squeeze=False + ) + complex_response_weight_left_fit = ( + frequency_response_weight_left_fit.complex.transpose(2, 0, 1) + ) + complex_response_weight_right_fit = ( + frequency_response_weight_right_fit.complex.transpose(2, 0, 1) + ) + complex_response_weight_fit_dict = { + "left": complex_response_weight_left_fit, + "right": complex_response_weight_right_fit, + } + + # Regression testing + ndarrays_regression.check( + complex_response_weight_fit_dict, + default_tolerance=dict(atol=1e-5, rtol=0), + ) + + +class TestConvertFrequencyResponseDataToArray: + @pytest.mark.parametrize( + "sys, omega", + [ + ( + control.TransferFunction([1], [0.5, 1]), + np.logspace(-2, 2, 100), + ), + ( + control.TransferFunction([1], [1, 2 * 0.5 * 1, 1**2]), + np.logspace(-1, 1, 100), + ), + ( + control.TransferFunction( + [ + [[1], [3]], + [[2], [1]], + ], + [ + [[1, 2 * 0.5 * 1, 1**2], [0.5, 1]], + [[1, 1], [1, 2 * 0.3 * 5, 5**2]], + ], + ), + np.logspace(-1, 1.5, 100), + ), + ], + ) + def test_convert_frequency_response_data_to_array(self, sys, omega): + # Desired frequency response array shape + complex_response_shape_true = (omega.size, sys.noutputs, sys.ninputs) + + # Evaluate frequency response + frequency_response = control.frequency_response(sys, omega) + + # Convert frequency response data + complex_response = ( + dkpy.uncertainty_characterization._convert_frequency_response_data_to_array( + frequency_response + ) + ) + + # Verify type + assert isinstance(complex_response, np.ndarray) + # Verify dimensions + assert complex_response.shape == complex_response_shape_true + # Verify complex datatype + assert np.all(np.iscomplex(complex_response)) + # Verify that the function does not affect responses in the correct format + assert np.all( + np.isclose( + complex_response, + dkpy.uncertainty_characterization._convert_frequency_response_data_to_array( + complex_response + ), + ), + ) + + @pytest.mark.parametrize( + "complex_response", + [ + np.ones(100, dtype=complex), + np.ones((100, 5), dtype=complex), + np.ones((100, 5, 5, 5), dtype=complex), + ], + ) + def test_convert_frequency_response_data_to_array_error(self, complex_response): + with pytest.raises(ValueError): + complex_response = dkpy.uncertainty_characterization._convert_frequency_response_data_to_array( + complex_response + ) + + +class TestConvertFrequencyResponseListToArray: + @pytest.mark.parametrize( + "sys_list, omega", + [ + ( + [ + control.TransferFunction([1], [0.3, 1]), + control.TransferFunction([1], [0.4, 1]), + control.TransferFunction([1], [0.6, 1]), + control.TransferFunction([1], [0.7, 1]), + ], + np.logspace(-2, 2, 100), + ), + ( + [ + control.TransferFunction([1], [1, 2 * 0.3 * 1.5, 1.5**2]), + control.TransferFunction([1], [1, 2 * 0.7 * 1.3, 1.3**2]), + control.TransferFunction([1], [1, 2 * 0.2 * 0.9, 0.9**2]), + control.TransferFunction([1], [1, 2 * 0.9 * 0.7, 0.7**2]), + ], + np.logspace(-1, 1, 100), + ), + ( + [ + control.TransferFunction( + [ + [[1], [3.25]], + [[1.9], [1]], + ], + [ + [[1, 2 * 0.4 * 1.1, 1.1**2], [0.6, 1]], + [[0.8, 1], [1, 2 * 0.4 * 4.2, 4.2**2]], + ], + ), + control.TransferFunction( + [ + [[1.1], [2.85]], + [[1.7], [0.9]], + ], + [ + [[1, 2 * 0.45 * 0.9, 0.9**2], [0.55, 1]], + [[0.95, 1], [1, 2 * 0.8 * 5.5, 5.5**2]], + ], + ), + control.TransferFunction( + [ + [[1], [3.05]], + [[1.7], [1.0]], + ], + [ + [[1, 2 * 0.5 * 1.3, 1.3**2], [0.42, 1]], + [[1, 1], [1, 2 * 0.9 * 5.25, 5.25**2]], + ], + ), + ], + np.logspace(-1, 1.5, 100), + ), + ], + ) + def test_convert_frequency_response_list_to_array(self, sys_list, omega): + # Desired frequency response array shape + complex_response_list_shape_true = ( + len(sys_list), + omega.size, + sys_list[0].noutputs, + sys_list[0].ninputs, + ) + + # Evaluate frequency response + frequency_response_list = control.frequency_response(sys_list, omega) + + # Convert frequency response data + complex_response_list = ( + dkpy.uncertainty_characterization._convert_frequency_response_list_to_array( + frequency_response_list + ) + ) + + # Verify type + assert isinstance(complex_response_list, np.ndarray) + # Verify dimensions + assert complex_response_list.shape == complex_response_list_shape_true + # Verify complex datatype + assert np.all(np.iscomplex(complex_response_list)) + # Verify that the function does not affect responses in the correct format + assert np.all( + np.isclose( + complex_response_list, + dkpy.uncertainty_characterization._convert_frequency_response_list_to_array( + complex_response_list + ), + ), + ) + + @pytest.mark.parametrize( + "complex_response_list", + [ + np.ones(100, dtype=complex), + np.ones((100, 5), dtype=complex), + np.ones((100, 5, 5), dtype=complex), + np.ones((100, 5, 5, 5, 5), dtype=complex), + ], + ) + def test_convert_frequency_response_list_to_array_error( + self, complex_response_list + ): + with pytest.raises(ValueError): + complex_response_list = dkpy.uncertainty_characterization._convert_frequency_response_list_to_array( + complex_response_list + ) diff --git a/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom0_sys_offnom_list0_omega0_diagonal_diagonal_.npz b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom0_sys_offnom_list0_omega0_diagonal_diagonal_.npz new file mode 100644 index 0000000..c0f49e7 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom0_sys_offnom_list0_omega0_diagonal_diagonal_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom10_sys_offnom_list10_omega10_diagonal_diagonal_.npz b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom10_sys_offnom_list10_omega10_diagonal_diagonal_.npz new file mode 100644 index 0000000..bf92a3d Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom10_sys_offnom_list10_omega10_diagonal_diagonal_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom11_sys_offnom_list11_omega11_diagonal_identity_.npz b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom11_sys_offnom_list11_omega11_diagonal_identity_.npz new file mode 100644 index 0000000..1df0f9f Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom11_sys_offnom_list11_omega11_diagonal_identity_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom12_sys_offnom_list12_omega12_identity_diagonal_.npz b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom12_sys_offnom_list12_omega12_identity_diagonal_.npz new file mode 100644 index 0000000..5c3c03c Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom12_sys_offnom_list12_omega12_identity_diagonal_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom13_sys_offnom_list13_omega13_scalar_identity_.npz b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom13_sys_offnom_list13_omega13_scalar_identity_.npz new file mode 100644 index 0000000..bed5220 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom13_sys_offnom_list13_omega13_scalar_identity_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom14_sys_offnom_list14_omega14_identity_scalar_.npz b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom14_sys_offnom_list14_omega14_identity_scalar_.npz new file mode 100644 index 0000000..848d598 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom14_sys_offnom_list14_omega14_identity_scalar_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom1_sys_offnom_list1_omega1_diagonal_identity_.npz b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom1_sys_offnom_list1_omega1_diagonal_identity_.npz new file mode 100644 index 0000000..7edb112 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom1_sys_offnom_list1_omega1_diagonal_identity_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom2_sys_offnom_list2_omega2_identity_diagonal_.npz b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom2_sys_offnom_list2_omega2_identity_diagonal_.npz new file mode 100644 index 0000000..a0a8ed1 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom2_sys_offnom_list2_omega2_identity_diagonal_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom3_sys_offnom_list3_omega3_scalar_identity_.npz b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom3_sys_offnom_list3_omega3_scalar_identity_.npz new file mode 100644 index 0000000..7edb112 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom3_sys_offnom_list3_omega3_scalar_identity_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom4_sys_offnom_list4_omega4_identity_scalar_.npz b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom4_sys_offnom_list4_omega4_identity_scalar_.npz new file mode 100644 index 0000000..a0a8ed1 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom4_sys_offnom_list4_omega4_identity_scalar_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom5_sys_offnom_list5_omega5_diagonal_diagonal_.npz b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom5_sys_offnom_list5_omega5_diagonal_diagonal_.npz new file mode 100644 index 0000000..2a7446e Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom5_sys_offnom_list5_omega5_diagonal_diagonal_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom6_sys_offnom_list6_omega6_diagonal_identity_.npz b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom6_sys_offnom_list6_omega6_diagonal_identity_.npz new file mode 100644 index 0000000..e7bf08e Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom6_sys_offnom_list6_omega6_diagonal_identity_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom7_sys_offnom_list7_omega7_identity_diagonal_.npz b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom7_sys_offnom_list7_omega7_identity_diagonal_.npz new file mode 100644 index 0000000..75abbdd Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom7_sys_offnom_list7_omega7_identity_diagonal_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom8_sys_offnom_list8_omega8_scalar_identity_.npz b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom8_sys_offnom_list8_omega8_scalar_identity_.npz new file mode 100644 index 0000000..e7bf08e Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom8_sys_offnom_list8_omega8_scalar_identity_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom9_sys_offnom_list9_omega9_identity_scalar_.npz b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom9_sys_offnom_list9_omega9_identity_scalar_.npz new file mode 100644 index 0000000..75abbdd Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_optimal_uncertainty_weight_response_sys_nom9_sys_offnom_list9_omega9_identity_scalar_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_uncertainty_residual_response_sys_nom0_sys_offnom_list0_omega0_.npz b/tests/test_uncertainty_characterization/test_compute_uncertainty_residual_response_sys_nom0_sys_offnom_list0_omega0_.npz new file mode 100644 index 0000000..8dd1770 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_uncertainty_residual_response_sys_nom0_sys_offnom_list0_omega0_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_uncertainty_residual_response_sys_nom1_sys_offnom_list1_omega1_.npz b/tests/test_uncertainty_characterization/test_compute_uncertainty_residual_response_sys_nom1_sys_offnom_list1_omega1_.npz new file mode 100644 index 0000000..99531e3 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_uncertainty_residual_response_sys_nom1_sys_offnom_list1_omega1_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_uncertainty_residual_response_sys_nom2_sys_offnom_list2_omega2_.npz b/tests/test_uncertainty_characterization/test_compute_uncertainty_residual_response_sys_nom2_sys_offnom_list2_omega2_.npz new file mode 100644 index 0000000..2211118 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_uncertainty_residual_response_sys_nom2_sys_offnom_list2_omega2_.npz differ diff --git a/tests/test_uncertainty_characterization/test_fit_overbounding_uncertainty_weight_sys_nom0_sys_offnom_list0_omega0_diagonal_diagonal_4_.npz b/tests/test_uncertainty_characterization/test_fit_overbounding_uncertainty_weight_sys_nom0_sys_offnom_list0_omega0_diagonal_diagonal_4_.npz new file mode 100644 index 0000000..8164121 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_fit_overbounding_uncertainty_weight_sys_nom0_sys_offnom_list0_omega0_diagonal_diagonal_4_.npz differ diff --git a/tests/test_uncertainty_characterization/test_fit_overbounding_uncertainty_weight_sys_nom1_sys_offnom_list1_omega1_diagonal_diagonal_4_.npz b/tests/test_uncertainty_characterization/test_fit_overbounding_uncertainty_weight_sys_nom1_sys_offnom_list1_omega1_diagonal_diagonal_4_.npz new file mode 100644 index 0000000..5b40d52 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_fit_overbounding_uncertainty_weight_sys_nom1_sys_offnom_list1_omega1_diagonal_diagonal_4_.npz differ diff --git a/tests/test_uncertainty_characterization/test_fit_overbounding_uncertainty_weight_sys_nom2_sys_offnom_list2_omega2_diagonal_diagonal_4_.npz b/tests/test_uncertainty_characterization/test_fit_overbounding_uncertainty_weight_sys_nom2_sys_offnom_list2_omega2_diagonal_diagonal_4_.npz new file mode 100644 index 0000000..b30ec1e Binary files /dev/null and b/tests/test_uncertainty_characterization/test_fit_overbounding_uncertainty_weight_sys_nom2_sys_offnom_list2_omega2_diagonal_diagonal_4_.npz differ diff --git a/tests/test_uncertainty_characterization/test_fit_overbounding_uncertainty_weight_sys_nom3_sys_offnom_list3_omega3_diagonal_diagonal_4_.npz b/tests/test_uncertainty_characterization/test_fit_overbounding_uncertainty_weight_sys_nom3_sys_offnom_list3_omega3_diagonal_diagonal_4_.npz new file mode 100644 index 0000000..9d67de8 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_fit_overbounding_uncertainty_weight_sys_nom3_sys_offnom_list3_omega3_diagonal_diagonal_4_.npz differ diff --git a/tests/test_uncertainty_characterization/test_fit_overbounding_uncertainty_weight_sys_nom3_sys_offnom_list3_omega3_diagonal_diagonal_fit_order3_.npz b/tests/test_uncertainty_characterization/test_fit_overbounding_uncertainty_weight_sys_nom3_sys_offnom_list3_omega3_diagonal_diagonal_fit_order3_.npz new file mode 100644 index 0000000..b30ec1e Binary files /dev/null and b/tests/test_uncertainty_characterization/test_fit_overbounding_uncertainty_weight_sys_nom3_sys_offnom_list3_omega3_diagonal_diagonal_fit_order3_.npz differ diff --git a/tests/test_uncertainty_characterization/test_fit_overbounding_uncertainty_weight_sys_nom4_sys_offnom_list4_omega4_diagonal_diagonal_fit_order4_.npz b/tests/test_uncertainty_characterization/test_fit_overbounding_uncertainty_weight_sys_nom4_sys_offnom_list4_omega4_diagonal_diagonal_fit_order4_.npz new file mode 100644 index 0000000..b30ec1e Binary files /dev/null and b/tests/test_uncertainty_characterization/test_fit_overbounding_uncertainty_weight_sys_nom4_sys_offnom_list4_omega4_diagonal_diagonal_fit_order4_.npz differ