diff --git a/doc/_static/example_7/7_sval_uncertainty_weight_left.png b/doc/_static/example_7/7_sval_uncertainty_weight_left.png new file mode 100644 index 0000000..888db41 Binary files /dev/null and b/doc/_static/example_7/7_sval_uncertainty_weight_left.png differ diff --git a/doc/_static/example_7/7_sval_uncertainty_weight_right.png b/doc/_static/example_7/7_sval_uncertainty_weight_right.png new file mode 100644 index 0000000..4d807b2 Binary files /dev/null and b/doc/_static/example_7/7_sval_uncertainty_weight_right.png differ diff --git a/doc/_static/example_7/7_uncertainty_weight.png b/doc/_static/example_7/7_uncertainty_weight.png deleted file mode 100644 index dccd187..0000000 Binary files a/doc/_static/example_7/7_uncertainty_weight.png and /dev/null differ diff --git a/doc/_static/example_8/8_magnitude_uncertainty_weight.png b/doc/_static/example_8/8_magnitude_uncertainty_weight.png deleted file mode 100644 index e3d1718..0000000 Binary files a/doc/_static/example_8/8_magnitude_uncertainty_weight.png and /dev/null differ diff --git a/doc/_static/example_8/8_sval_uncertainty_weight_left.png b/doc/_static/example_8/8_sval_uncertainty_weight_left.png new file mode 100644 index 0000000..ef0d974 Binary files /dev/null and b/doc/_static/example_8/8_sval_uncertainty_weight_left.png differ diff --git a/doc/examples.rst b/doc/examples.rst index a1326a7..c961238 100644 --- a/doc/examples.rst +++ b/doc/examples.rst @@ -204,17 +204,17 @@ by the variation in parameters. Next, the uncertainty residual can be computed for each off-nominal system and uncertainty model. `dkpy` implements six unstructured uncertainty models: -* Additive uncertainty ("A"); +* Additive uncertainty ("A") .. image:: _static/example_7/7_sval_residual_A.png -* Multiplicative input uncertainty ("I"); +* Multiplicative input uncertainty ("I") .. image:: _static/example_7/7_sval_residual_I.png -* Multiplicative output uncertainty ("O"); +* Multiplicative output uncertainty ("O") .. image:: _static/example_7/7_sval_residual_O.png -* Inverse additive uncertainty ("iA"); +* Inverse additive uncertainty ("iA") .. image:: _static/example_7/7_sval_residual_iA.png -* Inverse multiplicative input uncertainty ("iI"); +* Inverse multiplicative input uncertainty ("iI") .. image:: _static/example_7/7_sval_residual_iI.png -* Inverse multiplicative output uncertainty ("iO"). +* 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 @@ -239,7 +239,10 @@ 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 +* Left: + .. image:: _static/example_7/7_sval_uncertainty_weight_left.png +* Right: + .. image:: _static/example_7/7_sval_uncertainty_weight_right.png Multi-Model Uncertainty Characterization: Aircraft Actuator Model ----------------------------------------------------------------- @@ -269,7 +272,7 @@ weight is constrained to the identity matrix. Given that the right uncertainty weight is the identity matrix, a fit does not need to be performed for this weight as it will be neglected in the generalized plant. -.. image:: _static/example_8/8_magnitude_uncertainty_weight.png +.. image:: _static/example_8/8_sval_uncertainty_weight_left.png diff --git a/examples/7_uncertainty_characterization.py b/examples/7_uncertainty_characterization.py index eec2817..6e35a44 100644 --- a/examples/7_uncertainty_characterization.py +++ b/examples/7_uncertainty_characterization.py @@ -32,15 +32,26 @@ def example_uncertainty_characterization(): ) # Compute uncertainty weight frequency response - response_weight_left, response_weight_right = ( + complex_weight_left, complex_weight_right = ( dkpy.compute_uncertainty_weight_response( response_residuals_dict["inverse_additive"], "diagonal", "diagonal" ) ) # Fit overbounding stable and minimum-phase uncertainty weight system - weight_left = dkpy.fit_uncertainty_weight(response_weight_left, omega, [4, 5]) - weight_right = dkpy.fit_uncertainty_weight(response_weight_right, omega, [3, 5]) + weight_left = dkpy.fit_uncertainty_weight( + complex_weight_left, omega, [4, 5], "left", "diagonal" + ) + weight_right = dkpy.fit_uncertainty_weight( + complex_weight_right, omega, [3, 5], "right", "diagonal" + ) + + measure = dkpy.compute_uncertainty_measure_response( + complex_nominal, + complex_weight_left, + complex_weight_right, + "multiplicative_input", + ) # Plot: Magnitude response of nominal and off-nominal systems dkpy.plot_magnitude_response_uncertain_model_set( @@ -72,15 +83,20 @@ def example_uncertainty_characterization(): response_residuals_dict, omega ) - # Plot: Magnitude response of uncertainty weight frequency response and overbounding - # fit - dkpy.plot_magnitude_response_uncertainty_weight( - response_weight_left, - response_weight_right, + # Plot: Singular value response of left uncertainty weight + dkpy.plot_singular_value_response_uncertainty_weight( + complex_weight_left, omega, weight_left, + ) + + # Plot: Singular value response of right uncertainty weight + dkpy.plot_singular_value_response_uncertainty_weight( + complex_weight_right, + omega, weight_right, ) + plt.show() diff --git a/examples/8_aircraft_actuator_uncertainty_characterization.py b/examples/8_aircraft_actuator_uncertainty_characterization.py index 058b226..bc3b315 100644 --- a/examples/8_aircraft_actuator_uncertainty_characterization.py +++ b/examples/8_aircraft_actuator_uncertainty_characterization.py @@ -26,7 +26,7 @@ def example_aircraft_uncertainty_characterization(): ) # Compute the optimal uncertainty weights with a given structure - response_weight_left, response_weight_right = ( + complex_weight_left, complex_weight_right = ( dkpy.compute_uncertainty_weight_response( response_residual["multiplicative_input"], "scalar", @@ -35,7 +35,9 @@ def example_aircraft_uncertainty_characterization(): ) # Fit an overbounding LTI system to the optimal uncertainty weight response - weight_left = dkpy.fit_uncertainty_weight(response_weight_left, omega, 1) + weight_left = dkpy.fit_uncertainty_weight( + complex_weight_left, omega, 1, "left", "scalar" + ) # Plot: Singular value response of nominal and off-nominal systems dkpy.plot_singular_value_response_uncertain_model_set( @@ -48,15 +50,13 @@ def example_aircraft_uncertainty_characterization(): response_residual, omega, hz=True ) - # Plot: Magnitude response of uncertainty weight frequency response and overbounding - # fit - dkpy.plot_magnitude_response_uncertainty_weight( - response_weight_left, - response_weight_right, + # Plot: Singular value response of left uncertainty weight + dkpy.plot_singular_value_response_uncertainty_weight( + complex_weight_left, omega, - weight_left=weight_left, - hz=True, + weight_left, ) + plt.show() diff --git a/src/dkpy/__init__.py b/src/dkpy/__init__.py index 4a48a14..6a944bb 100644 --- a/src/dkpy/__init__.py +++ b/src/dkpy/__init__.py @@ -3,6 +3,7 @@ from .controller_synthesis import * from .dk_iteration import * from .d_scale_fit import * +from .lti_system_fit import * from .structured_singular_value import * from .uncertainty_structure import * from .uncertainty_characterization import * diff --git a/src/dkpy/lti_system_fit.py b/src/dkpy/lti_system_fit.py new file mode 100644 index 0000000..5b27af1 --- /dev/null +++ b/src/dkpy/lti_system_fit.py @@ -0,0 +1,739 @@ +"""LTI system fit utilities.""" + +from typing import Any, Dict, Union, Optional, Tuple, List + +import numpy as np +import control +import cvxpy + +__all__ = [ + "fit_magnitude_siso_ct", + "fit_power_siso_dt", +] + + +def fit_magnitude_siso_ct( + magnitude_fit: np.ndarray, + omega: np.ndarray, + order: int, + magnitude_upper_bound: Optional[np.ndarray] = None, + magnitude_lower_bound: Optional[np.ndarray] = None, + weight: Optional[np.ndarray] = None, + solver_params: Optional[Dict[str, Any]] = None, + tol_bisection: float = 1e-3, + max_iter_bisection: int = 500, + max_iter_bisection_init: int = 15, + nbr_power_constraint: int = 500, +) -> control.StateSpace: + """Fit a stable and minimum-phase biproper SISO system to magnitude. + + 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. + solver_param : Dict[str, Any] + Solver parameters for the optimization problem. These are keyword arguments for + `cvxpy.Problem.solve()` [#cvxpy_solver]_. + tol_bisection : float + Numerical tolerance for the bisection algorithm. + max_iter_bisection : int + Maximum allowable number of iterations in the bisection algorithm. + max_iter_bisection_init : int + Maximum number of iterations for the bisection algorithm initialization. + nbr_power_constraint : Optional[np.ndarray] + Number of frequencies to enforce non-negativity of power spectrum. + + Returns + ------- + control.StateSpace + Fitted state-space system. + + Notes + ----- + The algorithm fits the LTI system to the magnitude data as follows: + 1) Pre-warp the continuous-time frequencies to discrete-time frequencies using a + bilinear transformation. + 2) Fit a discrete-time power spectrum model to the square of the magnitude data. + 3) Extract the stable and minimum-phase spectral factor from the power spectrum + model. + 4) Convert the discrete-time spectral factor model to a continuous-time model + using the bilinear transformation used in step 1). + + References + ---------- + .. [#cxvpy_solver] https://www.cvxpy.org/tutorial/solvers/index.html + """ + + # Compute power spectrum + power_fit = np.abs(magnitude_fit) ** 2 + power_upper_bound = ( + np.abs(magnitude_upper_bound) ** 2 + if magnitude_upper_bound is not None + else None + ) + power_lower_bound = ( + np.abs(magnitude_lower_bound) ** 2 + if magnitude_lower_bound is not None + else None + ) + + # Pre-warp frequency (continuous- to discrete-time) with bilinear transformation + alpha = np.sqrt(omega[0] * omega[-1]) + theta = 2 * np.atan(omega / alpha) + + # Parse frequency-dependent weight + if weight is None: + weight = np.ones_like(omega) + + # Discrete-time power spectrum fit coefficients + num_power_coef, den_power_coef = fit_power_siso_dt( + power_fit, + theta, + order, + weight, + power_upper_bound, + power_lower_bound, + solver_params, + tol_bisection, + max_iter_bisection, + max_iter_bisection_init, + nbr_power_constraint, + ) + + # Discrete-time spectral factor model + num_factor_coef, den_factor_coef = _compute_spectral_factor_siso_dt( + num_power_coef, den_power_coef + ) + tf_factor_dt = control.TransferFunction(num_factor_coef, den_factor_coef, True) + + # Continous-time state-space + ss_factor_ct = _convert_discrete_to_continuous_bilinear(tf_factor_dt, alpha) + + return ss_factor_ct + + +def fit_power_siso_dt( + power_fit: np.ndarray, + theta: np.ndarray, + order: int, + weight: Optional[np.ndarray] = None, + power_upper_bound: Optional[np.ndarray] = None, + power_lower_bound: Optional[np.ndarray] = None, + solver_params: Optional[Dict[str, Any]] = None, + tol_bisection: float = 1e-5, + max_iter_bisection: int = 100, + max_iter_bisection_init: int = 15, + nbr_power_constraint: int = 500, +) -> Tuple[np.ndarray, np.ndarray]: + """Fit a biproper discrete-time SISO power spectrum transfer function. + + Parameters + ---------- + power_fit : np.ndarray + Power spectrum data to fit model. + theta : np.ndarray + Discrete-time angular frequencies (-). The frequencies range from [0, pi]. + order : int + Order of the power spectrum model. + weight : Optional[np.ndarray] + Frequency-dependent weight for fit accuracy. + power_upper_bound : np.ndarray + Power spectrum data to upper bound the model frequency response. + power_lower_bound : np.ndarray + Power spectrum data to lower bound the model frequency response. + solver_params: Optional[Dict[str, Any]] + Solver parameters for the optimization problem. These are keyword arguments for + `cvxpy.Problem.solve()` [#cvxpy_solver]_. + tol_bisection : float + Numerical tolerance for the bisection algorithm. + max_iter_bisection : int + Maximum number of iterations for the bisection algorithm. + max_iter_bisection_init : int + Maximum number of iterations for the bisection algorithm initialization. + nbr_power_constraint : Optional[np.ndarray] + Number of frequencies to enforce non-negativity of power spectrum. + + Returns + ------- + Tuple[np.ndarray, np.ndarray] + Numerator and denominator polynomial coefficients of the power spectrum. + + Notes + ----- + The algorithm fits a power spectrum of the form + b_{m} z^{-m} + ... + b_{0} + ... + b_{m} z^{m} + R(z) = ---------------------------------------------- + a_{m} z^{-m} + ... + a_{0} + ... + b_{m} z^{m} + + The algorithm uses a log-Chebyshev approximation method, which results in an + optimization problem that is solved by bisecting a bound on the fit error and + solving linear feasibility problems [#log_chebyshev]_. + + References + ---------- + .. [#cxvpy_solver] https://www.cvxpy.org/tutorial/solvers/index.html + + .. [#log_chebyshev] Shao-Po Wu, S. Boyd and L. Vandenberghe, "FIR filter design via + semidefinite programming and spectral factorization," Proceedings of 35th IEEE + Conference on Decision and Control, Kobe, Japan, 1996, pp. 271-276 vol.1, + doi: 10.1109/CDC.1996.574313. + """ + + # Parse arguments + solver_params = _parse_solver_param(solver_params) + weight = np.ones_like(theta) if weight is None else weight + + # Optimization variables and parameters + num_coef = cvxpy.Variable(shape=order + 1, name="num_coef") + den_coef = cvxpy.Variable(shape=order + 1, name="den_coef") + error_bound = cvxpy.Parameter(shape=(), name="error_bound") + + # Constraints + constraint_list = _generate_constraints_power_siso_dt( + num_coef, + den_coef, + error_bound, + power_fit, + power_upper_bound, + power_lower_bound, + theta, + order, + weight, + nbr_power_constraint, + ) + + # Optimization problem + objective = cvxpy.Minimize(1) + problem = cvxpy.Problem(objective, constraint_list) + + # Bisection solution + error_bound_max = _initialize_optimization_bisection_upper_bound( + problem, error_bound, solver_params, max_iter_bisection_init + ) + error_bound_min = 0 + _solve_optimization_bisection( + problem, + error_bound, + error_bound_max, + error_bound_min, + tol_bisection, + max_iter_bisection, + solver_params, + ) + + # Extract power spectrum coefficients + num_coef_fit = np.array(num_coef.value) + den_coef_fit = np.array(den_coef.value) + + # Construct power spectrum polynomials from unique coefficients + num_fit = np.concatenate([num_coef_fit[::-1], num_coef_fit[1:]]) + den_fit = np.concatenate([den_coef_fit[::-1], den_coef_fit[1:]]) + + return num_fit, den_fit + + +def _generate_constraints_power_siso_dt( + num_coef: cvxpy.Variable, + den_coef: cvxpy.Variable, + error_bound: cvxpy.Parameter, + power_fit: np.ndarray, + power_upper_bound: Optional[np.ndarray], + power_lower_bound: Optional[np.ndarray], + theta: np.ndarray, + order: int, + weight: np.ndarray, + nbr_power_constraints: int, +) -> List[cvxpy.Constraint]: + """Generate the constraints for the discrete-time SISO power spectrums. + + Parameters + ---------- + num_coef : cvxpy.Variable + Power spectrum numerator coefficient variable. + den_coef : cvxpy.Variable + Power spectrum denominator coefficient variable. + error_bound : cvxpy.Parameter + Power spectrum fit error bound that is minimized by the bisection algorithm. + power_fit: np.ndarray + Power spectrum data to fit model. + power_upper_bound : np.ndarray + Power spectrum data to upper bound the model frequency response. + power_lower_bound : np.ndarray + Power spectrum data to lower bound the model frequency response. + theta : np.ndarray + Discrete-time angular frequencies (-). The frequencies range from [0, pi]. + order : int + Order of the power spectrum model. + weight : np.ndarray + Frequency-dependent weight for fit accuracy. + nbr_power_constraint : int + Discrete-time frequencies to enforce non-negativity of power spectrum. + + Returns + ------- + List[cvxpy.Constraint] + List of constraints. + """ + # Fit error and magnitude bound constraints + constraint_bound_list = _construct_bound_constraints_siso_dt( + num_coef, + den_coef, + error_bound, + power_fit, + power_upper_bound, + power_lower_bound, + theta, + order, + weight, + ) + + # Spectral factorization constraint + constraint_nonneg_list = _construct_spectral_factorization_constraint_siso_dt( + num_coef, + den_coef, + order, + nbr_power_constraints, + ) + + # Normalization of the transfer function coefficients + constraint_normalization_list = [den_coef[-1] == 1] + + constraint_list = [] + constraint_list += ( + constraint_bound_list + constraint_nonneg_list + constraint_normalization_list + ) + + return constraint_list + + +def _construct_bound_constraints_siso_dt( + num_coef: cvxpy.Variable, + den_coef: cvxpy.Variable, + error_bound: cvxpy.Parameter, + power_fit: np.ndarray, + power_upper_bound: Optional[np.ndarray], + power_lower_bound: Optional[np.ndarray], + theta: np.ndarray, + order: int, + weight: np.ndarray, +) -> List[cvxpy.Constraint]: + """Construct the fit error and power spectrum bound linear constraints. + + The fit error and power spectrum bound linear constraints are used in the fit of + discrete-time SISO power spectrums. + + Parameters + ---------- + num_coef : cvxpy.Variable + Power spectrum numerator coefficient variable. + den_coef : cvxpy.Variable + Power spectrum denominator coefficient variable. + error_bound : cvxpy.Parameter + Fit error parameter that is minimized by the bisection algorithm. + power_fit: np.ndarray + Power spectrum data to fit model. + power_upper_bound : np.ndarray + Power spectrum data to upper bound the model frequency response. + power_lower_bound : np.ndarray + Power spectrum data to lower bound the model frequency response. + theta : np.ndarray + Discrete-time angular frequencies (-). The frequencies range from [0, pi]. + order : int + Order of the power spectrum model. + weight : np.ndarray + Frequency-dependent weight for fit accuracy. + + Returns + ------- + List[cvxpy.Constraint] + List of fit error and magnitude bound linear constraints. + """ + + # Initialize constraint lists + constraint_fit_upper_list = [] + constraint_fit_lower_list = [] + constraint_bound_upper_list = [] + constraint_bound_lower_list = [] + + # Check if the upper/lower bound magnitude is the same as the fit magnitude. If it + # is, the upper/lower fit error constraint may be removed as it is made redundant by + # the bound constraint. + is_power_fit_upper_bound = ( + np.allclose(power_fit, power_upper_bound) + if power_upper_bound is not None + else False + ) + is_power_fit_lower_bound = ( + np.allclose(power_fit, power_lower_bound) + if power_lower_bound is not None + else False + ) + + # Construct constraints at each frequency + for idx in range(theta.size): + # Constraint data + theta_idx = theta[idx] + power_idx = power_fit[idx] + weight_idx = weight[idx] + + # Constraint data matrices + cosine_row = np.array( + [2 * np.cos(idx_coef * theta_idx) for idx_coef in range(order + 1)] + ) + cosine_row[0] = 1 + + # Fit error upper bound constraint + if not is_power_fit_upper_bound: + constraint_fit_upper = ( + -(1 + error_bound / weight_idx) + * power_idx + * cvxpy.vdot(cosine_row, den_coef) + + cvxpy.vdot(cosine_row, num_coef) + <= 0 + ) + constraint_fit_upper_list.append(constraint_fit_upper) + + # Fit error lower bound constraint + if not is_power_fit_lower_bound: + constraint_fit_lower = ( + power_idx * cvxpy.vdot(cosine_row, den_coef) + - (1 + error_bound / weight_idx) * cvxpy.vdot(cosine_row, num_coef) + <= 0 + ) + constraint_fit_lower_list.append(constraint_fit_lower) + + # Upper bound constraint + if power_upper_bound is not None: + power_upper_idx = power_upper_bound[idx] + constraint_bound_upper = ( + -power_upper_idx * cvxpy.vdot(cosine_row, den_coef) + - cvxpy.vdot(cosine_row, num_coef) + <= 0 + ) + constraint_bound_upper_list.append(constraint_bound_upper) + + # Lower bound constraint + if power_lower_bound is not None: + power_lower_idx = power_lower_bound[idx] + constraint_bound_lower = ( + power_lower_idx * cvxpy.vdot(cosine_row, den_coef) + - cvxpy.vdot(cosine_row, num_coef) + <= 0 + ) + constraint_bound_lower_list.append(constraint_bound_lower) + + return ( + constraint_fit_upper_list + + constraint_fit_lower_list + + constraint_bound_upper_list + + constraint_bound_lower_list + ) + + +def _construct_spectral_factorization_constraint_siso_dt( + num_coef: cvxpy.Variable, + den_coef: cvxpy.Variable, + order: int, + nbr_theta: np.ndarray, +) -> List[cvxpy.Constraint]: + """Construct power spectrum spectral factorization constraints. + + The power spectrum spectral factorization constraints are used in the fit of + discrete-time SISO power spectrums. + + Parameters + ---------- + num_coef : cvxpy.Variable + Power spectrum numerator coefficient variable. + den_coef : cvxpy.Variable + Power spectrum denominator coefficient variable. + order : int + Order of the power spectrum model. + nbr_theta : int + Number of frequencies to enforce constraint. + + Returns + ------- + List[cvxpy.Constraint] + List of spectral factorization (power spectrum non-negativity) linear + 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 + cosine_row = np.array( + [2 * np.cos(idx_coef * theta_idx) for idx_coef in range(order + 1)] + ) + cosine_row[0] = 1 + + # Spectral factorization constraints + constraint_num = cvxpy.vdot(cosine_row, num_coef) >= 0 + constraint_den = cvxpy.vdot(cosine_row, den_coef) >= 0 + constraint_spectral_factorization_list += [constraint_num, constraint_den] + + return constraint_spectral_factorization_list + + +def _initialize_optimization_bisection_upper_bound( + problem: cvxpy.Problem, + objective: cvxpy.Parameter, + solver_params: Dict[str, Any], + max_iter_bisection_init: int = 15, +) -> float: + """Initialize the upper bound of a bisection optimization problem. + + Parameters + ---------- + problem : cvxpy.Problem + Optimization problem. + objective : cvxpy.Parameter + Objective function that the bisection algorithm minimizes. + solver_params : Dict[str, Any] + Solver parameters for the optimization problem. These are keyword arguments for + `cvxpy.Problem.solve()` [#cvxpy_solver]_. + max_iter_bisection_init : int + Maximum allowable number of iterations in the bisection initialization. + + Returns + ------- + float + Feasible upper bound for the bisection optimization algorithm. + + References + ---------- + .. [#cxvpy_solver] https://www.cvxpy.org/tutorial/solvers/index.html + """ + objective.value = 0.5 + + for _ in range(max_iter_bisection_init): + objective.value *= 2 + try: + problem.solve(**solver_params) + feasibility_status = problem.status + except cvxpy.SolverError: + break + + if feasibility_status == "optimal": + return objective.value + + raise ValueError( + "Unable to determine a upper bound on the objective to initialize " + "the bisection algorithm." + ) + + +def _solve_optimization_bisection( + problem: cvxpy.Problem, + objective: cvxpy.Parameter, + objective_max: float, + objective_min: float, + tol_bisection: float, + max_iter_bisection: int, + solver_params: Dict[str, Any], +) -> None: + """Solve an optimization problem using a bisection algorithm. + + Parameters + ---------- + problem : cvxpy.Problem + Optimization problem. + objective : cvxpy.Parameter + Objective function that the bisection algorithm minimizes. + objective_max : float + Initial upper bound on objective. + objective_min : float + Initial lower bound on objective. + tol_bisection : float + Numerical tolerance for the bisection algorithm. + max_iter_bisection : int + Maximum allowable number of iterations in the bisection algorithm. + solver_params : Dict[str, Any] + Solver parameters for the optimization problem. These are keyword arguments for + `cvxpy.Problem.solve()` [#cvxpy_solver]_. + + References + ---------- + .. [#cxvpy_solver] https://www.cvxpy.org/tutorial/solvers/index.html + """ + + # Initialize bisection algorithm + iter_bisection = 0 + feasibility_status = "infeasible" + + # Bisection minimization of objective + while ( + np.abs(objective_max - objective_min) >= tol_bisection + or feasibility_status != "optimal" + ): + # Increment bisection iteration count + iter_bisection += 1 + + # Stop bisection at maximum iterations + if iter_bisection >= max_iter_bisection: + break + + # Bisect objective function + objective.value = 0.5 * (objective_max + objective_min) + try: + problem.solve(**solver_params) + feasibility_status = problem.status + except cvxpy.SolverError: + objective_min = objective.value + + # Update upper and lower bisection bounds + if feasibility_status == "optimal": + objective_max = objective.value + else: + objective_min = objective.value + + +def _compute_spectral_factor_siso_dt( + num_power: np.ndarray, den_power: np.ndarray +) -> Tuple[np.ndarray, np.ndarray]: + """Compute the spectral factorization of a discrete-time SISO transfer function. + + Parameters + ---------- + num_power : np.ndarray + Power spectrum numerator polynomial. + den_power : np.ndarray + Power spectrum denominator polynomial. + + Returns + ------- + Tuple[np.ndarray, np.ndarray] + Spectral factor numerator and denominator polynomial. + + Notes + ----- + The algorithm takes a discrete-time power spectrum of the form + b_{m} z^{-m} + ... + b_{0} + ... + b_{m} z^{m} + R(z) = ---------------------------------------------- + a_{m} z^{-m} + ... + a_{0} + ... + b_{m} z^{m} + + and computes the spectral factor of the form + d_{m} z^{m} + ... + d_{0} + F(z) = ------------------------- + c_{m} z^{m} + ... + c_{0} + + The algorithm computes the spectral factor by computing the poles and zeros of the + power spectrum. Then, the DC gain, stable poles, and minimum-phase zeros are + extracted and used to construct the stable and minimum-phase spectral factor. + """ + + # Power spectrum poles and zeros + zeros_power = np.roots(num_power) + poles_power = np.roots(den_power) + + # Stable and minimum-phase power spectrum poles and zeros + zeros_power_stable = zeros_power[np.abs(zeros_power) < 1] + poles_power_stable = poles_power[np.abs(poles_power) < 1] + + # Spectral factorization + spectral_const_num = np.sqrt(num_power[0] / np.prod(-zeros_power_stable)) + spectral_const_den = np.sqrt(den_power[0] / np.prod(-poles_power_stable)) + num_factor = np.real( + spectral_const_num + * np.polynomial.polynomial.polyfromroots(zeros_power_stable)[::-1] + ) + den_factor = np.real( + spectral_const_den + * np.polynomial.polynomial.polyfromroots(poles_power_stable)[::-1] + ) + + return num_factor, den_factor + + +def _convert_discrete_to_continuous_bilinear( + sys_dt: Union[control.TransferFunction, control.StateSpace], + alpha: float, +) -> control.StateSpace: + """Convert a discrete-time system to continuous-time with a bilinear transformation. + + Parameters + ---------- + sys_dt : Union[control.TransferFunction, control.StateSpace] + Discrete-time system. + 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]) + In_Ad_inv = np.linalg.solve(In + Ad, In) + + # Continous-time matrices + Ac = alpha * (Ad - In) @ In_Ad_inv + Bc = alpha * (In - (Ad - In) @ In_Ad_inv) @ Bd + Cc = Cd @ In_Ad_inv + Dc = Dd - Cd @ In_Ad_inv @ Bd + + # Continuous-time system + sys_dt = control.StateSpace(Ac, Bc, Cc, Dc) + + return sys_dt + + +def _parse_solver_param(solver_params: Optional[Dict[str, Any]]) -> Dict[str, Any]: + """Parse solver parameters and ensure warm-start is enabled. + + Parameters + ---------- + solver_params: Optional[Dict[str, Any]] + Solver parameters for the optimization problem. These are keyword arguments for + `cvxpy.Problem.solve()` [#cvxpy_solver]_. + + Returns + ------- + Dict[str, Any] + Parsed solver parameters for the optimization problem. These are keyword + arguments for `cvxpy.Problem.solve()` [#cvxpy_solver]_. + + References + ---------- + .. [#cxvpy_solver] https://www.cvxpy.org/tutorial/solvers/index.html + """ + # Default solver parameters + solver_param_default = { + "solver": "CLARABEL", + "warm_start": True, + "verbose": False, + } + + if solver_params is None: + return solver_param_default + else: + # Ensure that `warm_start` is set to True as the bisection algorithm is greatly + # sped up if this holds. + solver_params["warm_start"] = True + return solver_params diff --git a/src/dkpy/uncertainty_characterization.py b/src/dkpy/uncertainty_characterization.py index f9787cd..9b6a7e0 100644 --- a/src/dkpy/uncertainty_characterization.py +++ b/src/dkpy/uncertainty_characterization.py @@ -3,16 +3,17 @@ __all__ = [ "compute_uncertainty_residual_response", "compute_uncertainty_weight_response", + "compute_uncertainty_measure_response", "fit_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", + "plot_singular_value_response_uncertainty_weight", + "plot_uncertainty_measure", ] -import warnings import control import numpy as np @@ -20,25 +21,18 @@ 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, Any, Literal from matplotlib.figure import Figure from matplotlib.axes import Axes from matplotlib.legend import Legend -from . import utilities +from . import lti_system_fit 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", - }, + uncertainty_model: Union[str, List[str]], tol_residual_existence: float = 1e-12, ) -> Dict[str, np.ndarray]: """Compute the residual response of unstructured uncertainty models. @@ -49,8 +43,15 @@ def compute_uncertainty_residual_response( 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. + uncertainty_model : Union[str, List[str]] + Uncertainty model identifiers to compute the residual response. The valid + uncertainty model identifiers are: + - "additive" + - "multiplicative_input" + - "multiplicative_output" + - "inverse_additive" + - "inverse_multiplicative_input", + - "inverse_multiplicative_output" tol_residual_existence : float Tolerance for the existence of an uncertainty residual. @@ -72,14 +73,14 @@ def compute_uncertainty_residual_response( >>> complex_response_nom, complex_response_offnom_list, omega = ( ... example_multimodel_uncertainty ... ) - >>> uncertainty_models = { + >>> 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, @@ -98,202 +99,82 @@ def compute_uncertainty_residual_response( # 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", + compute_residual_dispatcher = { + "additive": _compute_residual_additive, + "multiplicative_input": _compute_residual_multiplicative_input, + "multiplicative_output": _compute_residual_multiplicative_output, + "inverse_additive": _compute_residual_inverse_additive, + "inverse_multiplicative_input": _compute_residual_inverse_multiplicative_input, + "inverse_multiplicative_output": _compute_residual_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( + for model in uncertainty_model: + try: + compute_residual_model = compute_residual_dispatcher[model] + except KeyError: + raise KeyError( + 'The uncertainty model identifier must be "additive", ' + '"multiplicative_input", "multiplicative_output", "inverse_additive" ' + '"inverse_multiplicative_input", or "inverse_multiplicative_output" ' + f'(got "{model}").' + ) + complex_response_residual_list = compute_residual_model( 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 - ) + complex_response_residual_dict[model] = 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, +def _compute_residual_additive( + complex_nominal: np.ndarray, + complex_offnominal: np.ndarray, tol_residual_existence: Optional[float] = None, ) -> np.ndarray: - """Compute the additive uncertainty residual at a frequency. + """Compute the additive uncertainty residual frequency responses. 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. + complex_nominal : np.ndarray + Nominal model frequency response. + complex_offnominal : np.ndarray + Off-nominal model frequency responses. tol_residual_existence : float - Tolerance for the existence of a multiplicative input uncertainty residual. + Tolerance for the existence of an uncertainty residual. Returns ------- np.ndarray - Additive uncertainty residual at a given frequency. + Additive uncertainty residual frequency response for all off-nominal models. """ - complex_response_residual_freq = ( - complex_response_offnom_freq - complex_response_nom_freq - ) + complex_residual = complex_offnominal - complex_nominal - return complex_response_residual_freq + return complex_residual -def _compute_uncertainty_residual_multiplicative_input_freq( - complex_response_nom_freq: np.ndarray, - complex_response_offnom_freq: np.ndarray, +def _compute_residual_multiplicative_input( + complex_nominal: np.ndarray, + complex_offnominal: np.ndarray, tol_residual_existence: float = 1e-8, ) -> np.ndarray: - """Compute the multiplicative input uncertainty residual at a frequency. + """Compute the multiplicative input uncertainty residual frequency responses. 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. + complex_response_nom : np.ndarray + Nominal model frequency response. + complex_response_offnom : np.ndarray + Off-nominal model frequency responses. tol_residual_existence : float - Tolerance for the existence of a multiplicative input uncertainty residual. + Tolerance for the existence of an uncertainty residual. Returns ------- np.ndarray - Multiplicative input uncertainty residual at a given frequency. + Multiplicative input uncertainty residual frequency response for all + off-nominal models. Raises ------ @@ -303,19 +184,19 @@ def _compute_uncertainty_residual_multiplicative_input_freq( solution. """ - num_inputs = complex_response_nom_freq.shape[1] - num_outputs = complex_response_nom_freq.shape[0] + nbr_inputs = complex_nominal.shape[-1] + nbr_outputs = complex_nominal.shape[-2] - 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 + a = complex_nominal + b = complex_offnominal - complex_nominal + x, residues_lstsq, _, _ = scipy.linalg.lstsq(a, b) + complex_residual = x - if num_inputs >= num_outputs: - return complex_response_residual_freq + if nbr_inputs >= nbr_outputs: + return complex_residual else: if np.all(residues_lstsq <= tol_residual_existence): - return complex_response_residual_freq + return complex_residual else: raise ValueError( "A multiplicative input uncertainty residual does not exist for the " @@ -327,27 +208,27 @@ def _compute_uncertainty_residual_multiplicative_input_freq( ) -def _compute_uncertainty_residual_multiplicative_output_freq( - complex_response_nom_freq: np.ndarray, - complex_response_offnom_freq: np.ndarray, +def _compute_residual_multiplicative_output( + complex_nominal: np.ndarray, + complex_offnominal: np.ndarray, tol_residual_existence: float = 1e-8, ) -> np.ndarray: - """Compute the multiplicative output uncertainty residual at a frequency. + """Compute the multiplicative output uncertainty residual frequency responses. 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. + complex_nominal : np.ndarray + Nominal model frequency response. + complex_offnominal : np.ndarray + Off-nominal model frequency responses. tol_residual_existence : float - Tolerance for the existence of a multiplicative input uncertainty residual. + Tolerance for the existence of an uncertainty residual. Returns ------- np.ndarray - Multiplicative output uncertainty residual at a given frequency. + Multiplicative output uncertainty residual frequency response for all + off-nominal models. Raises ------ @@ -357,19 +238,19 @@ def _compute_uncertainty_residual_multiplicative_output_freq( solution. """ - num_inputs = complex_response_nom_freq.shape[1] - num_outputs = complex_response_nom_freq.shape[0] + nbr_inputs = complex_nominal.shape[-1] + nbr_outputs = complex_nominal.shape[-2] - 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 + a = np.moveaxis(complex_nominal, -1, -2) + b = np.moveaxis(complex_offnominal - complex_nominal, -1, -2) + x, residues_lstsq, _, _ = scipy.linalg.lstsq(a, b) + complex_residual = np.moveaxis(x, -1, -2) - if num_inputs <= num_outputs: - return complex_response_residual_freq + if nbr_inputs <= nbr_outputs: + return complex_residual else: if np.all(residues_lstsq <= tol_residual_existence): - return complex_response_residual_freq + return complex_residual else: raise ValueError( "A multiplicative output uncertainty residual does not exist for the " @@ -381,27 +262,27 @@ def _compute_uncertainty_residual_multiplicative_output_freq( ) -def _compute_uncertainty_residual_inverse_additive_freq( - complex_response_nom_freq: np.ndarray, - complex_response_offnom_freq: np.ndarray, +def _compute_residual_inverse_additive( + complex_nominal: np.ndarray, + complex_offnominal: np.ndarray, tol_residual_existence: float = 1e-8, ) -> np.ndarray: - """Compute the inverse additive uncertainty residual at a frequency. + """Compute the inverse additive uncertainty residual frequency responses. 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. + complex_nominal : np.ndarray + Nominal model frequency response. + complex_offnominal : np.ndarray + Off-nominal model frequency responses. tol_residual_existence : float - Tolerance for the existence of a inverse additive uncertainty residual. + Tolerance for the existence of an uncertainty residual. Returns ------- np.ndarray - Inverse additive uncertainty residual at a given frequency. + Inverse additive uncertainty residual frequency response for all + off-nominal models. Raises ------ @@ -411,24 +292,24 @@ def _compute_uncertainty_residual_inverse_additive_freq( solution. """ - num_inputs = complex_response_nom_freq.shape[1] - num_outputs = complex_response_nom_freq.shape[0] + nbr_inputs = complex_nominal.shape[-1] + nbr_outputs = complex_nominal.shape[-2] - 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 + a1 = complex_offnominal + b1 = complex_offnominal - complex_nominal + y, residues_lstsq_1, _, _ = scipy.linalg.lstsq(a1, b1) + a2 = np.moveaxis(complex_nominal, -1, -2) + b2 = np.moveaxis(y, -1, -2) + x, residues_lstsq_2, _, _ = scipy.linalg.lstsq(a2, b2) + complex_residual = np.moveaxis(x, -1, -2) - if num_inputs == num_outputs: - return complex_response_residual_freq + if nbr_inputs == nbr_outputs: + return complex_residual else: if np.all(residues_lstsq_1 <= tol_residual_existence) and np.all( residues_lstsq_2 <= tol_residual_existence ): - return complex_response_residual_freq + return complex_residual else: raise ValueError( "An inverse additive uncertainty residual does not exist for the " @@ -440,28 +321,27 @@ def _compute_uncertainty_residual_inverse_additive_freq( ) -def _compute_uncertainty_residual_inverse_multiplicative_input_freq( - complex_response_nom_freq: np.ndarray, - complex_response_offnom_freq: np.ndarray, +def _compute_residual_inverse_multiplicative_input( + complex_nominal: np.ndarray, + complex_offnominal: np.ndarray, tol_residual_existence: float = 1e-8, ) -> np.ndarray: - """Compute the inverse multiplicative input uncertainty residual at a frequency. + """Compute the inverse multiplicative input uncertainty residual frequency responses. 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. + complex_nominal : np.ndarray + Nominal model frequency response. + complex_offnominal : np.ndarray + Off-nominal model frequency responses. tol_residual_existence : float - Tolerance for the existence of an inverse multiplicative input uncertainty - residual. + Tolerance for the existence of an uncertainty residual. Returns ------- np.ndarray - Inverse multiplicative input uncertainty residual at a given frequency. + Inverse multiplicative input uncertainty residual frequency response for all + off-nominal models. Raises ------ @@ -471,19 +351,19 @@ def _compute_uncertainty_residual_inverse_multiplicative_input_freq( a solution. """ - num_inputs = complex_response_nom_freq.shape[1] - num_outputs = complex_response_nom_freq.shape[0] + nbr_inputs = complex_nominal.shape[-1] + nbr_outputs = complex_nominal.shape[-2] - 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 + a = complex_offnominal + b = complex_offnominal - complex_nominal + x, residues_lstsq, _, _ = scipy.linalg.lstsq(a, b) + complex_residual = x - if num_inputs >= num_outputs: - return complex_response_residual_freq + if nbr_inputs >= nbr_outputs: + return complex_residual else: if np.all(residues_lstsq <= tol_residual_existence): - return complex_response_residual_freq + return complex_residual else: raise ValueError( "An inverse multiplicative input uncertainty residual does not exist " @@ -495,27 +375,27 @@ def _compute_uncertainty_residual_inverse_multiplicative_input_freq( ) -def _compute_uncertainty_residual_inverse_multiplicative_output_freq( - complex_response_nom_freq: np.ndarray, - complex_response_offnom_freq: np.ndarray, +def _compute_residual_inverse_multiplicative_output( + complex_nominal: np.ndarray, + complex_offnominal: np.ndarray, tol_residual_existence: float = 1e-8, ) -> np.ndarray: - """Compute the inverse multiplicative output uncertainty residual at a frequency. + """Compute the inverse multiplicative output uncertainty residual frequency responses. 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. + complex_nominal : np.ndarray + Nominal model frequency response. + complex_offnominal : np.ndarray + Off-nominal model frequency responses. tol_residual_existence : float - Tolerance for the existence of a multiplicative input uncertainty residual. + Tolerance for the existence of an uncertainty residual. Returns ------- np.ndarray - Inverse multiplicative output uncertainty residual at a given frequency. + Inverse multiplicative output uncertainty residual frequency response for all + off-nominal models. Raises ------ @@ -525,19 +405,19 @@ def _compute_uncertainty_residual_inverse_multiplicative_output_freq( a solution. """ - num_inputs = complex_response_nom_freq.shape[1] - num_outputs = complex_response_nom_freq.shape[0] + nbr_inputs = complex_nominal.shape[-1] + nbr_outputs = complex_nominal.shape[-2] - 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 + a = np.moveaxis(complex_offnominal, -1, -2) + b = np.moveaxis(complex_offnominal - complex_nominal, -1, -2) + x, residues_lstsq, _, _ = scipy.linalg.lstsq(a, b) + complex_residual = np.moveaxis(x, -1, -2) - if num_inputs <= num_outputs: - return complex_response_residual_freq + if nbr_inputs <= nbr_outputs: + return complex_residual else: if np.all(residues_lstsq <= tol_residual_existence): - return complex_response_residual_freq + return complex_residual else: raise ValueError( "An inverse multiplicative output uncertainty residual does not exist " @@ -550,9 +430,9 @@ 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, + complex_residual: Union[np.ndarray, control.FrequencyResponseList], + 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. @@ -562,16 +442,14 @@ def compute_uncertainty_weight_response( Parameters ---------- - complex_response_residual_list : Union[np.ndarray, control.FrequencyResponseList] + complex_residual : 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] + weight_left_structure : Literal["full", "diagonal", "scalar", "identity"] + Structure of the left uncertainty weight. + weight_right_structure : Literal["full", "diagonal", "scalar", "identity"] + Structure of the right uncertainty weight. + solver_params : Dict[str, Any] Keyword arguments for the convex optimization solver. See [#cvxpy_solver]_ for more information. @@ -587,25 +465,25 @@ def compute_uncertainty_weight_response( residual frequency response where the left and right uncertainty weights are assumed to be diagonal. - >>> complex_response_nom, complex_response_offnom_list, omega = ( + >>> complex_nominal, complex_offnominal_list, omega = ( ... example_multimodel_uncertainty ... ) - >>> uncertainty_models = { + >>> 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, + ... ] + >>> complex_residual_dict = compute_uncertainty_residual_response( + ... complex_nominal, + ... complex_offnominal_list, ... uncertainty_models, ... ) - >>> complex_response_weight_left, complex_response_weight_right = ( + >>> complex_weight_left, complex_weight_right = ( ... dkpy.compute_uncertainty_weight_response( - ... complex_response_residual_dict["multiplicative_input"], + ... complex_residual_dict["multiplicative_input"], ... "diagonal", ... "diagonal", ... ) @@ -616,72 +494,12 @@ def compute_uncertainty_weight_response( .. [#uncertainty_characterization] G. J. Balas, A. K. Packard, and P. J. Seiler, “Uncertain Model Set Calculation from Frequency Domain Data,” Springer eBooks, pp. 89–105, Jan. 2009, doi: https://doi.org/10.1007/978-1-4419-0895-7_6. - """ # 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) + complex_residual = _convert_frequency_response_list_to_array(complex_residual) - # 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 + # Parse solver parameters solver_params = ( { "solver": cvxpy.CLARABEL, @@ -695,130 +513,148 @@ def _compute_optimal_weight_freq( 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] + # Frequency response parameters + nbr_frequency = complex_residual.shape[1] + nbr_left = complex_residual.shape[2] + nbr_right = complex_residual.shape[3] + nbr_offnom = complex_residual.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".' + weight_left_power_dispatcher = { + "full": cvxpy.Variable((nbr_left, nbr_left), hermitian=True), + "diagonal": cvxpy.Variable((nbr_left, nbr_left), diag=True), + "scalar": cvxpy.Variable() * scipy.sparse.eye_array(nbr_left), + "identity": cvxpy.Constant(value=np.eye(nbr_left)), + } + try: + weight_left_power = weight_left_power_dispatcher[weight_left_structure] + except KeyError: + raise KeyError( + '`weight_left_structure` must be "full", "diagonal", or "scalar" (got ' + f'"{weight_left_structure})".' ) # 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".' + weight_right_power_dispatcher = { + "full": cvxpy.Variable((nbr_right, nbr_right), hermitian=True), + "diagonal": cvxpy.Variable((nbr_right, nbr_right), diag=True), + "scalar": cvxpy.Variable() * scipy.sparse.eye_array(nbr_right), + "identity": cvxpy.Constant(value=np.eye(nbr_right)), + } + try: + weight_right_power = weight_right_power_dispatcher[weight_right_structure] + except KeyError: + raise KeyError( + '`weight_right_structure` must be "full", "diagonal", or "scalar" (got ' + f'"{weight_right_structure})".' ) - # 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, - ) + # Generate residual parameters + residual_offnom = cvxpy.Parameter( + shape=(nbr_offnom, nbr_left, nbr_right), + complex=True, + ) + + # Uncertainty set constraints over all frequencies + constraint_list = [] + for idx_offnom in range(nbr_offnom): constraint_matrix_freq = cvxpy.bmat( [ - [L, E_k], - [E_k.H, R], + [weight_left_power, residual_offnom[idx_offnom, :, :]], + [residual_offnom[idx_offnom, :, :].H, weight_right_power], ] ) 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) + constraint_list.append(constraint_freq) + + # Positive semidefiniteness constraints + constraint_list.append(weight_left_power >> 0) + constraint_list.append(weight_right_power >> 0) # Semidefinite program - objective = cvxpy.Minimize(cvxpy.trace(L) + cvxpy.trace(R)) - problem = cvxpy.Problem(objective, constraint_freq_list) - problem.solve(**solver_params) + objective = cvxpy.Minimize( + nbr_right * cvxpy.trace(weight_left_power) + + nbr_left * cvxpy.trace(weight_right_power) + ) + problem = cvxpy.Problem(objective, constraint_list) - # 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) + # Compute optimal uncertainty weights + complex_weight_left = [] + complex_weight_right = [] + for idx_freq in range(nbr_frequency): + # Residual of off-nominal models at given frequency + residual_offnom.value = complex_residual[:, idx_freq, :, :] + + # Solve optimal weight SDP + problem.solve(canon_backend=cvxpy.SCIPY_CANON_BACKEND, **solver_params) + + # Extract left weight + if weight_left_structure == "identity" or weight_left_structure == "full": + weight_left_power_opt = np.array(weight_left_power.value) + else: + weight_left_power_opt = np.array(weight_left_power.value.toarray()) + weight_left_opt = scipy.linalg.sqrtm(weight_left_power_opt) + + # Extract right weight + if weight_right_structure == "identity" or weight_right_structure == "full": + weight_right_power_opt = np.array(weight_right_power.value) + else: + weight_right_power_opt = np.array(weight_right_power.value.toarray()) + weight_right_opt = scipy.linalg.sqrtm(weight_right_power_opt) + + complex_weight_left.append(weight_left_opt) + complex_weight_right.append(weight_right_opt) + + # Generate uncertainty weight complex frequency reponses + complex_weight_left = np.array(complex_weight_left) + complex_weight_right = np.array(complex_weight_right) - return complex_response_weight_left_freq, complex_response_weight_right_freq + return complex_weight_left, complex_weight_right def fit_uncertainty_weight( - complex_response_uncertainty_weight: Union[ - np.ndarray, control.FrequencyResponseData - ], + complex_uncertainty_weight: Union[np.ndarray, control.FrequencyResponseData], omega: np.ndarray, order: Union[int, List[int], np.ndarray], + uncertainty_weight_type: Literal["left", "right"], + uncertainty_weight_structure: Literal["scalar", "diagonal", "full"], weight: Optional[np.ndarray] = None, - linear_solver_params: Optional[Dict[str, Any]] = None, + solver_params: Optional[Dict[str, Any]] = None, tol_bisection: float = 1e-3, max_iter_bisection: int = 500, - num_spec_constr: int = 500, + max_iter_bisection_init: int = 15, + nbr_power_constraint: int = 500, ) -> control.StateSpace: - """Fit an overbounding stable and minimum-phase state-space uncertainty weight to - frequency response data. + """Fit an overbounding stable and minimum-phase uncertainty weight. Parameters ---------- - complex_response_uncertainty_weight : Union[np.ndarray, control.FrequencyResponseData] - Uncertainty weight frequency response used for the overbounding fit. + complex_uncertainty_weight : Union[np.ndarray, control.FrequencyResponseData] + Frequency response of uncertainty weight. 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. + Order of the uncertainty weight model. + uncertainty_weight_type: Literal["left", "right"], + Identifier for the left or right uncertainty weight. + uncertainty_weight_structure : Literal["scalar", "diagonal", "full"] + Structure constraint for the uncertainty weight. + weight : Optional[np.ndrray] = None + Frequency-dependent weight for fit accuracy. + solver_params: Optional[Dict[str, Any]] + Solver parameters for the optimization problem. These are keyword arguments for + `cvxpy.Problem.solve()` [#cvxpy_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. + Maximum number of iterations for the bisection algorithm. + max_iter_bisection_init : int + Maximum number of iterations for the bisection algorithm initialization. + nbr_power_constraint : Optional[np.ndarray] + Number of frequencies to enforce non-negativity of power spectrum. Returns ------- control.StateSpace - Fitted overbounding uncertainty weight state-space systems. + Overbounding stable and minimum-phase uncertainty weight. Examples -------- @@ -828,14 +664,7 @@ def fit_uncertainty_weight( >>> 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", - ... } + >>> uncertainty_models = ["multiplicative_input"] >>> complex_response_residual_dict = compute_uncertainty_residual_response( ... complex_response_nom, ... complex_response_offnom_list, @@ -849,10 +678,10 @@ def fit_uncertainty_weight( ... ) ... ) >>> weight_left = dkpy.fit_uncertainty_weight( - ... complex_response_weight_left, omega, [4, 5] + ... complex_response_weight_left, omega, [4, 5], "left", "diagonal" ... ) >>> weight_right = dkpy.fit_uncertainty_weight( - ... complex_response_weight_right, omega, [3, 5] + ... complex_response_weight_right, omega, [3, 5], "right", "diagonal" ... ) References @@ -861,12 +690,12 @@ def fit_uncertainty_weight( """ # Convert frequency response data to expected type - complex_response_uncertainty_weight = _convert_frequency_response_data_to_array( - complex_response_uncertainty_weight + complex_uncertainty_weight = _convert_frequency_response_data_to_array( + complex_uncertainty_weight ) # Solver settings - linear_solver_params = ( + solver_params = ( { "solver": cvxpy.CLARABEL, "tol_gap_abs": 1e-9, @@ -875,56 +704,748 @@ def fit_uncertainty_weight( "tol_infeas_abs": 1e-9, "tol_infeas_rel": 1e-9, } - if linear_solver_params is None - else linear_solver_params + if solver_params is None + else 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) - ) + # Fit method dispatcher + fit_method_dispatcher = { + "scalar": _fit_uncertainty_weight_scalar, + "diagonal": _fit_uncertainty_weight_diagonal, + "full": _fit_uncertainty_weight_full, + } - 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 + # Select uncertainty weight fit method based on weight structure + try: + fit_uncertainty_weight_method = fit_method_dispatcher[ + uncertainty_weight_structure + ] + except KeyError: + raise KeyError( + "Exptected `weight_structure` to be `scalar`, `diagonal` or `full` (got " + f"`{uncertainty_weight_structure}`)." ) - 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 uncertainty weight + uncertainty_weight = fit_uncertainty_weight_method( + complex_uncertainty_weight, + omega, + order, + uncertainty_weight_type, + weight, + solver_params, + tol_bisection, + max_iter_bisection, + max_iter_bisection_init, + nbr_power_constraint, + ) + + return uncertainty_weight + + +def _fit_uncertainty_weight_scalar( + complex_uncertainty_weight: np.ndarray, + omega: np.ndarray, + order: int, + uncertainty_weight_type: Literal["left", "right"], + weight: Optional[np.ndarray] = None, + solver_params: Optional[Dict[str, Any]] = None, + tol_bisection: float = 1e-3, + max_iter_bisection: int = 500, + max_iter_bisection_init: int = 15, + nbr_power_constraint: int = 500, +) -> control.StateSpace: + """Fit an overbounding stable and minimum-phase scalar uncertainty weight. + + Parameters + ---------- + complex_uncertainty_weight : Union[np.ndarray, control.FrequencyResponseData] + Frequency response of uncertainty weight. + order : Union[int, List[int], np.ndarray] + Order of the uncertainty weight model. + uncertainty_weight_type: Literal["left", "right"], + Identifier for the left or right uncertainty weight. + weight : Optional[np.ndrray] = None + Frequency-dependent weight for fit accuracy. + solver_params: Optional[Dict[str, Any]] + Solver parameters for the optimization problem. These are keyword arguments for + `cvxpy.Problem.solve()` [#cvxpy_solver]_. + tol_bisection : float + Numerical tolerance for the bisection algorithm. + max_iter_bisection : int + Maximum number of iterations for the bisection algorithm. + max_iter_bisection_init : int + Maximum number of iterations for the bisection algorithm initialization. + nbr_power_constraint : Optional[np.ndarray] + Number of frequencies to enforce non-negativity of power spectrum. + + Returns + ------- + control.StateSpace + Overbounding stable and minimum-phase uncertainty weight. + + References + ---------- + .. [#cxvpy_solver] https://www.cvxpy.org/tutorial/solvers/index.html + """ + + # Auxiliary parameters + nbr_signals = complex_uncertainty_weight.shape[1] + + # Compute the magnitude of the scalar weight. Given that a scalar weight is assumed, + # the weight contains the same values on the diagonals and zeros elsewhere. + # Therefore, we take the first diagonal as they are identical along the diagonal. + magnitude_uncertainty_weight = np.abs(complex_uncertainty_weight[:, 0, 0]) + + # Fit the scalar uncertainty weight + uncertainty_weight_scalar = lti_system_fit.fit_magnitude_siso_ct( + magnitude_fit=magnitude_uncertainty_weight, + omega=omega, + order=order, + magnitude_upper_bound=None, + magnitude_lower_bound=magnitude_uncertainty_weight, + weight=weight, + solver_params=solver_params, + tol_bisection=tol_bisection, + max_iter_bisection=max_iter_bisection, + max_iter_bisection_init=max_iter_bisection_init, + nbr_power_constraint=nbr_power_constraint, + ) + + # Duplicate the scalar uncertainty weight along the diagonal + uncertainty_weight = control.append(*([uncertainty_weight_scalar] * nbr_signals)) + + return uncertainty_weight + + +def _fit_uncertainty_weight_diagonal( + complex_uncertainty_weight: np.ndarray, + omega: np.ndarray, + order: Union[int, List[int]], + uncertainty_weight_type: Literal["left", "right"], + weight: Optional[np.ndarray] = None, + solver_params: Optional[Dict[str, Any]] = None, + tol_bisection: float = 1e-3, + max_iter_bisection: int = 500, + max_iter_bisection_init: int = 15, + nbr_power_constraint: int = 500, +): + """Fit an overbounding stable and minimum-phase diagonal uncertainty weight. + + Parameters + ---------- + complex_uncertainty_weight : Union[np.ndarray, control.FrequencyResponseData] + Frequency response of uncertainty weight. + order : Union[int, List[int], np.ndarray] + Order of the uncertainty weight model. + uncertainty_weight_type: Literal["left", "right"], + Identifier for the left or right uncertainty weight. + weight : Optional[np.ndrray] = None + Frequency-dependent weight for fit accuracy. + solver_params: Optional[Dict[str, Any]] + Solver parameters for the optimization problem. These are keyword arguments for + `cvxpy.Problem.solve()` [#cvxpy_solver]_. + tol_bisection : float + Numerical tolerance for the bisection algorithm. + max_iter_bisection : int + Maximum number of iterations for the bisection algorithm. + max_iter_bisection_init : int + Maximum number of iterations for the bisection algorithm initialization. + nbr_power_constraint : Optional[np.ndarray] + Number of frequencies to enforce non-negativity of power spectrum. - # Fit the uncertainty weight to each SISO element - uncertainty_weight_element = utilities._fit_magnitude_log_chebyshev_siso( + Returns + ------- + control.StateSpace + Overbounding stable and minimum-phase uncertainty weight. + + References + ---------- + .. [#cxvpy_solver] https://www.cvxpy.org/tutorial/solvers/index.html + """ + # Auxiliary parameters + nbr_signals = complex_uncertainty_weight.shape[1] + + # Parse order into list format + order_list = [order] * nbr_signals if isinstance(order, int) else order + + # Fit the diagonal uncertainty weight elements + uncertainty_weight_list = [] + 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]) + order_idx = order_list[idx] + + # Fit the diagonal uncertainty weight element + uncertainty_weight_diagonal = lti_system_fit.fit_magnitude_siso_ct( + magnitude_fit=magnitude_uncertainty_weight, 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, + order=order_idx, + magnitude_upper_bound=None, + magnitude_lower_bound=magnitude_uncertainty_weight, + weight=weight, + solver_params=solver_params, tol_bisection=tol_bisection, max_iter_bisection=max_iter_bisection, - num_spec_constr=num_spec_constr, + max_iter_bisection_init=max_iter_bisection_init, + nbr_power_constraint=nbr_power_constraint, ) - uncertainty_weight_list.append(uncertainty_weight_element) + uncertainty_weight_list.append(uncertainty_weight_diagonal) - # Construct uncertainty weight fit from SISO elements + # Construct the diagonal uncertainty weight from the diagonal elements uncertainty_weight = control.append(*uncertainty_weight_list) return uncertainty_weight +def _fit_uncertainty_weight_full( + complex_uncertainty_weight: np.ndarray, + omega: np.ndarray, + order: int, + uncertainty_weight_type: Literal["left", "right"], + weight: Optional[np.ndarray] = None, + solver_params: Optional[Dict[str, Any]] = None, + tol_bisection: float = 1e-3, + max_iter_bisection: int = 500, + max_iter_bisection_init: int = 15, + nbr_power_constraint: int = 500, +): + """Fit an overbounding stable and minimum-phase full uncertainty weight. + + Parameters + ---------- + complex_uncertainty_weight : Union[np.ndarray, control.FrequencyResponseData] + Frequency response of uncertainty weight. + order : Union[int, List[int], np.ndarray] + Order of the uncertainty weight model. + uncertainty_weight_type: Literal["left", "right"], + Identifier for the left or right uncertainty weight. + weight : Optional[np.ndrray] = None + Frequency-dependent weight for fit accuracy. + solver_params: Optional[Dict[str, Any]] + Solver parameters for the optimization problem. These are keyword arguments for + `cvxpy.Problem.solve()` [#cvxpy_solver]_. + tol_bisection : float + Numerical tolerance for the bisection algorithm. + max_iter_bisection : int + Maximum number of iterations for the bisection algorithm. + max_iter_bisection_init : int + Maximum number of iterations for the bisection algorithm initialization. + nbr_power_constraint : Optional[np.ndarray] + Number of frequencies to enforce non-negativity of power spectrum. + + Returns + ------- + control.StateSpace + Overbounding stable and minimum-phase uncertainty weight. + + References + ---------- + .. [#cxvpy_solver] https://www.cvxpy.org/tutorial/solvers/index.html + """ + + raise NotImplementedError() + + +def compute_uncertainty_measure_response( + complex_nominal: Union[np.ndarray, control.FrequencyResponseData], + complex_weight_left: Union[np.ndarray, control.FrequencyResponseData], + complex_weight_right: Union[np.ndarray, control.FrequencyResponseData], + uncertainty_model: str, +) -> np.ndarray: + """Compute measure (size/volume) frequency response of uncertainty model. + + Parameters + ---------- + complex_nominal : Union[np.ndarray, control.FrequencyResponseData] + Nominal model complex frequency response. + complex_weight_left : Union[np.ndarray, control.FrequencyResponseData] + Left uncertainty weight complex frequency response. + complex_weight_right : Union[np.ndarray, control.FrequencyResponseData] + Right uncertainty weight complex frequency reponse. + uncertainty_model : str + Uncertainty model identifier to compute the residual response. The valid + uncertainty model identifiers are: + - "additive" + - "multiplicative_input" + - "multiplicative_output" + - "inverse_additive" + - "inverse_multiplicative_input", + - "inverse_multiplicative_output" + + + Returns + ------- + np.ndarray + Measure frequency response. + + + Examples + -------- + Compute the uncertainty measure response for a multiplicative input uncertainty + model. + + >>> complex_nominal, complex_offnominal_list, omega = ( + ... example_multimodel_uncertainty + ... ) + >>> uncertainty_models = ["multiplicative_input"] + >>> complex_residual_dict = compute_uncertainty_residual_response( + ... complex_nominal, + ... complex_offnominal_list, + ... uncertainty_models, + ... ) + >>> complex_weight_left, complex_weight_right = ( + ... dkpy.compute_uncertainty_weight_response( + ... complex_residual_dict["multiplicative_input"], + ... "diagonal", + ... "diagonal", + ... ) + ... ) + >>> measure = dkpy.compute_uncertainty_measure_response( + ... complex_nominal, + ... complex_weight_left, + ... complex_weight_right, + ... "multiplicative_input", + ... ) + """ + + complex_nominal = _convert_frequency_response_data_to_array(complex_nominal) + complex_weight_left = _convert_frequency_response_data_to_array(complex_weight_left) + complex_weight_right = _convert_frequency_response_data_to_array( + complex_weight_right + ) + + transform_to_additive_dispatcher = { + "additive": _transform_additive_to_additive, + "multiplicative_input": _transform_multiplicative_input_to_additive, + "multiplicative_output": _transform_multiplicative_output_to_additive, + "inverse_additive": _transform_inverse_additive_to_additive, + "inverse_multiplicative_input": _transform_inverse_multiplicative_input_to_additive, + "inverse_multiplicative_output": _transform_inverse_multiplicative_output_to_additive, + } + + # Compute equivalent additive uncertainty model frequency response + complex_uncertainty_add = transform_to_additive_dispatcher[uncertainty_model]( + complex_nominal, complex_weight_left, complex_weight_right + ) + complex_nominal_add = complex_uncertainty_add[0] + complex_weight_left_add = complex_uncertainty_add[1] + complex_weight_right_add = complex_uncertainty_add[2] + + # Compute minimal representation of additive uncertainty frequency response + complex_weight_add_min = _compute_minimal_weights_additive( + complex_weight_left_add, complex_weight_right_add + ) + complex_weight_left_add = complex_weight_add_min[0] + complex_weight_right_add = complex_weight_add_min[1] + + # Compute measure frequency response of minimal additive uncertainty response + measure = _compute_uncertainty_measure_additive( + complex_weight_left_add, complex_weight_right_add + ) + + return measure + + +def _compute_uncertainty_measure_additive( + complex_weight_left: np.ndarray, + complex_weight_right: np.ndarray, +) -> np.ndarray: + """Compute measure (size/volume) frequency response of additive uncertainty. + + Parameters + ---------- + complex_nominal : np.ndarray + Nominal model complex frequency response of additive uncertainty. + complex_weight_left : np.ndarray + Left uncertainty weight complex frequency response of additive uncertainty. + complex_weight_right : np.ndarray + Right uncertainty weight complex frequency response of additive uncertainty + + Returns + ------- + np.ndarray + Measure frequency response of additive uncertainty. + + Raises + ------ + ValueError + Left uncertainty weight is non-square. + ValueError + Right uncertainty weight is non-square. + """ + + # Check uncertainty weight dimensions + if complex_weight_left.shape[1] != complex_weight_left.shape[2]: + raise ValueError( + "Left uncertainty weight must be square (got " + f"{complex_weight_left.shape[1]} rows and {complex_weight_left.shape[2]} " + f"columns)." + ) + if complex_weight_right.shape[1] != complex_weight_right.shape[2]: + raise ValueError( + "Right uncertainty weight must be square (got " + f"{complex_weight_right.shape[1]} rows and {complex_weight_right.shape[2]} " + f"columns)." + ) + + # Auxiliary parameters + nbr_inputs = complex_weight_right.shape[1] + nbr_outputs = complex_weight_left.shape[1] + + # Measure computation + measure_weight_left = np.abs(np.linalg.det(complex_weight_left)) ** ( + 2 / nbr_outputs + ) + measure_weight_right = np.abs(np.linalg.det(complex_weight_right)) ** ( + 2 / nbr_inputs + ) + measure = measure_weight_left * measure_weight_right + + return measure + + +def _transform_additive_to_additive( + complex_nominal: np.ndarray, + complex_weight_left: np.ndarray, + complex_weight_right: np.ndarray, +) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """Transform additive uncertainty frequency response to additive form. + + Parameters + ---------- + complex_nominal : np.ndarray + Nominal model complex frequency response of additive uncertainty. + complex_weight_left : np.ndarray + Left uncertainty weight complex frequency response of additive uncertainty. + complex_weight_right : np.ndarray + Right uncertainty weight complex frequency response of additive uncertainty. + + Returns + ------- + Tuple[np.ndarray, np.ndarray, np.ndarray] + Nominal model, left uncertainty weight, and right uncertainty weight complex + frequency response of equivalent additive uncertainty. + """ + + # Nominal model + complex_nominal_add = complex_nominal + + # Left uncertainty weight + complex_weight_left_add = complex_weight_left + + # Right uncertainty weight + complex_weight_right_add = complex_weight_right + + return complex_nominal_add, complex_weight_left_add, complex_weight_right_add + + +def _transform_multiplicative_input_to_additive( + complex_nominal: np.ndarray, + complex_weight_left: np.ndarray, + complex_weight_right: np.ndarray, +) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """Transform multiplicative input uncertainty frequency response to additive form. + + Parameters + ---------- + complex_nominal : np.ndarray + Nominal model complex frequency response of multiplicative input uncertainty. + complex_weight_left : np.ndarray + Left uncertainty weight complex frequency response of multiplicative input + uncertainty. + complex_weight_right : np.ndarray + Right uncertainty weight complex frequency response of multiplicative input + uncertainty. + + Returns + ------- + Tuple[np.ndarray, np.ndarray, np.ndarray] + Nominal model, left uncertainty weight, and right uncertainty weight complex + frequency response of equivalent additive uncertainty. + """ + + # Nominal model + complexnominal_add = complex_nominal + + # Left uncertainty weight + complex_weight_left_add = complex_nominal @ complex_weight_left + + # Right uncertainty weight + complex_weight_right_add = complex_weight_right + + return complexnominal_add, complex_weight_left_add, complex_weight_right_add + + +def _transform_multiplicative_output_to_additive( + complex_nominal: np.ndarray, + complex_weight_left: np.ndarray, + complex_weight_right: np.ndarray, +) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """Transform multiplicative output uncertainty frequency response to additive form. + + Parameters + ---------- + complex_nominal : np.ndarray + Nominal model complex frequency response of multiplicative output uncertainty. + complex_weight_left : np.ndarray + Left uncertainty weight complex frequency response of multiplicative output + uncertainty. + complex_weight_right : np.ndarray + Right uncertainty weight complex frequency response of multiplicative output + uncertainty. + + Returns + ------- + Tuple[np.ndarray, np.ndarray, np.ndarray] + Nominal model, left uncertainty weight, and right uncertainty weight complex + frequency response of equivalent additive uncertainty. + """ + + # Nominal model + complex_nominal_add = complex_nominal + + # Left uncertainty weight + complex_weight_left_add = complex_weight_left + + # Right uncertainty weight + complex_weight_right_add = complex_weight_right @ complex_nominal + + return complex_nominal_add, complex_weight_left_add, complex_weight_right_add + + +def _transform_inverse_additive_to_additive( + complex_nominal: np.ndarray, + complex_weight_left: np.ndarray, + complex_weight_right: np.ndarray, +) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """Transform inverse additive uncertainty frequency response to additive form. + + Parameters + ---------- + complex_nominal : np.ndarray + Nominal model complex frequency response of inverse additive uncertainty. + complex_weight_left : np.ndarray + Left uncertainty weight complex frequency response of inverse additive + uncertainty. + complex_weight_right : np.ndarray + Right uncertainty weight complex frequency response of inverse additive + uncertainty. + + Returns + ------- + Tuple[np.ndarray, np.ndarray, np.ndarray] + Nominal model, left uncertainty weight, and right uncertainty weight complex + frequency response of equivalent additive uncertainty. + """ + + # Auxiliary variables + product = complex_weight_right @ complex_nominal @ complex_weight_left + product_herm = np.moveaxis(product.conj(), -1, -2) + eye_outer = np.eye(product.shape[-2])[None, :, :] + eye_inner = np.eye(product.shape[-1])[None, :, :] + + # Nominal model + a_nominal = eye_outer - product @ product_herm + b_nominal = complex_weight_right @ complex_nominal + x_nominal = np.linalg.solve(a_nominal, b_nominal) + complex_nominal_add = ( + complex_nominal + + complex_nominal @ complex_weight_left @ product_herm @ x_nominal + ) + + # Left uncertainty weight + # The transpose using np.moveaxis is required as the linear system solved is in the + # form X @ A = B whereas numpy requires C @ X = D. The original problem is converted + # to the equivalent problem A.T @ X.T = B.T. + a_left = np.moveaxis(scipy.linalg.sqrtm(eye_inner - product_herm @ product), -1, -2) + b_left = np.moveaxis(-(complex_nominal @ complex_weight_left), -1, -2) + complex_weight_left_add = np.moveaxis(np.linalg.solve(a_left, b_left), -1, -2) + + # Right uncertainty weight + a_right = scipy.linalg.sqrtm(eye_outer - product @ product_herm) + b_right = complex_weight_right @ complex_nominal + complex_weight_right_add = np.linalg.solve(a_right, b_right) + + return complex_nominal_add, complex_weight_left_add, complex_weight_right_add + + +def _transform_inverse_multiplicative_input_to_additive( + complex_nominal: np.ndarray, + complex_weight_left: np.ndarray, + complex_weight_right: np.ndarray, +) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """Transform inverse multiplicative input uncertainty frequency response to additive form. + + Parameters + ---------- + complex_nominal : np.ndarray + Nominal model complex frequency response of inverse multiplicative input + uncertainty. + complex_weight_left : np.ndarray + Left uncertainty weight complex frequency response of inverse multiplicative + input uncertainty. + complex_weight_right : np.ndarray + Right uncertainty weight complex frequency response of inverse multiplicative + input uncertainty. + + Returns + ------- + Tuple[np.ndarray, np.ndarray, np.ndarray] + Nominal model, left uncertainty weight, and right uncertainty weight complex + frequency response of equivalent additive uncertainty. + """ + # Auxiliary variables + product = complex_weight_right @ complex_weight_left + product_herm = np.moveaxis(product.conj(), -1, -2) + eye_outer = np.eye(product.shape[-2])[None, :, :] + eye_inner = np.eye(product.shape[-1])[None, :, :] + + # Nominal model + a_nominal = eye_outer - product @ product_herm + b_nominal = complex_weight_right + x_nominal = np.linalg.solve(a_nominal, b_nominal) + complex_nominal_add = ( + complex_nominal + + complex_nominal @ complex_weight_left @ product_herm @ x_nominal + ) + + # Left uncertainty weight + # The transpose using np.moveaxis is required as the linear system solved is in the + # form X @ A = B whereas numpy requires C @ X = D. The original problem is converted + # to the equivalent problem A.T @ X.T = B.T. + a_left = np.moveaxis(scipy.linalg.sqrtm(eye_inner - product_herm @ product), -1, -2) + b_left = np.moveaxis(-(complex_nominal @ complex_weight_left), -1, -2) + complex_weight_left_add = np.moveaxis(np.linalg.solve(a_left, b_left), -1, -2) + + # Right uncertainty weight + a_right = scipy.linalg.sqrtm(eye_outer - product @ product_herm) + b_right = complex_weight_right + complex_weight_right_add = np.linalg.solve(a_right, b_right) + + return complex_nominal_add, complex_weight_left_add, complex_weight_right_add + + +def _transform_inverse_multiplicative_output_to_additive( + complex_nominal: np.ndarray, + complex_weight_left: np.ndarray, + complex_weight_right: np.ndarray, +) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """Transform inverse multiplicative output uncertainty frequency response to additive form. + + Parameters + ---------- + complex_nominal : np.ndarray + Nominal model complex frequency response of inverse multiplicative output + uncertainty. + complex_weight_left : np.ndarray + Left uncertainty weight complex frequency response of inverse multiplicative + output uncertainty. + complex_weight_right : np.ndarray + Right uncertainty weight complex frequency response of inverse multiplicative + output uncertainty. + + Returns + ------- + Tuple[np.ndarray, np.ndarray, np.ndarray] + Nominal model, left uncertainty weight, and right uncertainty weight complex + frequency response of equivalent additive uncertainty. + """ + # Auxiliary variables + product = complex_weight_right @ complex_weight_left + product_herm = np.moveaxis(product.conj(), -1, -2) + eye_outer = np.eye(product.shape[-2])[None, :, :] + eye_inner = np.eye(product.shape[-1])[None, :, :] + + # Nominal model + a_nominal = eye_outer - product @ product_herm + b_nominal = complex_weight_right @ complex_nominal + x_nominal = np.linalg.solve(a_nominal, b_nominal) + complex_nominal_add = ( + complex_nominal + complex_weight_left @ product_herm @ x_nominal + ) + + # Left uncertainty weight + # The transpose using np.moveaxis is required as the linear system solved is in the + # form X @ A = B whereas numpy requires C @ X = D. The original problem is converted + # to the equivalent problem A.T @ X.T = B.T. + a_left = np.moveaxis(scipy.linalg.sqrtm(eye_inner - product_herm @ product), -1, -2) + b_left = np.moveaxis(-complex_weight_left, -1, -2) + complex_weight_left_add = np.moveaxis(np.linalg.solve(a_left, b_left), -1, -2) + + # Right uncertainty weight + a_right = scipy.linalg.sqrtm(eye_outer - product @ product_herm) + b_right = complex_weight_right @ complex_nominal + complex_weight_right_add = np.linalg.solve(a_right, b_right) + + return complex_nominal_add, complex_weight_left_add, complex_weight_right_add + + +def _compute_minimal_weights_additive( + complex_weight_left: np.ndarray, + complex_weight_right: np.ndarray, +) -> Tuple[np.ndarray, np.ndarray]: + """Compute minimal dimension additive uncertainty weight frequency response. + + Parameters + ---------- + complex_weight_left : np.ndarray + Left uncertainty weight complex response. + complex_weight_right : np.ndarray + Right uncertainty weight complex response. + + Returns + ------- + Tuple[np.ndarray, np.ndarray] + Minimal left and right uncertainty weight frequency response. + + Raises + ------ + ValueError + The number of perturbation outputs is less than the number of system outputs. + ValueError + The number of perturbation inputs is less than the number of system inputs. + """ + + # Number of inputs/outputs of system and perturbation (Delta block) + nbr_outputs_sys = complex_weight_left.shape[1] + nbr_inputs_sys = complex_weight_right.shape[2] + nbr_outputs_delta = complex_weight_left.shape[2] + nbr_inputs_delta = complex_weight_right.shape[1] + + # Compute minimial left uncertainty weight + if nbr_outputs_sys > nbr_outputs_delta: + raise ValueError( + "The number of perturbation outputs (number of left uncertainty weight " + "columns) must be greater than the number of system outputs (number of " + "left uncertainty weight rows) for additive uncertainty (got " + f"{nbr_outputs_delta} perturbations outputs and {nbr_outputs_sys} system " + "outputs)." + ) + elif nbr_outputs_sys < nbr_outputs_delta: + u_left, s_left, _ = np.linalg.svd(complex_weight_left, full_matrices=False) + complex_weight_left_min = s_left[:, None, :] * u_left # Scale columns of U + else: + complex_weight_left_min = complex_weight_left + + # Compute minimal right uncertainty weight + if nbr_inputs_sys > nbr_inputs_delta: + raise ValueError( + "The number of perturbation inputs (number of right uncertainty weight " + "rows) must be greater than the number of system inputs (number of " + "right uncertainty weight columns) for additive uncertainty (got " + f"{nbr_inputs_delta} perturbations inputs and {nbr_inputs_sys} system " + "inputs)." + ) + elif nbr_inputs_sys < nbr_inputs_delta: + _, s_right, vh_right = np.linalg.svd(complex_weight_right, full_matrices=False) + complex_weight_right_min = s_right[:, :, None] * vh_right # Scale rows of V^H + else: + complex_weight_right_min = complex_weight_right + + return complex_weight_left_min, complex_weight_right_min + + def _convert_frequency_response_data_to_array( frequency_response: Union[np.typing.ArrayLike, control.FrequencyResponseData], ) -> np.ndarray: @@ -1643,12 +2164,10 @@ def plot_singular_value_response_residual_comparison( return fig, ax, legend -def plot_magnitude_response_uncertainty_weight( - complex_response_weight_left: np.ndarray, - complex_response_weight_right: np.ndarray, +def plot_singular_value_response_uncertainty_weight( + complex_weight: np.ndarray, omega: np.ndarray, - weight_left: Optional[control.StateSpace] = None, - weight_right: Optional[control.StateSpace] = None, + weight_fit: Optional[control.LTI] = None, db: bool = True, hz: bool = False, frequency_log_scale: bool = True, @@ -1656,24 +2175,17 @@ def plot_magnitude_response_uncertainty_weight( plot_response_fit_kw: Dict[str, Any] = {}, subplot_kw: Dict[str, Any] = {}, ) -> Tuple[Figure, Union[Axes, np.ndarray], Legend]: - """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. + """Plot the singular value respone of an uncertainty weight. 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 + complex_weight : 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. - 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. + weight_fit : Optional[control.LTI] + Fitted uncertainty weight model. db : bool If True, plot the magnitude in units of dB. Otherwise, plot the magnitude in absolute units. @@ -1683,12 +2195,9 @@ def plot_magnitude_response_uncertainty_weight( 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. + 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. @@ -1727,93 +2236,43 @@ def plot_magnitude_response_uncertainty_weight( 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] + # Singular value weight response + sval_weight = np.linalg.svdvals(complex_weight) + sval_weight = control.mag2db(sval_weight) if db else sval_weight - # 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) + # Singular value weight fit response + if weight_fit is None: + sval_weight_fit = None + else: + frd_weight_fit = control.FrequencyResponseData(weight_fit, omega, squeeze=False) + response_weight_fit = frd_weight_fit.complex.transpose(2, 0, 1) + sval_weight_fit = np.linalg.svdvals(response_weight_fit) + sval_weight_fit = control.mag2db(sval_weight_fit) if db else sval_weight_fit # Initialize figure - fig, ax = plt.subplots( - max(num_left, num_right), 2, sharex=True, layout="constrained" - ) + fig, ax = plt.subplots(**subplot_kwargs) - # Plot left uncertainty weight frequency response - for idx_left in range(num_left): - ax[idx_left, 0].plot( + # Plot weight singular value response + for idx_sval in range(sval_weight.shape[-1]): + ax.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], + sval_weight[:, idx_sval], **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( + if sval_weight_fit is not None: + ax.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, :], + sval_weight_fit[:, idx_sval], **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() + ax.set_xlabel("$f$ (Hz)" if hz else r"$\omega$ (rad/s)") + ax.set_ylabel("Magnitude (dB)" if db else "Magnitude (-)") + ax.grid() + if frequency_log_scale: + ax.set_xscale("log") + handles, labels = ax.get_legend_handles_labels() legend_dict = dict(zip(labels, handles)) legend = fig.legend( labels=legend_dict.keys(), @@ -1823,3 +2282,86 @@ def plot_magnitude_response_uncertainty_weight( ) return fig, ax, legend + + +def plot_uncertainty_measure( + measure: np.ndarray, + omega: np.ndarray, + db: bool = True, + hz: bool = False, + frequency_log_scale: bool = True, + plot_kw: Dict[str, Any] = {}, + subplot_kw: Dict[str, Any] = {}, +) -> Tuple[Figure, Union[Axes, np.ndarray]]: + """Plot uncertainty measure frequency response. + + Parameters + ---------- + measure : np.ndarray + Uncertainty measure frequency response. + omega : np.narray + Angular frequency grid. + weight_fit : Optional[control.LTI] + Fitted uncertainty weight model. + 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. + + Returns + ------- + Tuple[Figure, Union[Axes, np.ndarray], Legend] + Matplotlib Figure object, Axes object (or np.ndarray of Axes objects), and + Legend object. + + 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 + + """ + # Plot keyword arguments + plot_kwargs = { + "color": "C0", + "marker": "", + "linestyle": "-", + } + plot_kwargs.update(plot_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) + + # Measure response + measure = control.mag2db(measure) if db else measure + + # Initialize figure + fig, ax = plt.subplots(**subplot_kwargs) + + # Plot weight singular value response + ax.plot( + omega / (2 * np.pi) if hz else omega, + measure, + **plot_kwargs, + ) + + # Plot settings + ax.set_xlabel("$f$ (Hz)" if hz else r"$\omega$ (rad/s)") + ax.set_ylabel("Magnitude (dB)" if db else "Magnitude (-)") + ax.grid() + if frequency_log_scale: + ax.set_xscale("log") + + return fig, ax diff --git a/tests/test_uncertainty_characterization.py b/tests/test_uncertainty_characterization.py index 390053c..a2f0e25 100644 --- a/tests/test_uncertainty_characterization.py +++ b/tests/test_uncertainty_characterization.py @@ -84,21 +84,14 @@ def test_compute_uncertainty_residual_response( """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( + frd_nominal = control.frequency_response(sys_nom, omega, squeeze=False) + frd_offnominal_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 + complex_nominal = frd_nominal.complex.transpose(2, 0, 1) + complex_offnominal = np.array( + [frd.complex.transpose(2, 0, 1) for frd in frd_offnominal_list] ) # Uncertainty residual computation @@ -111,8 +104,8 @@ def test_compute_uncertainty_residual_response( "inverse_multiplicative_output", ] frequency_response_residual_dict = dkpy.compute_uncertainty_residual_response( - complex_response_nom, - complex_response_offnom_list, + complex_nominal, + complex_offnominal, uncertainty_model_list, ) ndarrays_regression.check( @@ -121,7 +114,7 @@ def test_compute_uncertainty_residual_response( ) @pytest.mark.parametrize( - "complex_response_nom_freq, complex_response_offnom_freq", + "complex_response_nom, complex_response_offnom", [ ( np.array([[1, 0], [0, 1], [1, 1]]), @@ -129,20 +122,20 @@ def test_compute_uncertainty_residual_response( ), ], ) - def test_compute_uncertainty_residual_multiplicative_input_freq( - self, complex_response_nom_freq, complex_response_offnom_freq + def test_compute_residual_multiplicative_input_error( + self, complex_response_nom, complex_response_offnom ): """Test ValueError of - :func:`_compute_uncertainty_residual_multiplicative_input_freq`. + :func:`_compute_residual_multiplicative_input`. """ with pytest.raises(ValueError): - residual_freq = dkpy.uncertainty_characterization._compute_uncertainty_residual_multiplicative_input_freq( - complex_response_nom_freq, complex_response_offnom_freq + residual_freq = dkpy.uncertainty_characterization._compute_residual_multiplicative_input( + complex_response_nom, complex_response_offnom ) @pytest.mark.parametrize( - "complex_response_nom_freq, complex_response_offnom_freq", + "complex_response_nom, complex_response_offnom", [ ( np.array([[1, 0], [0, 1], [1, 1]]).T, @@ -150,19 +143,19 @@ def test_compute_uncertainty_residual_multiplicative_input_freq( ), ], ) - def test_compute_uncertainty_residual_multiplicative_output_freq( - self, complex_response_nom_freq, complex_response_offnom_freq + def test_compute_residual_multiplicative_output_error( + self, complex_response_nom, complex_response_offnom ): """Test ValueError of - :func:`_compute_uncertainty_residual_multiplicative_output_freq`. + :func:`_compute_residual_multiplicative_output`. """ with pytest.raises(ValueError): - residual_freq = dkpy.uncertainty_characterization._compute_uncertainty_residual_multiplicative_output_freq( - complex_response_nom_freq, complex_response_offnom_freq + residual_freq = dkpy.uncertainty_characterization._compute_residual_multiplicative_output( + complex_response_nom, complex_response_offnom ) @pytest.mark.parametrize( - "complex_response_nom_freq, complex_response_offnom_freq", + "complex_response_nom, complex_response_offnom", [ ( np.array([[1, 0], [0, 1], [1, 1]]), @@ -174,19 +167,21 @@ def test_compute_uncertainty_residual_multiplicative_output_freq( ), ], ) - def test_compute_uncertainty_residual_inverse_additive_freq( - self, complex_response_nom_freq, complex_response_offnom_freq + def test_compute_residual_inverse_additive_error( + self, complex_response_nom, complex_response_offnom ): """Test ValueError of - :func:`_compute_uncertainty_residual_inverse_additive_freq`. + :func:`_compute_residual_inverse_additive`. """ with pytest.raises(ValueError): - residual_freq = dkpy.uncertainty_characterization._compute_uncertainty_residual_inverse_additive_freq( - complex_response_nom_freq, complex_response_offnom_freq + residual_freq = ( + dkpy.uncertainty_characterization._compute_residual_inverse_additive( + complex_response_nom, complex_response_offnom + ) ) @pytest.mark.parametrize( - "complex_response_nom_freq, complex_response_offnom_freq", + "complex_response_nom, complex_response_offnom", [ ( np.array([[1, 0], [0, 1], [1, 1]]), @@ -194,19 +189,19 @@ def test_compute_uncertainty_residual_inverse_additive_freq( ), ], ) - def test_compute_uncertainty_residual_inverse_multiplicative_input_freq( - self, complex_response_nom_freq, complex_response_offnom_freq + def test_compute_residual_inverse_multiplicative_input_error( + self, complex_response_nom, complex_response_offnom ): """Test ValueError of - :func:`_compute_uncertainty_residual_inverse_multiplicative_input_freq`. + :func:`_compute_residual_inverse_multiplicative_input`. """ with pytest.raises(ValueError): - residual_freq = dkpy.uncertainty_characterization._compute_uncertainty_residual_inverse_multiplicative_input_freq( - complex_response_nom_freq, complex_response_offnom_freq + residual_freq = dkpy.uncertainty_characterization._compute_residual_inverse_multiplicative_input( + complex_response_nom, complex_response_offnom ) @pytest.mark.parametrize( - "complex_response_nom_freq, complex_response_offnom_freq", + "complex_response_nom, complex_response_offnom", [ ( np.array([[1, 0], [0, 1], [1, 1]]).T, @@ -214,15 +209,15 @@ def test_compute_uncertainty_residual_inverse_multiplicative_input_freq( ), ], ) - def test_compute_uncertainty_residual_inverse_multiplicative_output_freq( - self, complex_response_nom_freq, complex_response_offnom_freq + def test_compute_residual_inverse_multiplicative_output_error( + self, complex_response_nom, complex_response_offnom ): """Test ValueError of - :func:`_compute_uncertainty_residual_inverse_multiplicative_output_freq`. + :func:`_compute_residual_inverse_multiplicative_output`. """ with pytest.raises(ValueError): - residual_freq = dkpy.uncertainty_characterization._compute_uncertainty_residual_inverse_multiplicative_output_freq( - complex_response_nom_freq, complex_response_offnom_freq + residual_freq = dkpy.uncertainty_characterization._compute_residual_inverse_multiplicative_output( + complex_response_nom, complex_response_offnom ) @@ -641,15 +636,238 @@ def test_compute_uncertainty_weight_response( # Regression testing ndarrays_regression.check( complex_response_weight_dict, - default_tolerance=dict(atol=1e-5, rtol=0), + default_tolerance=dict(atol=1e-3, rtol=0), ) -class TestFitUncertaintyWeight: - """Test :func:`fit_uncertainty_weight`.""" +class TestTransformUncertaintyModelToAdditive: + @pytest.mark.parametrize( + "complex_nominal, complex_weight_left, complex_weight_right, " + "complex_nominal_add_exp, complex_weight_left_add_exp, " + "complex_weight_right_add_exp", + [(3 * np.eye(3), np.eye(3), np.eye(3), 3 * np.eye(3), np.eye(3), np.eye(3))], + ) + def test_transform_additive_to_additive( + self, + complex_nominal, + complex_weight_left, + complex_weight_right, + complex_nominal_add_exp, + complex_weight_left_add_exp, + complex_weight_right_add_exp, + ): + from dkpy.uncertainty_characterization import _transform_additive_to_additive + + complex_additive = _transform_additive_to_additive( + complex_nominal, complex_weight_left, complex_weight_right + ) + complex_nominal_add = complex_additive[0] + complex_weight_left_add = complex_additive[1] + complex_weight_right_add = complex_additive[2] + + assert np.allclose(complex_nominal_add, complex_nominal_add_exp) + assert np.allclose(complex_weight_left_add, complex_weight_left_add_exp) + assert np.allclose(complex_weight_right_add, complex_weight_right_add_exp) @pytest.mark.parametrize( - "sys_nom, sys_offnom_list, omega, weight_left_structure, weight_right_structure, fit_order", + "complex_nominal, complex_weight_left, complex_weight_right, " + "complex_nominal_add_exp, complex_weight_left_add_exp, " + "complex_weight_right_add_exp", + [ + ( + 3 * np.eye(3), + np.eye(3), + np.eye(3), + 3 * np.eye(3), + 3 * np.eye(3), + np.eye(3), + ) + ], + ) + def test_transform_multiplicative_input_to_additive( + self, + complex_nominal, + complex_weight_left, + complex_weight_right, + complex_nominal_add_exp, + complex_weight_left_add_exp, + complex_weight_right_add_exp, + ): + from dkpy.uncertainty_characterization import ( + _transform_multiplicative_input_to_additive, + ) + + complex_additive = _transform_multiplicative_input_to_additive( + complex_nominal, complex_weight_left, complex_weight_right + ) + complex_nominal_add = complex_additive[0] + complex_weight_left_add = complex_additive[1] + complex_weight_right_add = complex_additive[2] + + assert np.allclose(complex_nominal_add, complex_nominal_add_exp) + assert np.allclose(complex_weight_left_add, complex_weight_left_add_exp) + assert np.allclose(complex_weight_right_add, complex_weight_right_add_exp) + + @pytest.mark.parametrize( + "complex_nominal, complex_weight_left, complex_weight_right, " + "complex_nominal_add_exp, complex_weight_left_add_exp, " + "complex_weight_right_add_exp", + [ + ( + 3 * np.eye(3), + np.eye(3), + np.eye(3), + 3 * np.eye(3), + np.eye(3), + 3 * np.eye(3), + ) + ], + ) + def test_transform_multiplicative_output_to_additive( + self, + complex_nominal, + complex_weight_left, + complex_weight_right, + complex_nominal_add_exp, + complex_weight_left_add_exp, + complex_weight_right_add_exp, + ): + from dkpy.uncertainty_characterization import ( + _transform_multiplicative_output_to_additive, + ) + + complex_additive = _transform_multiplicative_output_to_additive( + complex_nominal, complex_weight_left, complex_weight_right + ) + complex_nominal_add = complex_additive[0] + complex_weight_left_add = complex_additive[1] + complex_weight_right_add = complex_additive[2] + + assert np.allclose(complex_nominal_add, complex_nominal_add_exp) + assert np.allclose(complex_weight_left_add, complex_weight_left_add_exp) + assert np.allclose(complex_weight_right_add, complex_weight_right_add_exp) + + @pytest.mark.parametrize( + "complex_nominal, complex_weight_left, complex_weight_right, " + "complex_nominal_add_exp, complex_weight_left_add_exp, " + "complex_weight_right_add_exp", + [ + ( + 3 * np.eye(3), + 0.5 * np.eye(3), + 0.5 * np.eye(3), + 6.85714 * np.eye(3), + -2.26779 * np.eye(3), + 2.26779 * np.eye(3), + ) + ], + ) + def test_transform_inverse_additive_to_additive( + self, + complex_nominal, + complex_weight_left, + complex_weight_right, + complex_nominal_add_exp, + complex_weight_left_add_exp, + complex_weight_right_add_exp, + ): + from dkpy.uncertainty_characterization import ( + _transform_inverse_additive_to_additive, + ) + + complex_additive = _transform_inverse_additive_to_additive( + complex_nominal, complex_weight_left, complex_weight_right + ) + complex_nominal_add = complex_additive[0] + complex_weight_left_add = complex_additive[1] + complex_weight_right_add = complex_additive[2] + + assert np.allclose(complex_nominal_add, complex_nominal_add_exp) + assert np.allclose(complex_weight_left_add, complex_weight_left_add_exp) + assert np.allclose(complex_weight_right_add, complex_weight_right_add_exp) + + @pytest.mark.parametrize( + "complex_nominal, complex_weight_left, complex_weight_right, " + "complex_nominal_add_exp, complex_weight_left_add_exp, " + "complex_weight_right_add_exp", + [ + ( + 3 * np.eye(3), + 0.5 * np.eye(3), + 0.5 * np.eye(3), + 3.2 * np.eye(3), + -1.54919 * np.eye(3), + 0.51640 * np.eye(3), + ) + ], + ) + def test_transform_inverse_multiplicative_input_to_additive( + self, + complex_nominal, + complex_weight_left, + complex_weight_right, + complex_nominal_add_exp, + complex_weight_left_add_exp, + complex_weight_right_add_exp, + ): + from dkpy.uncertainty_characterization import ( + _transform_inverse_multiplicative_input_to_additive, + ) + + complex_additive = _transform_inverse_multiplicative_input_to_additive( + complex_nominal, complex_weight_left, complex_weight_right + ) + complex_nominal_add = complex_additive[0] + complex_weight_left_add = complex_additive[1] + complex_weight_right_add = complex_additive[2] + + assert np.allclose(complex_nominal_add, complex_nominal_add_exp) + assert np.allclose(complex_weight_left_add, complex_weight_left_add_exp) + assert np.allclose(complex_weight_right_add, complex_weight_right_add_exp) + + @pytest.mark.parametrize( + "complex_nominal, complex_weight_left, complex_weight_right, " + "complex_nominal_add_exp, complex_weight_left_add_exp, " + "complex_weight_right_add_exp", + [ + ( + 3 * np.eye(3), + 0.5 * np.eye(3), + 0.5 * np.eye(3), + 3.2 * np.eye(3), + -0.51640 * np.eye(3), + 1.54919 * np.eye(3), + ) + ], + ) + def test_transform_inverse_multiplicative_output_to_additive( + self, + complex_nominal, + complex_weight_left, + complex_weight_right, + complex_nominal_add_exp, + complex_weight_left_add_exp, + complex_weight_right_add_exp, + ): + from dkpy.uncertainty_characterization import ( + _transform_inverse_multiplicative_output_to_additive, + ) + + complex_additive = _transform_inverse_multiplicative_output_to_additive( + complex_nominal, complex_weight_left, complex_weight_right + ) + complex_nominal_add = complex_additive[0] + complex_weight_left_add = complex_additive[1] + complex_weight_right_add = complex_additive[2] + + assert np.allclose(complex_nominal_add, complex_nominal_add_exp) + assert np.allclose(complex_weight_left_add, complex_weight_left_add_exp) + assert np.allclose(complex_weight_right_add, complex_weight_right_add_exp) + + +class TestComputeUncertaintyMeasureResponse: + @pytest.mark.parametrize( + "nominal, offnominal, omega, uncertainty_model", [ ( control.TransferFunction([1], [0.5, 1]), @@ -660,22 +878,200 @@ class TestFitUncertaintyWeight: control.TransferFunction([1], [0.7, 1]), ], np.logspace(-2, 2, 100), - "diagonal", - "diagonal", - 4, + "additive", ), ( - control.TransferFunction([1], [1, 2 * 0.5 * 1, 1**2]), + control.TransferFunction([1], [0.5, 1]), [ - 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]), + 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(-1, 1, 100), - "diagonal", - "diagonal", - 4, + np.logspace(-2, 2, 100), + "multiplicative_input", + ), + ( + 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), + "multiplicative_output", + ), + ( + 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), + "inverse_additive", + ), + ( + 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), + "inverse_multiplicative_input", + ), + ( + 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), + "inverse_multiplicative_output", + ), + ( + 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), + "additive", + ), + ( + 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), + "multiplicative_input", + ), + ( + 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), + "multiplicative_output", ), ( control.TransferFunction( @@ -721,6 +1117,181 @@ class TestFitUncertaintyWeight: ), ], np.logspace(-1, 1.5, 100), + "inverse_additive", + ), + ( + 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), + "inverse_multiplicative_input", + ), + ( + 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), + "inverse_multiplicative_output", + ), + ], + ) + def test_compute_uncertainty_measure_response( + self, + ndarrays_regression, + nominal, + offnominal, + omega, + uncertainty_model, + ): + """Regression test for :func:`compute_uncertainty_measure_response`.""" + + # Frequency response of systems + frd_nominal = control.frequency_response(nominal, omega, squeeze=False) + frd_offnominal_list = control.frequency_response( + offnominal, omega, squeeze=False + ) + + # Uncertainty residual computation + uncertainty_model_list = [uncertainty_model] + complex_residual_dict = dkpy.compute_uncertainty_residual_response( + frd_nominal, + frd_offnominal_list, + uncertainty_model_list, + ) + + measure_dict = {} + for weight_structure in ["scalar", "diagonal", "full"]: + # Uncertainty weight computation + complex_weight_left, complex_weight_right = ( + dkpy.compute_uncertainty_weight_response( + complex_residual_dict[uncertainty_model], + weight_structure, + weight_structure, + ) + ) + + # Uncertainty measure computation + measure = dkpy.compute_uncertainty_measure_response( + frd_nominal, + complex_weight_left, + complex_weight_right, + uncertainty_model, + ) + + measure_dict[weight_structure] = measure + + ndarrays_regression.check( + measure_dict, + default_tolerance=dict(atol=1e-2, rtol=1e-2), + ) + + +class TestFitUncertaintyWeight: + """Test :func:`fit_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, @@ -771,7 +1342,7 @@ class TestFitUncertaintyWeight: np.logspace(-1, 1.5, 100), "diagonal", "diagonal", - [4, 4], + 4, ), ( control.TransferFunction( @@ -819,7 +1390,7 @@ class TestFitUncertaintyWeight: np.logspace(-1, 1.5, 100), "diagonal", "diagonal", - np.array([4.0, 4.0]), + [4, 4], ), ], ) @@ -871,10 +1442,18 @@ def test_fit_uncertainty_weight( # Overbounding transfer function fit weight_left_fit = dkpy.fit_uncertainty_weight( - complex_response_weight_left, omega, fit_order + complex_response_weight_left, + omega, + fit_order, + "left", + weight_left_structure, ) weight_right_fit = dkpy.fit_uncertainty_weight( - complex_response_weight_right, omega, fit_order + complex_response_weight_right, + omega, + fit_order, + "right", + weight_right_structure, ) # Uncertainty weight fit frequency response @@ -898,7 +1477,7 @@ def test_fit_uncertainty_weight( # Regression testing ndarrays_regression.check( complex_response_weight_fit_dict, - default_tolerance=dict(atol=1e-5, rtol=0), + default_tolerance=dict(atol=1e-2, rtol=1e-2), ) diff --git a/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal0_offnominal0_omega0_additive_.npz b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal0_offnominal0_omega0_additive_.npz new file mode 100644 index 0000000..3f1b3bd Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal0_offnominal0_omega0_additive_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal10_offnominal10_omega10_inverse_multiplicative_input_.npz b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal10_offnominal10_omega10_inverse_multiplicative_input_.npz new file mode 100644 index 0000000..6a31355 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal10_offnominal10_omega10_inverse_multiplicative_input_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal11_offnominal11_omega11_inverse_multiplicative_output_.npz b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal11_offnominal11_omega11_inverse_multiplicative_output_.npz new file mode 100644 index 0000000..546b286 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal11_offnominal11_omega11_inverse_multiplicative_output_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal1_offnominal1_omega1_multiplicative_input_.npz b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal1_offnominal1_omega1_multiplicative_input_.npz new file mode 100644 index 0000000..c47b6dc Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal1_offnominal1_omega1_multiplicative_input_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal2_offnominal2_omega2_multiplicative_output_.npz b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal2_offnominal2_omega2_multiplicative_output_.npz new file mode 100644 index 0000000..91e59f3 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal2_offnominal2_omega2_multiplicative_output_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal3_offnominal3_omega3_inverse_additive_.npz b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal3_offnominal3_omega3_inverse_additive_.npz new file mode 100644 index 0000000..80fdef4 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal3_offnominal3_omega3_inverse_additive_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal4_offnominal4_omega4_inverse_multiplicative_input_.npz b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal4_offnominal4_omega4_inverse_multiplicative_input_.npz new file mode 100644 index 0000000..77bf3bc Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal4_offnominal4_omega4_inverse_multiplicative_input_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal5_offnominal5_omega5_inverse_multiplicative_output_.npz b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal5_offnominal5_omega5_inverse_multiplicative_output_.npz new file mode 100644 index 0000000..4882b12 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal5_offnominal5_omega5_inverse_multiplicative_output_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal6_offnominal6_omega6_additive_.npz b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal6_offnominal6_omega6_additive_.npz new file mode 100644 index 0000000..a92e12a Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal6_offnominal6_omega6_additive_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal7_offnominal7_omega7_multiplicative_input_.npz b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal7_offnominal7_omega7_multiplicative_input_.npz new file mode 100644 index 0000000..e294181 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal7_offnominal7_omega7_multiplicative_input_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal8_offnominal8_omega8_multiplicative_output_.npz b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal8_offnominal8_omega8_multiplicative_output_.npz new file mode 100644 index 0000000..7f2a674 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal8_offnominal8_omega8_multiplicative_output_.npz differ diff --git a/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal9_offnominal9_omega9_inverse_additive_.npz b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal9_offnominal9_omega9_inverse_additive_.npz new file mode 100644 index 0000000..a1ce2a4 Binary files /dev/null and b/tests/test_uncertainty_characterization/test_compute_uncertainty_measure_response_nominal9_offnominal9_omega9_inverse_additive_.npz differ diff --git a/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom0_sys_offnom_list0_omega0_diagonal_diagonal_4_.npz b/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom0_sys_offnom_list0_omega0_diagonal_diagonal_4_.npz index 8164121..5f9e61f 100644 Binary files a/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom0_sys_offnom_list0_omega0_diagonal_diagonal_4_.npz and b/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom0_sys_offnom_list0_omega0_diagonal_diagonal_4_.npz differ diff --git a/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom1_sys_offnom_list1_omega1_diagonal_diagonal_4_.npz b/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom1_sys_offnom_list1_omega1_diagonal_diagonal_4_.npz index 5b40d52..820fff5 100644 Binary files a/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom1_sys_offnom_list1_omega1_diagonal_diagonal_4_.npz and b/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom1_sys_offnom_list1_omega1_diagonal_diagonal_4_.npz differ diff --git a/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom2_sys_offnom_list2_omega2_diagonal_diagonal_4_.npz b/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom2_sys_offnom_list2_omega2_diagonal_diagonal_4_.npz index b30ec1e..beb51f0 100644 Binary files a/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom2_sys_offnom_list2_omega2_diagonal_diagonal_4_.npz and b/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom2_sys_offnom_list2_omega2_diagonal_diagonal_4_.npz differ diff --git a/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom3_sys_offnom_list3_omega3_diagonal_diagonal_fit_order3_.npz b/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom3_sys_offnom_list3_omega3_diagonal_diagonal_fit_order3_.npz index b30ec1e..beb51f0 100644 Binary files a/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom3_sys_offnom_list3_omega3_diagonal_diagonal_fit_order3_.npz and b/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom3_sys_offnom_list3_omega3_diagonal_diagonal_fit_order3_.npz differ diff --git a/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom4_sys_offnom_list4_omega4_diagonal_diagonal_fit_order4_.npz b/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom4_sys_offnom_list4_omega4_diagonal_diagonal_fit_order4_.npz deleted file mode 100644 index b30ec1e..0000000 Binary files a/tests/test_uncertainty_characterization/test_fit_uncertainty_weight_sys_nom4_sys_offnom_list4_omega4_diagonal_diagonal_fit_order4_.npz and /dev/null differ