Fix the spectral factorization (positive semidefiniteness) constraint for the power spectrum fit.

Author: timeverettadams

examples/7_uncertainty_characterization.py

@@ -40,10 +40,10 @@ def example_uncertainty_characterization():
 
     # Fit overbounding stable and minimum-phase uncertainty weight system
     weight_left = dkpy.fit_uncertainty_weight(
-        response_weight_left, omega, [4, 5], "diagonal"
+        response_weight_left, omega, [8, 8], "diagonal"
     )
     weight_right = dkpy.fit_uncertainty_weight(
-        response_weight_right, omega, [3, 5], "diagonal"
+        response_weight_right, omega, [8, 8], "diagonal"
     )
 
     # Plot: Magnitude response of nominal and off-nominal systems

src/dkpy/lti_system_fit.py

@@ -25,7 +25,7 @@ def fit_spectral_factor_mimo_ct(
     tol_bisection: float = 1e-5,
     max_iter_bisection: int = 100,
     max_iter_bisection_init: int = 15,
-    omega_power_constraint: Optional[np.ndarray] = None,
+    nbr_power_constraint: int = 500,
 ) -> control.StateSpace:
     """Fit a stable and minimum-phase biproper MIMO spectral factor to power spectrum.
 
@@ -52,8 +52,8 @@ def fit_spectral_factor_mimo_ct(
         Maximum number of iterations for the bisection algorithm.
     max_iter_bisection_init : int
         Maximum number of iterations for the bisection algorithm initialization.
-    omega_power_constraint : Optional[np.ndarray]
-        Angular frequencies to enforce non-negativity of power spectrum.
+    nbr_power_constraint : Optional[np.ndarray]
+        Number of frequencies to enforce non-negativity of power spectrum.
 
     Returns
     -------
@@ -81,11 +81,6 @@ def fit_spectral_factor_mimo_ct(
     # Pre-warp frequency (continuous- to discrete-time) with bilinear transformation
     alpha = np.sqrt(omega[0] * omega[-1])
     theta = 2 * np.atan(omega / alpha)
-    theta_power_constraint = (
-        2 * np.atan(omega_power_constraint / alpha)
-        if omega_power_constraint is not None
-        else None
-    )
 
     # Parse frequency-dependent weight
     if weight is None:
@@ -103,7 +98,7 @@ def fit_spectral_factor_mimo_ct(
         tol_bisection,
         max_iter_bisection,
         max_iter_bisection_init,
-        theta_power_constraint,
+        nbr_power_constraint,
     )
 
     # Discrete-time spectral factor fit model
@@ -133,7 +128,7 @@ def fit_power_mimo_dt(
     tol_bisection: float = 1e-5,
     max_iter_bisection: int = 100,
     max_iter_bisection_init: int = 15,
-    theta_power_constraint: Optional[np.ndarray] = None,
+    nbr_power_constraint: int = 500,
 ) -> Tuple[np.ndarray, np.ndarray]:
     """Fit a biproper discrete-time MIMO power spectrum transfer function.
 
@@ -160,9 +155,8 @@ def fit_power_mimo_dt(
         Maximum number of iterations for the bisection algorithm.
     max_iter_bisection_init : int
         Maximum number of iterations for the bisection algorithm initialization.
-    theta_power_constraint : Optional[np.ndarray]
-        Discrete-time frequencies to enforce non-negativity of power spectrum. The
-        frequencies range from [0, pi].
+    nbr_power_constraint : Optional[np.ndarray]
+        Number of frequencies to enforce non-negativity of power spectrum.
 
     Returns
     -------
@@ -177,9 +171,6 @@ def fit_power_mimo_dt(
     # Parse arguments
     solver_params = _parse_solver_param(solver_params)
     weight = np.ones_like(theta) if weight is None else weight
-    theta_power_constraint = (
-        theta if theta_power_constraint is None else theta_power_constraint
-    )
 
     # Auxiliary parameters
     nbr_signals = power_fit.shape[1]
@@ -210,7 +201,7 @@ def fit_power_mimo_dt(
         theta,
         order,
         weight,
-        theta_power_constraint,
+        nbr_power_constraint,
     )
 
     # Optimization objective
@@ -253,7 +244,7 @@ def _generate_constraints_power_mimo_dt(
     theta: np.ndarray,
     order: int,
     weight: np.ndarray,
-    theta_power_constraint: np.ndarray,
+    nbr_power_constraint: int,
 ) -> List[cvxpy.Constraint]:
     """Generate the constraints for the fit of discrete-time MIMO power spectrums.
 
@@ -281,9 +272,8 @@ def _generate_constraints_power_mimo_dt(
         Order of the power spectrum numerator and denominator polynomials.
     weight : np.ndarray,
         Frequency-dependent weight for fit accuracy.
-    theta_power_constraint : Optional[np.ndarray]
-        Discrete-time frequencies to enforce positive semidefiniteness of power
-        spectrum. The frequencies range from [0, pi].
+    nbr_power_constraint : Optional[np.ndarray]
+        Number of frequencies to enforce positive semidefiniteness of power spectrum.
 
     Returns
     -------
@@ -306,7 +296,7 @@ def _generate_constraints_power_mimo_dt(
 
     # Power spectrum positive definiteness constraints
     constraint_power_psd = _construct_power_psd_constraints_mimo_dt(
-        num_coef_sym_list, num_coef_skew_list, den_coef, theta_power_constraint, order
+        num_coef_sym_list, num_coef_skew_list, den_coef, order, nbr_power_constraint
     )
 
     # Numerator variable structure constraints
@@ -472,8 +462,8 @@ def _construct_power_psd_constraints_mimo_dt(
     num_coef_sym_list: List[cvxpy.Variable],
     num_coef_skew_list: List[cvxpy.Variable],
     den_coef: cvxpy.Variable,
-    theta: np.ndarray,
     order: int,
+    nbr_theta: int,
 ) -> List[cvxpy.Constraint]:
     """Construct the power spectrum positive semidefiniteness constraints.
 
@@ -503,6 +493,9 @@ def _construct_power_psd_constraints_mimo_dt(
     # Parameters
     nbr_signals = num_coef_sym_list[0].shape[0]
 
+    # Frequency points
+    theta = np.linspace(0, np.pi, nbr_theta)
+
     # Augmented variables
     num_coef_sym_stack = cvxpy.vstack(num_coef_sym_list)
     num_coef_skew_stack = cvxpy.vstack(num_coef_skew_list)
@@ -704,7 +697,7 @@ def fit_magnitude_siso_ct(
     tol_bisection: float = 1e-3,
     max_iter_bisection: int = 500,
     max_iter_bisection_init: int = 15,
-    omega_power_constraint: Optional[np.ndarray] = None,
+    nbr_power_constraint: int = 500,
 ) -> control.StateSpace:
     """Fit a stable and minimum-phase biproper SISO system to magnitude.
 
@@ -734,8 +727,8 @@ def fit_magnitude_siso_ct(
         Maximum allowable number of iterations in the bisection algorithm.
     max_iter_bisection_init : int
         Maximum number of iterations for the bisection algorithm initialization.
-    omega_power_constraint : Optional[np.ndarray]
-        Angular frequencies to enforce non-negativity of power spectrum.
+    nbr_power_constraint : Optional[np.ndarray]
+        Number of frequencies to enforce non-negativity of power spectrum.
 
     Returns
     -------
@@ -774,11 +767,6 @@ def fit_magnitude_siso_ct(
     # Pre-warp frequency (continuous- to discrete-time) with bilinear transformation
     alpha = np.sqrt(omega[0] * omega[-1])
     theta = 2 * np.atan(omega / alpha)
-    theta_power_constraint = (
-        2 * np.atan(omega_power_constraint / alpha)
-        if omega_power_constraint is not None
-        else None
-    )
 
     # Parse frequency-dependent weight
     if weight is None:
@@ -796,7 +784,7 @@ def fit_magnitude_siso_ct(
         tol_bisection,
         max_iter_bisection,
         max_iter_bisection_init,
-        theta_power_constraint,
+        nbr_power_constraint,
     )
 
     # Discrete-time spectral factor model
@@ -822,7 +810,7 @@ def fit_power_siso_dt(
     tol_bisection: float = 1e-5,
     max_iter_bisection: int = 100,
     max_iter_bisection_init: int = 15,
-    theta_power_constraint: Optional[np.ndarray] = None,
+    nbr_power_constraint: int = 500,
 ) -> Tuple[np.ndarray, np.ndarray]:
     """Fit a biproper discrete-time SISO power spectrum transfer function.
 
@@ -849,9 +837,8 @@ def fit_power_siso_dt(
         Maximum number of iterations for the bisection algorithm.
     max_iter_bisection_init : int
         Maximum number of iterations for the bisection algorithm initialization.
-    theta_power_constraint : Optional[np.ndarray]
-        Discrete-time frequencies to enforce non-negativity of power spectrum. The
-        frequencies range from [0, pi].
+    nbr_power_constraint : Optional[np.ndarray]
+        Number of frequencies to enforce non-negativity of power spectrum.
 
     Returns
     -------
@@ -882,9 +869,6 @@ def fit_power_siso_dt(
     # Parse arguments
     solver_params = _parse_solver_param(solver_params)
     weight = np.ones_like(theta) if weight is None else weight
-    theta_power_constraint = (
-        theta if theta_power_constraint is None else theta_power_constraint
-    )
 
     # Optimization variables and parameters
     num_coef = cvxpy.Variable(shape=order + 1, name="num_coef")
@@ -902,7 +886,7 @@ def fit_power_siso_dt(
         theta,
         order,
         weight,
-        theta_power_constraint,
+        nbr_power_constraint,
     )
 
     # Optimization problem
@@ -945,7 +929,7 @@ def _generate_constraints_power_siso_dt(
     theta: np.ndarray,
     order: int,
     weight: np.ndarray,
-    theta_power_constraints: np.ndarray,
+    nbr_power_constraints: int,
 ) -> List[cvxpy.Constraint]:
     """Generate the constraints for the discrete-time SISO power spectrums.
 
@@ -969,9 +953,8 @@ def _generate_constraints_power_siso_dt(
         Order of the power spectrum model.
     weight : np.ndarray
         Frequency-dependent weight for fit accuracy.
-    theta_power_constraint : Optional[np.ndarray]
-        Discrete-time frequencies to enforce non-negativity of power spectrum. The
-        frequencies range from [0, pi].
+    nbr_power_constraint : int
+        Discrete-time frequencies to enforce non-negativity of power spectrum.
 
     Returns
     -------
@@ -996,7 +979,7 @@ def _generate_constraints_power_siso_dt(
         num_coef,
         den_coef,
         order,
-        theta_power_constraints,
+        nbr_power_constraints,
     )
 
     # Normalization of the transfer function coefficients
@@ -1135,7 +1118,10 @@ def _construct_bound_constraints_siso_dt(
 
 
 def _construct_spectral_factorization_constraint_siso_dt(
-    num_coef: cvxpy.Variable, den_coef: cvxpy.Variable, order: int, theta: np.ndarray
+    num_coef: cvxpy.Variable,
+    den_coef: cvxpy.Variable,
+    order: int,
+    nbr_theta: np.ndarray,
 ) -> List[cvxpy.Constraint]:
     """Construct power spectrum spectral factorization constraints.
 
@@ -1150,6 +1136,8 @@ def _construct_spectral_factorization_constraint_siso_dt(
         Power spectrum denominator coefficient variable.
     order : int
         Order of the power spectrum model.
+    nbr_theta : int
+        Number of frequencies to enforce constraint.
 
     Returns
     -------
@@ -1158,6 +1146,9 @@ def _construct_spectral_factorization_constraint_siso_dt(
         constraints.
     """
 
+    # Frequencies to enforce power spectrum non-negativity
+    theta = np.linspace(0, np.pi, nbr_theta)
+
     constraint_spectral_factorization_list = []
     for theta_idx in theta:
         # Constraint data matrices

src/dkpy/uncertainty_characterization.py

@@ -12,7 +12,6 @@
     "plot_magnitude_response_uncertainty_weight",
 ]
 
-import warnings
 
 import control
 import numpy as np
@@ -787,7 +786,7 @@ def fit_uncertainty_weight(
     tol_bisection: float = 1e-3,
     max_iter_bisection: int = 500,
     max_iter_bisection_init: int = 15,
-    nbr_omega_power_constraint: int = 500,
+    nbr_power_constraint: int = 500,
 ) -> control.StateSpace:
     """Fit an overbounding stable and minimum-phase uncertainty weight.
 
@@ -888,9 +887,6 @@ def fit_uncertainty_weight(
             f"`{weight_structure}`."
         )
 
-    # Thing
-    omega_power_constraint = np.linspace(0, np.pi, nbr_omega_power_constraint)
-
     # Fit uncertainty weight
     uncertainty_weight = fit_uncertainty_weight_method(
         complex_uncertainty_weight,
@@ -901,7 +897,7 @@ def fit_uncertainty_weight(
         tol_bisection,
         max_iter_bisection,
         max_iter_bisection_init,
-        omega_power_constraint,
+        nbr_power_constraint,
     )
 
     return uncertainty_weight
@@ -916,7 +912,7 @@ def _fit_uncertainty_weight_scalar(
     tol_bisection: float = 1e-3,
     max_iter_bisection: int = 500,
     max_iter_bisection_init: int = 15,
-    omega_power_constraint: Optional[np.ndarray] = None,
+    nbr_power_constraint: int = 500,
 ) -> control.StateSpace:
     # Auxiliary parameters
     nbr_signals = complex_uncertainty_weight.shape[1]
@@ -938,7 +934,7 @@ def _fit_uncertainty_weight_scalar(
         tol_bisection=tol_bisection,
         max_iter_bisection=max_iter_bisection,
         max_iter_bisection_init=max_iter_bisection_init,
-        omega_power_constraint=omega_power_constraint,
+        nbr_power_constraint=nbr_power_constraint,
     )
 
     # Duplicate the scalar uncertainty weight along the diagonal
@@ -956,7 +952,7 @@ def _fit_uncertainty_weight_diagonal(
     tol_bisection: float = 1e-3,
     max_iter_bisection: int = 500,
     max_iter_bisection_init: int = 15,
-    omega_power_constraint: Optional[np.ndarray] = None,
+    nbr_power_constraint: int = 500,
 ):
     # Auxiliary parameters
     nbr_signals = complex_uncertainty_weight.shape[1]
@@ -969,7 +965,6 @@ def _fit_uncertainty_weight_diagonal(
     for idx in range(nbr_signals):
         # Compute the magnitude of the diagonal uncertainty weight element
         magnitude_uncertainty_weight = np.abs(complex_uncertainty_weight[:, idx, idx])
-        print(magnitude_uncertainty_weight)
         order_idx = order_list[idx]
 
         # Fit the diagonal uncertainty weight element
@@ -984,7 +979,7 @@ def _fit_uncertainty_weight_diagonal(
             tol_bisection=tol_bisection,
             max_iter_bisection=max_iter_bisection,
             max_iter_bisection_init=max_iter_bisection_init,
-            omega_power_constraint=omega_power_constraint,
+            nbr_power_constraint=nbr_power_constraint,
         )
         uncertainty_weight_list.append(uncertainty_weight_diagonal)
 
@@ -994,6 +989,7 @@ def _fit_uncertainty_weight_diagonal(
     return uncertainty_weight
 
 
+# TODO: Fit spectral factorization constraint generation
 def _fit_uncertainty_weight_full(
     complex_uncertainty_weight: np.ndarray,
     omega: np.ndarray,