Implement multi-model uncertainty quantification tools

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"],
+    )

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

examples/7_uncertainty_characterization.py

@@ -0,0 +1,88 @@
+"""Multi-model uncertainty characterization from frequency response data."""
+
+import numpy as np
+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()
+    complex_response_nom = eg["complex_response_nominal"]
+    complex_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_nom_offnom(
+        complex_response_nom,
+        complex_response_offnom_list,
+        omega,
+    )
+
+    # Plot: Phase response of nominal and off-nominal systems
+    fig, _ = dkpy.plot_phase_response_nom_offnom(
+        complex_response_nom,
+        complex_response_offnom_list,
+        omega,
+    )
+
+    # Plot: Singular value response of nominal and off-nominal systems
+    fig, ax = dkpy.plot_singular_value_response_nom_offnom(
+        complex_response_nom,
+        complex_response_offnom_list,
+        omega,
+    )
+
+    # Uncertainty models
+    uncertainty_models = {"A", "I", "O", "iA", "iI", "iO"}
+    complex_response_residuals_dict = dkpy.compute_uncertainty_residual_response(
+        complex_response_nom,
+        complex_response_offnom_list,
+        uncertainty_models,
+    )
+
+    # Plot: Singular value response of uncerainty residuals
+    figure_dict = dkpy.plot_singular_value_response_uncertainty_residual(
+        complex_response_residuals_dict, omega
+    )
+
+    # Plot: Comparison of singular value response of uncerainty residuals for each
+    # uncertainty model
+    fig, _ = dkpy.plot_singular_value_response_uncertainty_residual_comparison(
+        complex_response_residuals_dict, omega
+    )
+    fig.savefig("doc/_static/example_7/7_sval_max_residual.png")
+
+    # Compute uncertainty weight frequency response
+    complex_response_weight_left, complex_response_weight_right = (
+        dkpy.compute_optimal_uncertainty_weight_response(
+            complex_response_residuals_dict["iA"], "diagonal", "diagonal"
+        )
+    )
+
+    print("fit weight")
+    # Fit overbounding stable and minimum-phase uncertainty weight system
+    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]
+    )
+
+    # Plot: Magnitude response of uncertainty weight and overbounding uncertainty weight
+    # fit
+    fig, _ = dkpy.plot_magnitude_response_uncertainty_weight(
+        complex_response_weight_left,
+        complex_response_weight_right,
+        omega,
+        weight_left,
+        weight_right,
+    )
+    plt.show()
+
+
+if __name__ == "__main__":
+    example_uncertainty_characterization()

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 *