Add ability to specify unstructured/full uncertainty weight in optimal uncertainty weight response.

Author: timeverettadams

src/dkpy/uncertainty_characterization.py

@@ -20,7 +20,7 @@
 import scipy
 from matplotlib import pyplot as plt
 
-from typing import List, Optional, Union, Tuple, Dict, Callable, Set, Any
+from typing import List, Optional, Union, Tuple, Dict, Callable, Set, Any, Literal
 from matplotlib.figure import Figure
 from matplotlib.axes import Axes
 from matplotlib.legend import Legend
@@ -551,8 +551,8 @@ def _compute_uncertainty_residual_inverse_multiplicative_output_freq(
 
 def compute_uncertainty_weight_response(
     complex_response_residual_list: Union[np.ndarray, control.FrequencyResponseList],
-    weight_left_structure: str,
-    weight_right_structure: str,
+    weight_left_structure: Literal["full", "diagonal", "scalar", "identity"],
+    weight_right_structure: Literal["full", "diagonal", "scalar", "identity"],
     solver_params: Optional[Dict[str, Any]] = None,
 ) -> Tuple[np.ndarray, np.ndarray]:
     """Compute the optimal uncertainty weight frequency response.
@@ -566,11 +566,9 @@ def compute_uncertainty_weight_response(
         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".
+        Structure of the left uncertainty weight.
     weight_right_structure : str
-        Structure of the right uncertainty weight. Valid structures include: "diagonal",
-        "scalar", and "identity".
+        Structure of the right uncertainty weight.
     solver_param : Dict[str, Any]
         Keyword arguments for the convex optimization solver. See [#cvxpy_solver]_ for
         more information.
@@ -650,8 +648,8 @@ def compute_uncertainty_weight_response(
 
 def _compute_optimal_weight_freq(
     complex_residual_offnom_set_freq: np.ndarray,
-    weight_left_structure: str,
-    weight_right_structure: str,
+    weight_left_structure: Literal["full", "diagonal", "scalar", "identity"],
+    weight_right_structure: Literal["full", "diagonal", "scalar", "identity"],
     solver_params: Optional[Dict[str, Any]] = None,
 ) -> Tuple[np.ndarray, np.ndarray]:
     """Compute the optimal uncertainty weight at a given frequency.
@@ -661,14 +659,12 @@ def _compute_optimal_weight_freq(
     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".
+        Structure of the left uncertainty weight.
     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.
+        Structure of the right uncertainty weight.
+    solver_params : Dict[str, Any]
+        Solver parameters for the optimization problem. These are keyword arguments for
+        `cvxpy.Problem.solve()` [#cvxpy_solve]_.
 
     Returns
     -------
@@ -696,49 +692,51 @@ def _compute_optimal_weight_freq(
     )
 
     # 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]
+    nbr_left = complex_residual_offnom_set_freq.shape[1]
+    nbr_right = complex_residual_offnom_set_freq.shape[2]
+    nbr_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)
+    if weight_left_structure == "full":
+        L = cvxpy.Variable((nbr_left, nbr_left))
+    elif weight_left_structure == "diagonal":
+        L = cvxpy.Variable((nbr_left, nbr_left), diag=True)
     elif weight_left_structure == "scalar":
         L_scalar = cvxpy.Variable()
-        L = L_scalar * scipy.sparse.eye_array(num_left)
+        L = L_scalar * scipy.sparse.eye_array(nbr_left)
     elif weight_left_structure == "identity":
         L = cvxpy.Parameter(
-            shape=(num_left, num_left), value=np.eye(num_left), diag=True
+            shape=(nbr_left, nbr_left), value=np.eye(nbr_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".'
+            "Invalid `weight_left_structure`. Expected `full`, `diagonal`, `scalar`, "
+            f"or `identity`. (Got {weight_left_structure})."
         )
 
     # Generate right weight variable
-    if weight_right_structure == "diagonal":
-        R = cvxpy.Variable((num_right, num_right), diag=True)
+    if weight_right_structure == "full":
+        R = cvxpy.Variable((nbr_right, nbr_right))
+    elif weight_right_structure == "diagonal":
+        R = cvxpy.Variable((nbr_right, nbr_right), diag=True)
     elif weight_right_structure == "scalar":
         R_scalar = cvxpy.Variable()
-        R = R_scalar * scipy.sparse.eye_array(num_right)
+        R = R_scalar * scipy.sparse.eye_array(nbr_right)
     elif weight_right_structure == "identity":
         R = cvxpy.Parameter(
-            shape=(num_right, num_right), value=np.eye(num_right), diag=True
+            shape=(nbr_right, nbr_right), value=np.eye(nbr_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".'
+            "Invalid `weight_right_structure`. Expected `full`, `diagonal`, `scalar`, "
+            f"or `identity`. (Got {weight_right_structure})."
         )
 
     # Generate optimal uncertainty weight constraints
     constraint_freq_list = []
-    for idx_offnom in range(num_offnom):
+    for idx_offnom in range(nbr_offnom):
         E_k = cvxpy.Parameter(
-            shape=(num_left, num_right),
+            shape=(nbr_left, nbr_right),
             value=complex_residual_offnom_set_freq[idx_offnom, :, :],
             complex=True,
         )
@@ -761,19 +759,20 @@ def _compute_optimal_weight_freq(
     problem.solve(**solver_params)
 
     # Extract left weight
-    if weight_left_structure == "identity":
+    if weight_left_structure == "identity" or weight_left_structure == "full":
         L_value = np.array(L.value)
     else:
         L_value = np.array(L.value.toarray())
-    complex_response_weight_left_freq = np.sqrt(L_value)
+    weight_left_opt = np.linalg.cholesky(L_value)
+
     # Extract right weight
-    if weight_right_structure == "identity":
+    if weight_right_structure == "identity" or weight_right_structure == "full":
         R_value = np.array(R.value)
     else:
         R_value = np.array(R.value.toarray())
-    complex_response_weight_right_freq = np.sqrt(R_value)
+    weight_right_opt = np.linalg.cholesky(R_value)
 
-    return complex_response_weight_left_freq, complex_response_weight_right_freq
+    return weight_left_opt, weight_right_opt
 
 
 def fit_uncertainty_weight(

tests/test_uncertainty_characterization.py

@@ -241,8 +241,8 @@ class TestComputeUncertaintyWeightResponse:
                     control.TransferFunction([1], [0.7, 1]),
                 ],
                 np.logspace(-2, 2, 100),
-                "diagonal",
-                "diagonal",
+                "full",
+                "full",
             ),
             (
                 control.TransferFunction([1], [0.5, 1]),
@@ -254,18 +254,6 @@ class TestComputeUncertaintyWeightResponse:
                 ],
                 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",
             ),
             (
@@ -301,8 +289,8 @@ class TestComputeUncertaintyWeightResponse:
                     control.TransferFunction([1], [1, 2 * 0.9 * 0.7, 0.7**2]),
                 ],
                 np.logspace(-1, 1, 100),
-                "diagonal",
-                "diagonal",
+                "full",
+                "full",
             ),
             (
                 control.TransferFunction([1], [1, 2 * 0.5 * 1, 1**2]),
@@ -314,18 +302,6 @@ class TestComputeUncertaintyWeightResponse:
                 ],
                 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",
             ),
             (
@@ -396,8 +372,8 @@ class TestComputeUncertaintyWeightResponse:
                     ),
                 ],
                 np.logspace(-1, 1.5, 100),
-                "diagonal",
-                "diagonal",
+                "full",
+                "full",
             ),
             (
                 control.TransferFunction(
@@ -444,53 +420,6 @@ class TestComputeUncertaintyWeightResponse:
                 ],
                 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",
             ),
             (