diff --git a/CITATION.cff b/CITATION.cff index 734d5cd..c7ad650 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -5,10 +5,14 @@ authors: given-names: "Steven" orcid: "https://orcid.org/0000-0003-4930-9634" affiliation: "McGill University" + - family-names: "Adams" + given-names: "Timothy Everett" + orcid: "https://orcid.org/0009-0008-5291-0686" + affiliation: "McGill University" - family-names: "Forbes" given-names: "James Richard" orcid: "https://orcid.org/0000-0002-1987-9268" affiliation: "McGill University" title: "decargroup/dkpy" -version: v0.1.9 +version: v0.1.10 url: "https://github.com/decargroup/dkpy" diff --git a/README.rst b/README.rst index 60098a3..bb19896 100644 --- a/README.rst +++ b/README.rst @@ -102,7 +102,11 @@ Example # Synthesize a controller omega = np.logspace(-3, 3, 61) - block_structure = np.array([[1, 1], [1, 1], [2, 2]]) + block_structure = [ + dkpy.ComplexFullBlock(1, 1), + dkpy.ComplexFullBlock(1, 1), + dkpy.ComplexFullBlock(2, 2), + ] K, N, mu, d_scale_fit_info, info = dk_iter.synthesize( eg["P"], eg["n_y"], @@ -131,12 +135,12 @@ If you use this software in your research, please cite it as below or see .. code-block:: bibtex - @software{dahdah_dkpy_2024, + @software{dahdah_dkpy_2025, title={{decargroup/dkpy}}, doi={10.5281/zenodo.14511244}, url={https://github.com/decargroup/dkpy}, publisher={Zenodo}, - author={Steven Dahdah and James Richard Forbes}, - version = {{v0.1.9}}, - year={2024}, + author={Steven Dahdah and Timothy Everett Adams and James Richard Forbes}, + version = {{v0.1.10}}, + year={2025}, } diff --git a/doc/conf.py b/doc/conf.py index 688441e..18484e5 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -7,10 +7,10 @@ # https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information project = "dkpy" -copyright = "2024, Steven Dahdah and James Richard Forbes" -author = "Steven Dahdah and James Richard Forbes" -version = "0.1.9" -release = "0.1.9" +copyright = "2024, Steven Dahdah, Timothy Everett Adams and James Richard Forbes" +author = "Steven Dahdah, Timothy Everett Adams and James Richard Forbes" +version = "0.1.10" +release = "0.1.10" # -- General configuration --------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration diff --git a/doc/dkpy.rst b/doc/dkpy.rst index 5d9170f..9371a18 100644 --- a/doc/dkpy.rst +++ b/doc/dkpy.rst @@ -1,3 +1,38 @@ +``dkpy`` architecture overview +============================== + +``dkpy`` implements each step in the D-K iteration algorithm as an abstract +class. This ensures that each step satisfies the interface required by the +algorithm while leaving the numerical implementation up to the user. Sensible +implementations of each of these steps are provided in ``dkpy``. However, the +user can create custom implementations via the abstract classes. Therefore, +anyone aiming to extend or customize ``dkpy`` should familiarize themselves +with them. The abstract classes are provided below. + +.. autosummary:: + :toctree: _autosummary/ + + dkpy.ControllerSynthesis + dkpy.StructuredSingularValue + dkpy.DScaleFit + dkpy.DkIteration + +The steps of the D-K iteration algorithm are as follows [SP06]_. +1. Controller synthesis (:class:`ControllerSynthesis`): Synthesize an + H-infinity controller for the scaled problem with fixed fitted D-scales. +2. Structured singular value and D-scale computation + (:class:`StructuredSingularValue`): Compute the structured singular value + and the D-scales over a discrete grid of frequencies with a fixed + controller. +3. D-scale fit (:class:`DScaleFit`): Fit the magnitude of each D-scale to a + stable minimum-phase LTI system. + +The D-K iteration algorithm, specified by a :class:`DkIteration` object, loops +through these three steps until the robustness criteria are satisfied. +Therefore, the algorithm is specified using implementations of each of the +abstract classes :class:`ControllerSynthesis`, +:class:`StructuredSingularValue`, and :class:`DScaleFit`. + D-K iteration methods ===================== @@ -69,17 +104,23 @@ selecting the order in :func:`DScaleFit.fit`. dkpy.DScaleFitSlicot -Extending ``dkpy`` -================== +Uncertainty block structure +=========================== -The abstract classes defining the structure of ``dkpy`` are presented below. -Anyone aiming to extend or customize ``dkpy`` should familiarize themselves -with them. +The uncertainty block structure is specified via an +:class:`UncertaintyBlockStructure` object, which encodes the block diagonal +uncertainty structure. The :class:`UncertaintyBlockStructure` object is +composed of individual uncertainty blocks that satisfy the interface in +:class:`UncertaintyBlock`, which are provided below. .. autosummary:: :toctree: _autosummary/ - dkpy.DkIteration - dkpy.ControllerSynthesis - dkpy.StructuredSingularValue - dkpy.DScaleFit + dkpy.RealDiagonalBlock + dkpy.ComplexDiagonalBlock + dkpy.ComplexFullBlock + +The :class:`UncertaintyBlockStructure` object can also be specified using the +MATLAB two-column array format (see `MATLAB documentation +`) for users that are +more comfortable with this notation. diff --git a/examples/1_example_dk_iter_fixed_order.py b/examples/1_example_dk_iter_fixed_order.py index ba9340a..318dfa4 100644 --- a/examples/1_example_dk_iter_fixed_order.py +++ b/examples/1_example_dk_iter_fixed_order.py @@ -19,13 +19,23 @@ def example_dk_iter_fixed_order(): ) omega = np.logspace(-3, 3, 61) - block_structure = np.array([[1, 1], [1, 1], [2, 2]]) + # Alternative MATLAB descr. + # uncertainty_structure = dkpy.UncertaintyBlockStructure( + # [[1, 1], [1, 1], [2, 2]] + # ) + uncertainty_structure = dkpy.UncertaintyBlockStructure( + [ + dkpy.ComplexFullBlock(1, 1), + dkpy.ComplexFullBlock(1, 1), + dkpy.ComplexFullBlock(2, 2), + ] + ) K, N, mu, iter_results, info = dk_iter.synthesize( eg["P"], eg["n_y"], eg["n_u"], omega, - block_structure, + uncertainty_structure, ) print(f"mu={mu}") diff --git a/examples/2_example_dk_iter_list_order.py b/examples/2_example_dk_iter_list_order.py index 2bcbe2d..96d8484 100644 --- a/examples/2_example_dk_iter_list_order.py +++ b/examples/2_example_dk_iter_list_order.py @@ -23,13 +23,23 @@ def example_dk_iter_list_order(): ) omega = np.logspace(-3, 3, 61) - block_structure = np.array([[1, 1], [1, 1], [2, 2]]) + # Alternative MATLAB descr. + # uncertainty_structure = dkpy.UncertaintyBlockStructure( + # [[1, 1], [1, 1], [2, 2]] + # ) + uncertainty_structure = dkpy.UncertaintyBlockStructure( + [ + dkpy.ComplexFullBlock(1, 1), + dkpy.ComplexFullBlock(1, 1), + dkpy.ComplexFullBlock(2, 2), + ] + ) K, N, mu, iter_results, info = dk_iter.synthesize( eg["P"], eg["n_y"], eg["n_u"], omega, - block_structure, + uncertainty_structure, ) print(f"mu={mu}") diff --git a/examples/3_example_dk_iter_auto_order.py b/examples/3_example_dk_iter_auto_order.py index 1770546..9ba2b17 100644 --- a/examples/3_example_dk_iter_auto_order.py +++ b/examples/3_example_dk_iter_auto_order.py @@ -44,13 +44,23 @@ def example_dk_iter_auto_order(): ) omega = np.logspace(-3, 3, 61) - block_structure = np.array([[1, 1], [1, 1], [2, 2]]) + # Alternative MATLAB descr. + # uncertainty_structure = dkpy.UncertaintyBlockStructure( + # [[1, 1], [1, 1], [2, 2]] + # ) + uncertainty_structure = dkpy.UncertaintyBlockStructure( + [ + dkpy.ComplexFullBlock(1, 1), + dkpy.ComplexFullBlock(1, 1), + dkpy.ComplexFullBlock(2, 2), + ] + ) K, N, mu, iter_results, info = dk_iter.synthesize( eg["P"], eg["n_y"], eg["n_u"], omega, - block_structure, + uncertainty_structure, ) print(f"mu={mu}") diff --git a/examples/4_example_dk_iter_interactive.py b/examples/4_example_dk_iter_interactive.py index d713ff3..368e9c3 100644 --- a/examples/4_example_dk_iter_interactive.py +++ b/examples/4_example_dk_iter_interactive.py @@ -37,13 +37,23 @@ def example_dk_iter_interactive(): ) omega = np.logspace(-3, 3, 61) - block_structure = np.array([[1, 1], [1, 1], [2, 2]]) + # Alternative MATLAB descr. + # uncertainty_structure = dkpy.UncertaintyBlockStructure( + # [[1, 1], [1, 1], [2, 2]] + # ) + uncertainty_structure = dkpy.UncertaintyBlockStructure( + [ + dkpy.ComplexFullBlock(1, 1), + dkpy.ComplexFullBlock(1, 1), + dkpy.ComplexFullBlock(2, 2), + ] + ) K, N, mu, iter_results, info = dk_iter.synthesize( eg["P"], eg["n_y"], eg["n_u"], omega, - block_structure, + uncertainty_structure, ) print(f"mu={mu}") diff --git a/examples/5_example_dk_iteration_custom.py b/examples/5_example_dk_iteration_custom.py index 0c33e6b..8f01287 100644 --- a/examples/5_example_dk_iteration_custom.py +++ b/examples/5_example_dk_iteration_custom.py @@ -34,7 +34,7 @@ def _get_fit_order( D_omega, P, K, - block_structure, + uncertainty_structure, ): print(f"Iteration {self.my_iter_count} with mu of {np.max(mu_omega)}") if self.my_iter_count < 3: @@ -70,13 +70,23 @@ def example_dk_iter_custom(): ) omega = np.logspace(-3, 3, 61) - block_structure = np.array([[1, 1], [1, 1], [2, 2]]) + # Alternative MATLAB descr. + # uncertainty_structure = dkpy.UncertaintyBlockStructure( + # [[1, 1], [1, 1], [2, 2]] + # ) + uncertainty_structure = dkpy.UncertaintyBlockStructure( + [ + dkpy.ComplexFullBlock(1, 1), + dkpy.ComplexFullBlock(1, 1), + dkpy.ComplexFullBlock(2, 2), + ] + ) K, N, mu, iter_results, info = dk_iter.synthesize( eg["P"], eg["n_y"], eg["n_u"], omega, - block_structure, + uncertainty_structure, ) print(f"mu={mu}") diff --git a/pyproject.toml b/pyproject.toml index ce8267b..93923f9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,11 +4,11 @@ build-backend = "hatchling.build" [project] name = "dkpy" -version = "0.1.9" +version = "0.1.10" dependencies = [ "numpy>=1.21.0", "scipy>=1.7.0", - "control>=0.10.0", + "control>=0.10.2", "slycot>=0.6.0", "cvxpy>=1.5.0", "joblib>=1.4.0", diff --git a/requirements.txt b/requirements.txt index 932fbf3..a7b80b6 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,7 +1,7 @@ # Dependencies numpy scipy -control +control>=0.10.2 slycot cvxpy joblib diff --git a/src/dkpy/__init__.py b/src/dkpy/__init__.py index a34174b..de6b21e 100644 --- a/src/dkpy/__init__.py +++ b/src/dkpy/__init__.py @@ -4,4 +4,5 @@ from .dk_iteration import * from .d_scale_fit import * from .structured_singular_value import * +from .uncertainty_structure import * from .utilities import * diff --git a/src/dkpy/controller_synthesis.py b/src/dkpy/controller_synthesis.py index 0b4647c..4e25551 100644 --- a/src/dkpy/controller_synthesis.py +++ b/src/dkpy/controller_synthesis.py @@ -71,7 +71,7 @@ class HinfSynSlicot(ControllerSynthesis): >>> P, n_y, n_u = example_scherer1997_p907 >>> K, N, gamma, info = dkpy.HinfSynSlicot().synthesize(P, n_y, n_u) >>> gamma - 9.5080 + 9.51 """ def __init__(self): @@ -107,7 +107,7 @@ class HinfSynLmi(ControllerSynthesis): >>> P, n_y, n_u = example_scherer1997_p907 >>> K, N, gamma, info = dkpy.HinfSynLmi().synthesize(P, n_y, n_u) >>> gamma - 9.5081 + 9.51 H-infinity controller synthesis with CLARABEL @@ -407,7 +407,7 @@ class HinfSynLmiBisection(ControllerSynthesis): ... }, ... ).synthesize(P, n_y, n_u) >>> gamma - 9.5093 + 9.51 """ def __init__( diff --git a/src/dkpy/d_scale_fit.py b/src/dkpy/d_scale_fit.py index c3281c7..cecb99d 100644 --- a/src/dkpy/d_scale_fit.py +++ b/src/dkpy/d_scale_fit.py @@ -6,7 +6,7 @@ ] import abc -from typing import Optional, Tuple, Union +from typing import Tuple, Union, List, Optional import warnings import control @@ -14,7 +14,7 @@ import scipy.linalg import slycot -from . import utilities +from . import uncertainty_structure class DScaleFit(metaclass=abc.ABCMeta): @@ -26,7 +26,13 @@ def fit( omega: np.ndarray, D_omega: np.ndarray, order: Union[int, np.ndarray] = 0, - block_structure: Optional[np.ndarray] = None, + block_structure: Optional[ + Union[ + List[uncertainty_structure.UncertaintyBlock], + List[List[int]], + np.ndarray, + ] + ] = None, ) -> Tuple[control.StateSpace, control.StateSpace]: """Fit D-scale magnitudes. @@ -39,10 +45,8 @@ def fit( dimension. order : Union[int, np.ndarray] Transfer function order to fit. Can be specified per-entry. - block_structure : np.ndarray - 2D array with 2 columns and as many rows as uncertainty blocks - in Delta. The columns represent the number of rows and columns in - each uncertainty block. See [#mussv]_. + block_structure : Optional[Union[List[uncertainty_structure.UncertaintyBlock], List[List[int], np.ndarray]] + Uncertainty block structure description. Returns ------- @@ -53,7 +57,7 @@ def fit( ------ ValueError If ``order`` is an array but its dimensions are inconsistent with - ``block_structure``. + ``uncertainty_structure``. References ---------- @@ -70,7 +74,11 @@ class DScaleFitSlicot(DScaleFit): Compute ``mu`` and ``D`` at each frequency and fit a transfer matrix to ``D`` >>> P, n_y, n_u, K = example_skogestad2006_p325 - >>> block_structure = np.array([[1, 1], [1, 1], [2, 2]]) + >>> block_structure = [ + ... dkpy.ComplexFullBlock(1, 1), + ... dkpy.ComplexFullBlock(1, 1), + ... dkpy.ComplexFullBlock(2, 2), + ... ] >>> omega = np.logspace(-3, 3, 61) >>> N = P.lft(K) >>> N_omega = N(1j * omega) @@ -79,26 +87,6 @@ class DScaleFitSlicot(DScaleFit): ... block_structure, ... ) >>> D, D_inv = dkpy.DScaleFitSlicot().fit(omega, D_omega, 2, block_structure) - >>> print(control.ss2tf(D[0, 0])) - : sys[3641]$indexed$converted - Inputs (1): ['u[0]'] - Outputs (1): ['y[0]'] - - - 0.0157 s^2 + 0.257 s + 0.1391 - ----------------------------- - s^2 + 0.9658 s + 0.01424 - - >>> print(control.ss2tf(D[1, 1])) - : sys[3641]$indexed$converted - Inputs (1): ['u[1]'] - Outputs (1): ['y[1]'] - - - 0.01573 s^2 + 0.2574 s + 0.1394 - ------------------------------- - s^2 + 0.9663 s + 0.01425 - """ def fit( @@ -106,19 +94,30 @@ def fit( omega: np.ndarray, D_omega: np.ndarray, order: Union[int, np.ndarray] = 0, - block_structure: Optional[np.ndarray] = None, + block_structure: Optional[ + Union[ + List[uncertainty_structure.UncertaintyBlock], + List[List[int]], + np.ndarray, + ] + ] = None, ) -> Tuple[control.StateSpace, control.StateSpace]: # Get mask if block_structure is None: mask = -1 * np.ones((D_omega.shape[0], D_omega.shape[1]), dtype=int) else: - mask = _mask_from_block_structure(block_structure) + block_structure = ( + uncertainty_structure._convert_block_structure_representation( + block_structure + ) + ) + mask = _generate_d_scale_mask(block_structure) # Check order dimensions orders = order if isinstance(order, np.ndarray) else order * np.ones_like(mask) if orders.shape != mask.shape: raise ValueError( "`order` must be an integer or an array whose dimensions are " - "consistent with `block_structure`." + "consistent with `uncertainty_structure`." ) # Transfer matrix tf_array = np.zeros((D_omega.shape[0], D_omega.shape[1]), dtype=object) @@ -149,48 +148,47 @@ def fit( ) sys = control.StateSpace(A, B, C, D, dt=0) tf_array[row, col] = control.ss2tf(sys) - tf = utilities._tf_combine(tf_array) + tf = control.combine_tf(tf_array) ss = control.tf2ss(tf) ss_inv = _invert_biproper_ss(ss) return ss, ss_inv -def _mask_from_block_structure(block_structure: np.ndarray) -> np.ndarray: - """Create a mask from a specified block structure. +def _generate_d_scale_mask( + block_structure: List[uncertainty_structure.UncertaintyBlock], +) -> np.ndarray: + """Create a mask for the D-scale fit for the block structure. Entries known to be zero are set to 0. Entries known to be one are set to 1. Entries to be fit numerically are set to -1. Parameters ---------- - block_structure : np.ndarray - 2D array with 2 columns and as many rows as uncertainty blocks - in Delta. The columns represent the number of rows and columns in - each uncertainty block. See [#mussv]_. + block_structure : List[uncertainty_structure.UncertaintyBlock] + Uncertainty block structure description. Returns ------- np.ndarray Array of integers indicating zero, one, and unknown elements in the block structure. - - References - ---------- - .. [#mussv] https://www.mathworks.com/help/robust/ref/mussv.html """ + num_blocks = len(block_structure) X_lst = [] - for i in range(block_structure.shape[0]): - if block_structure[i, 0] <= 0: + for i in range(num_blocks): + # Uncertainty block + block = block_structure[i] + if not block.is_complex: raise NotImplementedError("Real perturbations are not yet supported.") - if block_structure[i, 1] <= 0: + if block.is_diagonal: raise NotImplementedError("Diagonal perturbations are not yet supported.") - if block_structure[i, 0] != block_structure[i, 1]: + if not block.is_square: raise NotImplementedError("Nonsquare perturbations are not yet supported.") # Set last scaling to identity - if i == block_structure.shape[0] - 1: - X_lst.append(np.eye(block_structure[i, 0], dtype=int)) + if i == num_blocks - 1: + X_lst.append(np.eye(block.num_inputs, dtype=int)) else: - X_lst.append(-1 * np.eye(block_structure[i, 0], dtype=int)) + X_lst.append(-1 * np.eye(block.num_inputs, dtype=int)) X = scipy.linalg.block_diag(*X_lst) return X diff --git a/src/dkpy/dk_iteration.py b/src/dkpy/dk_iteration.py index d875a0e..a522ea3 100644 --- a/src/dkpy/dk_iteration.py +++ b/src/dkpy/dk_iteration.py @@ -13,7 +13,7 @@ import abc import logging -from typing import Any, Dict, List, Optional, Tuple, Union +from typing import Any, Dict, List, Optional, Tuple, Union, List import control import numpy as np @@ -24,6 +24,7 @@ controller_synthesis, d_scale_fit, structured_singular_value, + uncertainty_structure, utilities, ) @@ -46,7 +47,9 @@ def __init__( mu_fit_omega: np.ndarray, D_fit_omega: np.ndarray, D_fit: control.StateSpace, - block_structure: np.ndarray, + block_structure: Union[ + List[uncertainty_structure.UncertaintyBlock], List[List[int]], np.ndarray + ], ): """Instantiate :class:`IterResult`. @@ -64,14 +67,8 @@ def __init__( Fit D-scale magnitude at each frequency. D_fit : control.StateSpace Fit D-scale state-space representation. - block_structure : np.ndarray - 2D array with 2 columns and as many rows as uncertainty blocks - in Delta. The columns represent the number of rows and columns in - each uncertainty block. See [#mussv]_. - - References - ---------- - .. [#mussv] https://www.mathworks.com/help/robust/ref/mussv.html + uncertainty_structure : Union[List[uncertainty_structure.UncertaintyBlock], List[List[int]], np.ndarray] + Uncertainty block structure representation. """ self.omega = omega self.mu_omega = mu_omega @@ -79,7 +76,11 @@ def __init__( self.mu_fit_omega = mu_fit_omega self.D_fit_omega = D_fit_omega self.D_fit = D_fit - self.block_structure = block_structure + self.block_structure = ( + uncertainty_structure._convert_block_structure_representation( + block_structure + ) + ) @classmethod def create_from_fit( @@ -91,7 +92,9 @@ def create_from_fit( K: control.StateSpace, D_fit: control.StateSpace, D_fit_inv: control.StateSpace, - block_structure: np.ndarray, + block_structure: Union[ + List[uncertainty_structure.UncertaintyBlock], List[List[int]], np.ndarray + ], ) -> "IterResult": """Instantiate :class:`IterResult` from fit D-scales. @@ -111,10 +114,8 @@ def create_from_fit( Fit D-scale magnitude at each frequency. D_fit_inv : control.StateSpace Fit inverse D-scale magnitude at each frequency. - block_structure : np.ndarray - 2D array with 2 columns and as many rows as uncertainty blocks - in Delta. The columns represent the number of rows and columns in - each uncertainty block. See [#mussv]_. + block_structure : Union[List[uncertainty_structure.UncertaintyBlock], List[List[int]], np.ndarray] + Uncertainty block structure representation. Returns ------- @@ -126,7 +127,11 @@ def create_from_fit( Create a ``IterResult`` object from fit data >>> P, n_y, n_u, K = example_skogestad2006_p325 - >>> block_structure = np.array([[1, 1], [1, 1], [2, 2]]) + >>> block_structure = [ + ... dkpy.ComplexFullBlock(1, 1), + ... dkpy.ComplexFullBlock(1, 1), + ... dkpy.ComplexFullBlock(2, 2), + ... ] >>> omega = np.logspace(-3, 3, 61) >>> N = P.lft(K) >>> N_omega = N(1j * omega) @@ -145,11 +150,8 @@ def create_from_fit( ... D_inv, ... block_structure, ... ) - - References - ---------- - .. [#mussv] https://www.mathworks.com/help/robust/ref/mussv.html """ + # Compute ``mu(omega)`` based on fit D-scales N = P.lft(K) scaled_cl = (D_fit * N * D_fit_inv)(1j * omega) @@ -204,7 +206,9 @@ def synthesize( n_y: int, n_u: int, omega: np.ndarray, - block_structure: np.ndarray, + block_structure: Union[ + List[uncertainty_structure.UncertaintyBlock], List[List[int]], np.ndarray + ], ) -> Tuple[ control.StateSpace, control.StateSpace, @@ -228,10 +232,8 @@ def synthesize( Number of controller outputs. omega : np.ndarray Angular frequencies to evaluate D-scales (rad/s). - block_structure : np.ndarray - 2D array with 2 columns and as many rows as uncertainty blocks - in Delta. The columns represent the number of rows and columns in - each uncertainty block. See [#mussv]_. + block_structure : Optional[List[uncertainty_structure.UncertaintyBlock], List[List[int]], np.ndarray] + Uncertainty block structure representation. Returns ------- @@ -251,11 +253,11 @@ def synthesize( object. :func:`plot_D` Plot D-scale fit from an :class:`IterResult` object. - - References - ---------- - .. [#mussv] https://www.mathworks.com/help/robust/ref/mussv.html """ + # Convert uncertainty block structure representation + block_structure = uncertainty_structure._convert_block_structure_representation( + block_structure + ) # Solution information info = {} d_scale_fit_info = [] @@ -341,7 +343,7 @@ def _get_fit_order( D_omega: np.ndarray, P: control.StateSpace, K: control.StateSpace, - block_structure: np.ndarray, + block_structure: List[uncertainty_structure.UncertaintyBlock], ) -> Optional[Union[int, np.ndarray]]: """Get D-scale fit order. @@ -359,19 +361,13 @@ def _get_fit_order( Generalized plant. K : control.StateSpace Controller. - block_structure : np.ndarray - 2D array with 2 columns and as many rows as uncertainty blocks - in Delta. The columns represent the number of rows and columns in - each uncertainty block. See [#mussv]_. + block_structure : List[uncertainty_structure.UncertaintyBlock] + Uncertainty block structure representation. Returns ------- Optional[Union[int, np.ndarray]] D-scale fit order. If ``None``, iteration ends. - - References - ---------- - .. [#mussv] https://www.mathworks.com/help/robust/ref/mussv.html """ raise NotImplementedError() @@ -413,7 +409,11 @@ def __init__( ... fit_order=4, ... ) >>> omega = np.logspace(-3, 3, 61) - >>> block_structure = np.array([[1, 1], [1, 1], [2, 2]]) + >>> block_structure = [ + ... dkpy.ComplexFullBlock(1, 1), + ... dkpy.ComplexFullBlock(1, 1), + ... dkpy.ComplexFullBlock(2, 2), + ... ] >>> K, N, mu, d_scale_fit_info, info = dk_iter.synthesize( ... eg["P"], ... eg["n_y"], @@ -422,7 +422,7 @@ def __init__( ... block_structure, ... ) >>> mu - 1.0360 + 1.03 """ super().__init__( controller_synthesis, @@ -440,7 +440,7 @@ def _get_fit_order( D_omega: np.ndarray, P: control.StateSpace, K: control.StateSpace, - block_structure: np.ndarray, + block_structure: List[uncertainty_structure.UncertaintyBlock], ) -> Optional[Union[int, np.ndarray]]: if iteration < self.n_iterations: return self.fit_order @@ -481,7 +481,11 @@ def __init__( ... fit_orders=[4, 4, 4], ... ) >>> omega = np.logspace(-3, 3, 61) - >>> block_structure = np.array([[1, 1], [1, 1], [2, 2]]) + >>> block_structure = [ + ... dkpy.ComplexFullBlock(1, 1), + ... dkpy.ComplexFullBlock(1, 1), + ... dkpy.ComplexFullBlock(2, 2), + ... ] >>> K, N, mu, d_scale_fit_info, info = dk_iter.synthesize( ... eg["P"], ... eg["n_y"], @@ -490,7 +494,7 @@ def __init__( ... block_structure, ... ) >>> mu - 1.0360 + 1.03 """ super().__init__( controller_synthesis, @@ -507,7 +511,7 @@ def _get_fit_order( D_omega: np.ndarray, P: control.StateSpace, K: control.StateSpace, - block_structure: np.ndarray, + block_structure: List[uncertainty_structure.UncertaintyBlock], ) -> Optional[Union[int, np.ndarray]]: if iteration < len(self.fit_orders): return self.fit_orders[iteration] @@ -560,7 +564,11 @@ def __init__( ... max_fit_order=4, ... ) >>> omega = np.logspace(-3, 3, 61) - >>> block_structure = np.array([[1, 1], [1, 1], [2, 2]]) + >>> block_structure = [ + ... dkpy.ComplexFullBlock(1, 1), + ... dkpy.ComplexFullBlock(1, 1), + ... dkpy.ComplexFullBlock(2, 2), + ... ] >>> K, N, mu, d_scale_fit_info, info = dk_iter.synthesize( ... eg["P"], ... eg["n_y"], @@ -569,7 +577,7 @@ def __init__( ... block_structure, ... ) >>> mu - 1.0360 + 1.03 """ super().__init__( controller_synthesis, @@ -589,7 +597,7 @@ def _get_fit_order( D_omega: np.ndarray, P: control.StateSpace, K: control.StateSpace, - block_structure: np.ndarray, + block_structure: List[uncertainty_structure.UncertaintyBlock], ) -> Optional[Union[int, np.ndarray]]: # Check termination conditions if (self.max_iterations is not None) and (iteration >= self.max_iterations): @@ -668,7 +676,11 @@ def __init__( ... max_fit_order=4, ... ) >>> omega = np.logspace(-3, 3, 61) - >>> block_structure = np.array([[1, 1], [1, 1], [2, 2]]) + >>> block_structure = [ + ... dkpy.ComplexFullBlock(1, 1), + ... dkpy.ComplexFullBlock(1, 1), + ... dkpy.ComplexFullBlock(2, 2), + ... ] >>> K, N, mu, d_scale_fit_info, info = dk_iter.synthesize( ... eg["P"], ... eg["n_y"], @@ -686,14 +698,14 @@ def __init__( def _get_fit_order( self, - iteration, - omega, - mu_omega, - D_omega, - P, - K, - block_structure, - ): + iteration: int, + omega: np.ndarray, + mu_omega: np.ndarray, + D_omega: np.ndarray, + P: control.StateSpace, + K: control.StateSpace, + block_structure: List[uncertainty_structure.UncertaintyBlock], + ) -> Optional[Union[int, np.ndarray]]: d_info = [] for fit_order in range(self.max_fit_order + 1): D_fit, D_fit_inv = self.d_scale_fit.fit( @@ -773,7 +785,11 @@ def plot_mu( Create a ``IterResult`` object from fit data and plot ``mu`` >>> P, n_y, n_u, K = example_skogestad2006_p325 - >>> block_structure = np.array([[1, 1], [1, 1], [2, 2]]) + >>> block_structure = [ + ... dkpy.ComplexFullBlock(1, 1), + ... dkpy.ComplexFullBlock(1, 1), + ... dkpy.ComplexFullBlock(2, 2), + ... ] >>> omega = np.logspace(-3, 3, 61) >>> N = P.lft(K) >>> N_omega = N(1j * omega) @@ -864,7 +880,11 @@ def plot_D( Create a ``IterResult`` object from fit data and plot ``D`` >>> P, n_y, n_u, K = example_skogestad2006_p325 - >>> block_structure = np.array([[1, 1], [1, 1], [2, 2]]) + >>> block_structure = [ + ... dkpy.ComplexFullBlock(1, 1), + ... dkpy.ComplexFullBlock(1, 1), + ... dkpy.ComplexFullBlock(2, 2), + ... ] >>> omega = np.logspace(-3, 3, 61) >>> N = P.lft(K) >>> N_omega = N(1j * omega) @@ -885,7 +905,8 @@ def plot_D( ... ) >>> fig, ax = dkpy.plot_D(d_scale_fit_info) """ - mask = d_scale_fit._mask_from_block_structure(d_scale_info.block_structure) + block_structure = d_scale_info.block_structure + mask = d_scale_fit._generate_d_scale_mask(block_structure) # Create figure if not provided if ax is None: fig, ax = plt.subplots( @@ -962,6 +983,10 @@ def _augment_d_scales( Tuple[control.StateSpace, control.StateSpace] Augmented D-scales and inverse D-scales. """ - D_aug = control.append(D, utilities._tf_eye(n_y)) - D_aug_inv = control.append(D_inv, utilities._tf_eye(n_u)) + D_aug = control.ss( + control.append(D, control.TransferFunction([1], [1]) * np.eye(n_y)) + ) + D_aug_inv = control.ss( + control.append(D_inv, control.TransferFunction([1], [1]) * np.eye(n_u)) + ) return (D_aug, D_aug_inv) diff --git a/src/dkpy/structured_singular_value.py b/src/dkpy/structured_singular_value.py index 06463a5..eba8c7b 100644 --- a/src/dkpy/structured_singular_value.py +++ b/src/dkpy/structured_singular_value.py @@ -7,7 +7,7 @@ import abc import warnings -from typing import Any, Dict, Optional, Tuple +from typing import Any, Dict, Optional, Tuple, Union, List import cvxpy import joblib @@ -15,6 +15,7 @@ import scipy.linalg from . import utilities +from . import uncertainty_structure class StructuredSingularValue(metaclass=abc.ABCMeta): @@ -24,7 +25,9 @@ class StructuredSingularValue(metaclass=abc.ABCMeta): def compute_ssv( self, N_omega: np.ndarray, - block_structure: np.ndarray, + block_structure: Union[ + List[uncertainty_structure.UncertaintyBlock], List[List[int]], np.ndarray + ], ) -> Tuple[np.ndarray, np.ndarray, Dict[str, Any]]: """Compute structured singular value. @@ -32,22 +35,14 @@ def compute_ssv( ---------- N_omega : np.ndarray Closed-loop transfer function evaluated at each frequency. - block_structure : np.ndarray - 2D array with 2 columns and as many rows as uncertainty blocks - in Delta. The columns represent the number of rows and columns in - each uncertainty block. See [#mussv]_. - - Returns + block_structure : Union[List[uncertainty_structure.UncertaintyBlock], List[List[int]], np.ndarray] + Uncertainty block structure representation. Returns ------- Tuple[np.ndarray, np.ndarray, Dict[str, Any]] Structured singular value at each frequency, D-scales at each frequency, and solution information. If the structured singular value cannot be computed, the first two elements of the tuple are ``None``, but solution information is still returned. - - References - ---------- - .. [#mussv] https://www.mathworks.com/help/robust/ref/mussv.html """ raise NotImplementedError() @@ -63,7 +58,11 @@ class SsvLmiBisection(StructuredSingularValue): Structured singular value computation from [SP06]_ >>> P, n_y, n_u, K = example_skogestad2006_p325 - >>> block_structure = np.array([[1, 1], [1, 1], [2, 2]]) + >>> block_structure = [ + ... dkpy.ComplexFullBlock(1, 1), + ... dkpy.ComplexFullBlock(1, 1), + ... dkpy.ComplexFullBlock(2, 2), + ... ] >>> omega = np.logspace(-3, 3, 61) >>> N = P.lft(K) >>> N_omega = N(1j * omega) @@ -72,7 +71,7 @@ class SsvLmiBisection(StructuredSingularValue): ... block_structure, ... ) >>> float(np.max(mu_omega)) - 5.7726 + 5.77 """ def __init__( @@ -139,7 +138,9 @@ def __init__( def compute_ssv( self, N_omega: np.ndarray, - block_structure: np.ndarray, + block_structure: Union[ + List[uncertainty_structure.UncertaintyBlock], List[List[int]], np.ndarray + ], ) -> Tuple[np.ndarray, np.ndarray, Dict[str, Any]]: # Solver settings solver_params = ( @@ -165,10 +166,14 @@ def compute_ssv( constant_objective = True else: raise ValueError("`objective` must be `'minimize'` or `'constant'`.") + # Convert uncertainty block structure representation + block_structure = uncertainty_structure._convert_block_structure_representation( + block_structure + ) def _ssv_at_omega( N_omega: np.ndarray, - ) -> Tuple[float, np.ndarray, Dict[str, Any]]: + ) -> Tuple[Union[float, None], Union[np.ndarray, None], Dict[str, Any]]: """Compute the structured singular value at a given frequency. Split into its own function to allow parallelization over @@ -195,7 +200,7 @@ def _ssv_at_omega( info = {} # Get an optimization variable that shares its block structure with # the D-scalings - X = _variable_from_block_structure(block_structure) + X = _generate_ssv_variable(block_structure) # Set objective function if constant_objective: objective = cvxpy.Minimize(1) @@ -290,8 +295,7 @@ def _ssv_at_omega( info["size_metrics"] = [p.size_metrics for p in problems] info["results"] = results info["iterations"] = n_iterations - return None, None, info - # Save info + return None, None, info # Save info info["status"] = "Bisection succeeded." info["gammas"] = gammas info["solver_stats"] = [p.solver_stats for p in problems] @@ -319,58 +323,62 @@ def _ssv_at_omega( return mu, D_scales, info -def _variable_from_block_structure(block_structure: np.ndarray) -> cvxpy.Variable: - """Get optimization variable with specified block structure. +def _generate_ssv_variable( + block_structure: List[uncertainty_structure.UncertaintyBlock], +) -> cvxpy.Variable: + """Get structured singular value optimization variable for the block structure. Parameters ---------- - block_structure : np.ndarray - 2D array with 2 columns and as many rows as uncertainty blocks - in Delta. The columns represent the number of rows and columns in - each uncertainty block. See [#mussv]_. + block_structure : List[uncertainty_structure.UncertaintyBlock] + Uncertainty block structure representation. Returns ------- cvxpy.Variable CVXPY variable with specified block structure. - - References - ---------- - .. [#mussv] https://www.mathworks.com/help/robust/ref/mussv.html """ + num_blocks = len(block_structure) X_lst = [] - for i in range(block_structure.shape[0]): + for i in range(num_blocks): row = [] - for j in range(block_structure.shape[0]): + for j in range(num_blocks): + # Uncertainty blocks + block_i = block_structure[i] + block_j = block_structure[j] if i == j: # If on the block diagonal, insert variable - if block_structure[i, 0] <= 0: + if (not block_i.is_complex) and (not block_i.is_diagonal): + raise NotImplementedError( + "Real full perturbations are not supported." + ) + if (not block_i.is_complex) and (block_i.is_diagonal): raise NotImplementedError( - "Real perturbations are not yet supported." + "Real diagonal perturbations are not yet supported." ) - if block_structure[i, 1] <= 0: + if (block_i.is_complex) and (block_i.is_diagonal): raise NotImplementedError( - "Diagonal perturbations are not yet supported." + "Complex diagonal perturbations are not yet supported." ) - if block_structure[i, 0] != block_structure[i, 1]: + if (block_i.is_complex) and (not block_i.is_square): raise NotImplementedError( "Nonsquare perturbations are not yet supported." ) - if i == block_structure.shape[0] - 1: + if i == num_blocks - 1: # Last scaling is always identity - row.append(np.eye(block_structure[i, 0])) + row.append(np.eye(block_i.num_inputs)) else: # Every other scaling is either a scalar or a scalar # multiplied by identity - if block_structure[i, 0] == 1: + if block_i.num_inputs == 1: xi = cvxpy.Variable((1, 1), complex=True, name=f"x{i}") row.append(xi) else: xi = cvxpy.Variable(1, complex=True, name=f"x{i}") - row.append(xi * np.eye(block_structure[i, 0])) + row.append(xi * np.eye(block_i.num_inputs)) else: # If off the block diagonal, insert zeros - row.append(np.zeros((block_structure[i, 0], block_structure[j, 1]))) + row.append(np.zeros((block_i.num_inputs, block_j.num_inputs))) X_lst.append(row) X = cvxpy.bmat(X_lst) return X diff --git a/src/dkpy/uncertainty_structure.py b/src/dkpy/uncertainty_structure.py new file mode 100644 index 0000000..4d925f1 --- /dev/null +++ b/src/dkpy/uncertainty_structure.py @@ -0,0 +1,257 @@ +"""Uncertainty block.""" + +import numpy as np +import abc +from typing import List, Union + + +class UncertaintyBlock(metaclass=abc.ABCMeta): + """Generic uncertainty block.""" + + @property + @abc.abstractmethod + def num_inputs(self) -> int: + """Get number of inputs of the uncertainty block.""" + raise NotImplementedError() + + @property + @abc.abstractmethod + def num_outputs(self) -> int: + """Get number of output of the uncertainty block.""" + raise NotImplementedError() + + @property + @abc.abstractmethod + def is_diagonal(self) -> bool: + """Get boolean for diagonal uncertainty.""" + raise NotImplementedError() + + @property + @abc.abstractmethod + def is_complex(self) -> bool: + """Get boolean for complex uncertainty.""" + raise NotImplementedError() + + @property + @abc.abstractmethod + def is_square(self) -> bool: + """Get boolean for square uncertainty.""" + raise NotImplementedError() + + +class RealDiagonalBlock(UncertaintyBlock): + """Real-valued diagonal uncertainty block.""" + + def __init__(self, num_channels: int): + """Instantiate :class:`RealDiagonalBlock`. + + Parameters + ---------- + num_channels : int + Number of inputs/outputs to the uncertainty block. + + Raises + ------ + ValueError + If ``num_channels`` is not greater than zero. + """ + if num_channels <= 0: + raise ValueError("`num_channels` must be greater than 0.") + + self._num_channels = num_channels + + @property + def num_inputs(self) -> int: + return self._num_channels + + @property + def num_outputs(self) -> int: + """Get number of output of the uncertainty block.""" + return self._num_channels + + @property + def is_diagonal(self) -> bool: + """Get boolean for diagonal uncertainty.""" + return True + + @property + def is_complex(self) -> bool: + """Get boolean for complex uncertainty.""" + return False + + @property + def is_square(self) -> bool: + """Get boolean for square uncertainty.""" + return True + + +class ComplexDiagonalBlock(UncertaintyBlock): + """Complex-valued diagonal uncertainty block.""" + + def __init__(self, num_channels: int): + """Instantiate :class:`ComplexDiagonalBlock`. + + Parameters + ---------- + num_channels : int + Number of inputs/outputs to the uncertainty block. + + Raises + ------ + ValueError + If ``num_channels`` is not greater than zero. + """ + if num_channels <= 0: + raise ValueError("`num_channels` must be greater than 0.") + + self._num_channels = num_channels + + @property + def num_inputs(self) -> int: + return self._num_channels + + @property + def num_outputs(self) -> int: + """Get number of output of the uncertainty block.""" + return self._num_channels + + @property + def is_diagonal(self) -> bool: + """Get boolean for diagonal uncertainty.""" + return True + + @property + def is_complex(self) -> bool: + """Get boolean for complex uncertainty.""" + return True + + @property + def is_square(self) -> bool: + """Get boolean for square uncertainty.""" + return True + + +class ComplexFullBlock(UncertaintyBlock): + """Complex-valued full uncertainty block.""" + + def __init__(self, num_inputs: int, num_outputs: int): + """Instantiate :class:`ComplexFullBlock`. + + Parameters + ---------- + num_inputs : int + Number of inputs (columns) to the uncertainty block. + num_outputs : int + Number of outputs (rows) to the uncertainty block. + + Raises + ------ + ValueError + If ``num_inputs`` is not greater than zero. + ValueError + If ``num_outputs`` is not greater than zero. + """ + + if num_inputs <= 0: + raise ValueError("`num_inputs` must be greater than 0.") + if num_outputs <= 0: + raise ValueError("`num_outputs` must be greater than 0.") + + self._num_inputs = num_inputs + self._num_outputs = num_outputs + + @property + def num_inputs(self) -> int: + return self._num_inputs + + @property + def num_outputs(self) -> int: + """Get number of output of the uncertainty block.""" + return self._num_outputs + + @property + def is_diagonal(self) -> bool: + """Get boolean for diagonal uncertainty.""" + return False + + @property + def is_complex(self) -> bool: + """Get boolean for complex uncertainty.""" + return True + + @property + def is_square(self) -> bool: + """Get boolean for square uncertainty.""" + return self._num_inputs == self._num_outputs + + +def _convert_block_structure_representation( + block_structure: Union[List[UncertaintyBlock], List[List[int]], np.ndarray], +) -> List[UncertaintyBlock]: + """ + Convert the uncertainty block description into the object-oriented description. + + The MATLAB uncertainty block structure description is given by + - block_structure_matlab[i,:] = [-r 0]: i-th block is an r-by-r repeated, diagonal + real scalar perturbation. + - block_structure_matlab[i,:] = [r 0]: i-th block is an r-by-r repeated, diagonal + complex scalar perturbation. + - block_structure_matlab[i,:] = [r c]: i-th block is an r-by-c complex full-block + perturbation. + For additional reference, see [#matlab_block_struct]_. + + Parameters + ---------- + block_structure: Union[List[UncertaintyBlock], List[int]], + Uncertainty block structure representation. + + Returns + ------- + List[UncertaintyBlock] + Object-oriented uncertainty block structure representation. + + Raises + ------ + ValueError + Block structure cannot be converted to an array. + ValueError + Block structure does not have the correct number of columns for the + MATLAB uncertainty structure representation. + ValueError + Uncertainty block does not correspond to MATLAB standard description. + + References + ---------- + .. [#matlab_block_struct] https://www.mathworks.com/help/robust/ref/mussv.html + """ + + # Block structure is in OOP representation + is_oop_representation = all( + isinstance(block, UncertaintyBlock) for block in block_structure + ) + if is_oop_representation: + return block_structure + # Block structure is in MATLAB form + block_structure_matlab = np.array(block_structure) + if block_structure_matlab.shape[1] != 2: + raise ValueError( + "The `block_structure` array does not have the correct dimemsions " + "for the MATLAB uncertainty specification. It must have 2 columns, " + f"but it has {block_structure_matlab.shape[1]}." + ) + block_structure_oop = [] + for block_matlab in block_structure_matlab: + if (block_matlab[0] < 0) and (block_matlab[1] == 0): + block_structure_oop.append(RealDiagonalBlock(np.abs(block_matlab[0]))) + elif (block_matlab[0] > 0) and (block_matlab[1] == 0): + block_structure_oop.append(ComplexDiagonalBlock(block_matlab[0])) + elif (block_matlab[0] > 0) and (block_matlab[1] > 0): + block_structure_oop.append( + ComplexFullBlock(block_matlab[1], block_matlab[0]) + ) + else: + raise ValueError( + "The uncertainty block array does not conform to the MATLAB standard." + ) + + return block_structure_oop diff --git a/src/dkpy/utilities.py b/src/dkpy/utilities.py index 46ca391..6ab2f6f 100644 --- a/src/dkpy/utilities.py +++ b/src/dkpy/utilities.py @@ -5,15 +5,10 @@ "example_skogestad2006_p325", "_ensure_tf", "_tf_close_coeff", - "_tf_combine", - "_tf_split", - "_tf_eye", - "_tf_zeros", - "_tf_ones", "_auto_lmi_strictness", ] -from typing import Any, Dict, List, Union +from typing import Any, Dict, Union import control import cvxpy @@ -229,236 +224,6 @@ def _ensure_tf( return tf -def _tf_combine( - tf_array: List[List[Union[ArrayLike, control.TransferFunction]]], -) -> control.TransferFunction: - """Combine array-like of transfer functions into MIMO transfer function. - - Parameters - ---------- - tf_array : List[List[Union[ArrayLike, control.TransferFunction]]] - Transfer matrix represented as a two-dimensional array or list-of-lists - containing ``TransferFunction`` objects. The ``TransferFunction`` - objects can have multiple outputs and inputs, as long as the dimensions - are compatible. - - Returns - ------- - control.TransferFunction - Transfer matrix represented as a single MIMO ``TransferFunction`` - object. - - Raises - ------ - ValueError - If timesteps of transfer functions do not match. - ValueError - If ``tf_array`` has incorrect dimensions. - - Examples - -------- - Combine two transfer functions - - >>> s = control.TransferFunction.s - >>> dkpy._tf_combine([ - ... [1 / (s + 1)], - ... [s / (s + 2)], - ... ]) - TransferFunction([[array([1])], [array([1, 0])]], [[array([1, 1])], [array([1, 2])]]) - - Combine NumPy arrays with transfer functions - - >>> dkpy._tf_combine([ - ... [np.eye(2), np.zeros((2, 1))], - ... [np.zeros((1, 2)), control.TransferFunction([1], [1, 0])], - ... ]) - TransferFunction([[array([1.]), array([0.]), array([0.])], - [array([0.]), array([1.]), array([0.])], - [array([0.]), array([0.]), array([1])]], - [[array([1.]), array([1.]), array([1.])], - [array([1.]), array([1.]), array([1.])], - [array([1.]), array([1.]), array([1, 0])]]) - """ - # Find common timebase or raise error - dt_list = [] - try: - for row in tf_array: - for tf in row: - dt_list.append(getattr(tf, "dt", None)) - except OSError: - raise ValueError("`tf_array` has too few dimensions.") - dt_set = set(dt_list) - dt_set.discard(None) - if len(dt_set) > 1: - raise ValueError(f"Timesteps of transfer functions are mismatched: {dt_set}") - elif len(dt_set) == 0: - dt = None - else: - dt = dt_set.pop() - # Convert all entries to transfer function objects - ensured_tf_array = [] - for row in tf_array: - ensured_row = [] - for tf in row: - ensured_row.append(_ensure_tf(tf, dt)) - ensured_tf_array.append(ensured_row) - # Iterate over - num = [] - den = [] - for row_index, row in enumerate(ensured_tf_array): - for j_out in range(row[0].noutputs): - num_row = [] - den_row = [] - for col in row: - if col.noutputs != row[0].noutputs: - raise ValueError( - f"Mismatched number of transfer function outputs in row {row_index}." - ) - for j_in in range(col.ninputs): - num_row.append(col.num[j_out][j_in]) - den_row.append(col.den[j_out][j_in]) - num.append(num_row) - den.append(den_row) - G_tf = control.TransferFunction(num, den, dt=dt) - return G_tf - - -def _tf_split(tf: control.TransferFunction) -> np.ndarray: - """Split MIMO transfer function into NumPy array of SISO tranfer functions. - - Parameters - ---------- - tf : control.TransferFunction - MIMO transfer function to split. - - Returns - ------- - np.ndarray - NumPy array of SISO transfer functions. - - Examples - -------- - Split a MIMO transfer function - - >>> G = control.TransferFunction( - ... [ - ... [[87.8], [-86.4]], - ... [[108.2], [-109.6]], - ... ], - ... [ - ... [[1, 1], [1, 1]], - ... [[1, 1], [1, 1]], - ... ], - ... ) - >>> dkpy._tf_split(G) - array([[TransferFunction(array([87.8]), array([1, 1])), - TransferFunction(array([-86.4]), array([1, 1]))], - [TransferFunction(array([108.2]), array([1, 1])), - TransferFunction(array([-109.6]), array([1, 1]))]], dtype=object) - """ - tf_split_lst = [] - for i_out in range(tf.noutputs): - row = [] - for i_in in range(tf.ninputs): - row.append( - control.TransferFunction( - tf.num[i_out][i_in], - tf.den[i_out][i_in], - dt=tf.dt, - ) - ) - tf_split_lst.append(row) - tf_split = np.array(tf_split_lst, dtype=object) - return tf_split - - -def _tf_eye( - n: int, - dt: Union[None, bool, float] = None, -) -> control.TransferFunction: - """Transfer function identity matrix. - - Parameters - ---------- - n : int - Dimension. - dt : Union[None, bool, float] - Timestep (s). Based on the ``control`` package, ``True`` indicates a - discrete-time system with unspecified timestep, ``0`` indicates a - continuous-time system, and ``None`` indicates a continuous- or - discrete-time system with unspecified timestep. - - Returns - ------- - control.TransferFunction - Identity transfer matrix. - """ - num = np.eye(n).reshape(n, n, 1) - den = np.ones((n, n, 1)) - eye = control.TransferFunction(num, den, dt=dt) - return eye - - -def _tf_zeros( - m: int, - n: int, - dt: Union[None, bool, float] = None, -) -> control.TransferFunction: - """Transfer function matrix of zeros. - - Parameters - ---------- - m : int - First dimension. - n : int - Second dimension. - dt : Union[None, bool, float] - Timestep (s). Based on the ``control`` package, ``True`` indicates a - discrete-time system with unspecified timestep, ``0`` indicates a - continuous-time system, and ``None`` indicates a continuous- or - discrete-time system with unspecified timestep. - - Returns - ------- - control.TransferFunction - Identity transfer matrix. - """ - num = np.zeros((m, n, 1)) - den = np.ones((m, n, 1)) - zeros = control.TransferFunction(num, den, dt=dt) - return zeros - - -def _tf_ones( - m: int, - n: int, - dt: Union[None, bool, float] = None, -) -> control.TransferFunction: - """Transfer matrix of ones. - - Parameters - ---------- - m : int - First dimension. - n : int - Second dimension. - dt : Union[None, bool, float] - Timestep (s). Based on the ``control`` package, ``True`` indicates a - discrete-time system with unspecified timestep, ``0`` indicates a - continuous-time system, and ``None`` indicates a continuous- or - discrete-time system with unspecified timestep. - - Returns - ------- - control.TransferFunction - Identity transfer matrix. - """ - num = np.ones((m, n, 1)) - den = np.ones((m, n, 1)) - zeros = control.TransferFunction(num, den, dt=dt) - return zeros - - def _auto_lmi_strictness( solver_params: Dict[str, Any], scale: float = 10, diff --git a/tests/test_d_scale_fit.py b/tests/test_d_scale_fit.py index 53017ae..2bc3a96 100644 --- a/tests/test_d_scale_fit.py +++ b/tests/test_d_scale_fit.py @@ -111,6 +111,90 @@ class TestTfFitSlicot: None, 1e-2, ), + ( + np.logspace(-2, 2, 100), + control.TransferFunction( + [ + [ + [1, 1], + [0], + ], + [ + [0], + [1], + ], + ], + [ + [ + [1, 10], + [1], + ], + [ + [1], + [1], + ], + ], + ), + 1, + [dkpy.ComplexFullBlock(1, 1), dkpy.ComplexFullBlock(1, 1)], + 1e-2, + ), + ( + np.logspace(-2, 2, 100), + control.TransferFunction( + [ + [ + [1, 1], + [0], + ], + [ + [0], + [1], + ], + ], + [ + [ + [1, 10], + [1], + ], + [ + [1], + [1], + ], + ], + ), + np.diag([1, 0]), + [dkpy.ComplexFullBlock(1, 1), dkpy.ComplexFullBlock(1, 1)], + 1e-2, + ), + ( + np.logspace(-2, 2, 100), + control.TransferFunction( + [ + [ + [1], + [0], + ], + [ + [0], + [1], + ], + ], + [ + [ + [1], + [1], + ], + [ + [1], + [1], + ], + ], + ), + 1, + [dkpy.ComplexFullBlock(2, 2)], + 1e-2, + ), ( np.logspace(-2, 2, 100), control.TransferFunction( @@ -281,14 +365,12 @@ def test_tf_fit_slicot_error(self, omega, tf, order, block_structure): ) -class TestMaskFromBlockStructure: - """Test :func:`_mask_from_block_strucure`.""" - +class TestGenerateDScaleMask: @pytest.mark.parametrize( "block_structure, mask_exp", [ ( - np.array([[1, 1], [1, 1]]), + [dkpy.ComplexFullBlock(1, 1), dkpy.ComplexFullBlock(1, 1)], np.array( [ [-1, 0], @@ -298,7 +380,7 @@ class TestMaskFromBlockStructure: ), ), ( - np.array([[2, 2], [1, 1]]), + [dkpy.ComplexFullBlock(2, 2), dkpy.ComplexFullBlock(1, 1)], np.array( [ [-1, 0, 0], @@ -309,7 +391,7 @@ class TestMaskFromBlockStructure: ), ), ( - np.array([[1, 1], [2, 2]]), + [dkpy.ComplexFullBlock(1, 1), dkpy.ComplexFullBlock(2, 2)], np.array( [ [-1, 0, 0], @@ -321,9 +403,9 @@ class TestMaskFromBlockStructure: ), ], ) - def test_mask_from_block_structure(self, block_structure, mask_exp): + def test_generate_d_scale_mask(self, block_structure, mask_exp): """Test :func:`_mask_from_block_strucure`.""" - mask = dkpy.d_scale_fit._mask_from_block_structure(block_structure) + mask = dkpy.d_scale_fit._generate_d_scale_mask(block_structure) np.testing.assert_allclose(mask_exp, mask) diff --git a/tests/test_dk_iteration.py b/tests/test_dk_iteration.py index d42025a..5a3c14f 100644 --- a/tests/test_dk_iteration.py +++ b/tests/test_dk_iteration.py @@ -14,29 +14,25 @@ class TestAugmentDScales: "D, D_inv, n_y, n_u, D_aug_exp, D_aug_inv_exp", [ ( - dkpy.utilities._tf_eye(2), - dkpy.utilities._tf_eye(2), + control.TransferFunction([1], [1]) * np.eye(2), + control.TransferFunction([1], [1]) * np.eye(2), 2, 1, - dkpy.utilities._tf_eye(4), - dkpy.utilities._tf_eye(3), + control.TransferFunction([1], [1]) * np.eye(4), + control.TransferFunction([1], [1]) * np.eye(3), ), ( - 0.5 * dkpy.utilities._tf_eye(2), - 0.8 * dkpy.utilities._tf_eye(2), + control.TransferFunction([0.5], [1]) * np.eye(2), + control.TransferFunction([0.8], [1]) * np.eye(2), 2, 1, - control.ss2tf( - control.append( - 0.5 * dkpy.utilities._tf_eye(2), - dkpy.utilities._tf_eye(2), - ) + control.append( + control.TransferFunction([0.5], [1]) * np.eye(2), + control.TransferFunction([1], [1]) * np.eye(2), ), - control.ss2tf( - control.append( - 0.8 * dkpy.utilities._tf_eye(2), - dkpy.utilities._tf_eye(1), - ) + control.append( + control.TransferFunction([0.8], [1]) * np.eye(2), + control.TransferFunction([1], [1]) * np.eye(1), ), ), ], diff --git a/tests/test_structured_singular_value.py b/tests/test_structured_singular_value.py index 25ae08a..6d34c2e 100644 --- a/tests/test_structured_singular_value.py +++ b/tests/test_structured_singular_value.py @@ -7,14 +7,12 @@ import dkpy -class TestVariableFromBlockStructure: - """Test :func:`_variable_from_block_structure`.""" - +class TestGenerateSsvVariable: @pytest.mark.parametrize( "block_structure, variable_exp", [ ( - np.array([[1, 1], [2, 2]]), + [dkpy.ComplexFullBlock(1, 1), dkpy.ComplexFullBlock(2, 2)], cvxpy.bmat( [ [ @@ -29,7 +27,11 @@ class TestVariableFromBlockStructure: ), ), ( - np.array([[1, 1], [2, 2], [1, 1]]), + [ + dkpy.ComplexFullBlock(1, 1), + dkpy.ComplexFullBlock(2, 2), + dkpy.ComplexFullBlock(1, 1), + ], cvxpy.bmat( [ [ @@ -52,10 +54,10 @@ class TestVariableFromBlockStructure: ), ], ) - def test_variable_from_block_structure(self, block_structure, variable_exp): + def test_generate_ssv_variable(self, block_structure, variable_exp): """Test :func:`_variable_from_block_structure`.""" - variable = dkpy.structured_singular_value._variable_from_block_structure( - block_structure, + variable = dkpy.structured_singular_value._generate_ssv_variable( + block_structure ) assert variable.ndim == variable_exp.ndim assert variable.shape == variable_exp.shape diff --git a/tests/test_uncertainty_structure.py b/tests/test_uncertainty_structure.py new file mode 100644 index 0000000..45b576a --- /dev/null +++ b/tests/test_uncertainty_structure.py @@ -0,0 +1,151 @@ +"""Test :mod:`uncertainty_structure`.""" + +import numpy as np +import cvxpy +import pytest + +import dkpy + + +class TestRealDiagonalBlock: + """Test :class:`RealDiagonalBlock`.""" + + def test_real_diagonal_block(self): + """Test :class:`RealDiagonalBlock`.""" + block = dkpy.RealDiagonalBlock(5) + + assert block.num_inputs == 5 + assert block.num_outputs == 5 + assert block.is_diagonal + assert not block.is_complex + + +class TestComplexDiagonalBlock: + """Test :class:`ComplexDiagonalBlock`.""" + + def test_complex_diagonal_block(self): + """Test :class:`ComplexDiagonalBlock`.""" + block = dkpy.ComplexDiagonalBlock(5) + assert block.num_inputs == 5 + assert block.num_outputs == 5 + assert block.is_diagonal + assert block.is_complex + + +class TestComplexFullBlock: + """Test :class:`ComplexFullBlock`.""" + + def test_complex_full_block(self): + """Test :class:`ComplexFullBlock`.""" + block = dkpy.ComplexFullBlock(5, 10) + assert block.num_inputs == 5 + assert block.num_outputs == 10 + assert not block.is_diagonal + assert block.is_complex + + +class TestUncertaintyBlockStructure: + """Test :class:`UncertaintyBlockStructure""" + + @pytest.mark.parametrize( + "block_structure, num_inputs_list, num_outputs_list," + " is_diagonal_list, is_complex_list", + [ + ( + np.array([[2, 2], [4, 4], [-3, 0]]), + [2, 4, 3], + [2, 4, 3], + [False, False, True], + [True, True, False], + ), + ( + np.array([[-3, 0], [6, 0], [1, 2]]), + [3, 6, 2], + [3, 6, 1], + [True, True, False], + [False, True, True], + ), + ( + np.array([[3, 6], [-1, 0], [5, 1], [10, 0]]), + [6, 1, 1, 10], + [3, 1, 5, 10], + [False, True, False, True], + [True, False, True, True], + ), + ( + [[2, 2], [4, 4], [-3, 0]], + [2, 4, 3], + [2, 4, 3], + [False, False, True], + [True, True, False], + ), + ( + [[-3, 0], [6, 0], [1, 2]], + [3, 6, 2], + [3, 6, 1], + [True, True, False], + [False, True, True], + ), + ( + [[3, 6], [-1, 0], [5, 1], [10, 0]], + [6, 1, 1, 10], + [3, 1, 5, 10], + [False, True, False, True], + [True, False, True, True], + ), + ( + [ + dkpy.ComplexFullBlock(2, 2), + dkpy.ComplexFullBlock(4, 4), + dkpy.RealDiagonalBlock(3), + ], + [2, 4, 3], + [2, 4, 3], + [False, False, True], + [True, True, False], + ), + ( + [ + dkpy.RealDiagonalBlock(3), + dkpy.ComplexDiagonalBlock(6), + dkpy.ComplexFullBlock(2, 1), + ], + [3, 6, 2], + [3, 6, 1], + [True, True, False], + [False, True, True], + ), + ( + [ + dkpy.ComplexFullBlock(6, 3), + dkpy.RealDiagonalBlock(1), + dkpy.ComplexFullBlock(1, 5), + dkpy.ComplexDiagonalBlock(10), + ], + [6, 1, 1, 10], + [3, 1, 5, 10], + [False, True, False, True], + [True, False, True, True], + ), + ], + ) + def test_convert_block_structure_representation( + self, + block_structure, + num_inputs_list, + num_outputs_list, + is_diagonal_list, + is_complex_list, + ): + block_structure = ( + dkpy.uncertainty_structure._convert_block_structure_representation( + block_structure + ) + ) + + for idx_block in range(len(block_structure)): + block = block_structure[idx_block] + assert block.num_inputs == num_inputs_list[idx_block] + assert block.num_outputs == num_outputs_list[idx_block] + assert block.is_diagonal == is_diagonal_list[idx_block] + assert block.is_complex == is_complex_list[idx_block] diff --git a/tests/test_utilities.py b/tests/test_utilities.py index 8d0bdcd..fa70659 100644 --- a/tests/test_utilities.py +++ b/tests/test_utilities.py @@ -201,487 +201,6 @@ def test_error_ensure(self, arraylike_or_tf, dt, exception): dkpy._ensure_tf(arraylike_or_tf, dt) -class TestTfCombineSplit: - """Test :func:`_tf_combine` and :func:`_tf_split`.""" - - @pytest.mark.parametrize( - "tf_array, tf", - [ - # Continuous-time - ( - [ - [control.TransferFunction([1], [1, 1])], - [control.TransferFunction([2], [1, 0])], - ], - control.TransferFunction( - [ - [[1]], - [[2]], - ], - [ - [[1, 1]], - [[1, 0]], - ], - ), - ), - # Discrete-time - ( - [ - [control.TransferFunction([1], [1, 1], dt=1)], - [control.TransferFunction([2], [1, 0], dt=1)], - ], - control.TransferFunction( - [ - [[1]], - [[2]], - ], - [ - [[1, 1]], - [[1, 0]], - ], - dt=1, - ), - ), - # Scalar - ( - [ - [2], - [control.TransferFunction([2], [1, 0])], - ], - control.TransferFunction( - [ - [[2]], - [[2]], - ], - [ - [[1]], - [[1, 0]], - ], - ), - ), - # Matrix - ( - [ - [np.eye(3)], - [ - control.TransferFunction( - [ - [[2], [0], [3]], - [[1], [2], [3]], - ], - [ - [[1], [1], [1]], - [[1], [1], [1]], - ], - ) - ], - ], - control.TransferFunction( - [ - [[1], [0], [0]], - [[0], [1], [0]], - [[0], [0], [1]], - [[2], [0], [3]], - [[1], [2], [3]], - ], - [ - [[1], [1], [1]], - [[1], [1], [1]], - [[1], [1], [1]], - [[1], [1], [1]], - [[1], [1], [1]], - ], - ), - ), - # Inhomogeneous - ( - [ - [np.eye(3)], - [ - control.TransferFunction( - [ - [[2], [0]], - [[1], [2]], - ], - [ - [[1], [1]], - [[1], [1]], - ], - ), - control.TransferFunction( - [ - [[3]], - [[3]], - ], - [ - [[1]], - [[1]], - ], - ), - ], - ], - control.TransferFunction( - [ - [[1], [0], [0]], - [[0], [1], [0]], - [[0], [0], [1]], - [[2], [0], [3]], - [[1], [2], [3]], - ], - [ - [[1], [1], [1]], - [[1], [1], [1]], - [[1], [1], [1]], - [[1], [1], [1]], - [[1], [1], [1]], - ], - ), - ), - # Discrete-time - ( - [ - [2], - [control.TransferFunction([2], [1, 0], dt=0.1)], - ], - control.TransferFunction( - [ - [[2]], - [[2]], - ], - [ - [[1]], - [[1, 0]], - ], - dt=0.1, - ), - ), - ], - ) - def test_combine(self, tf_array, tf): - """Test combining transfer functions.""" - tf_combined = dkpy._tf_combine(tf_array) - assert dkpy._tf_close_coeff(tf_combined, tf) - - @pytest.mark.parametrize( - "tf_array, tf", - [ - ( - np.array( - [ - [control.TransferFunction([1], [1, 1])], - ], - dtype=object, - ), - control.TransferFunction( - [ - [[1]], - ], - [ - [[1, 1]], - ], - ), - ), - ( - np.array( - [ - [control.TransferFunction([1], [1, 1])], - [control.TransferFunction([2], [1, 0])], - ], - dtype=object, - ), - control.TransferFunction( - [ - [[1]], - [[2]], - ], - [ - [[1, 1]], - [[1, 0]], - ], - ), - ), - ( - np.array( - [ - [control.TransferFunction([1], [1, 1], dt=1)], - [control.TransferFunction([2], [1, 0], dt=1)], - ], - dtype=object, - ), - control.TransferFunction( - [ - [[1]], - [[2]], - ], - [ - [[1, 1]], - [[1, 0]], - ], - dt=1, - ), - ), - ( - np.array( - [ - [control.TransferFunction([2], [1], dt=0.1)], - [control.TransferFunction([2], [1, 0], dt=0.1)], - ], - dtype=object, - ), - control.TransferFunction( - [ - [[2]], - [[2]], - ], - [ - [[1]], - [[1, 0]], - ], - dt=0.1, - ), - ), - ], - ) - def test_split(self, tf_array, tf): - """Test splitting transfer functions.""" - tf_split = dkpy._tf_split(tf) - # Test entry-by-entry - for i in range(tf_split.shape[0]): - for j in range(tf_split.shape[1]): - assert dkpy._tf_close_coeff( - tf_split[i, j], - tf_array[i, j], - ) - # Test combined - assert dkpy._tf_close_coeff( - dkpy._tf_combine(tf_split), - dkpy._tf_combine(tf_array), - ) - - @pytest.mark.parametrize( - "tf_array, exception", - [ - # Wrong timesteps - ( - [ - [control.TransferFunction([1], [1, 1], 0.1)], - [control.TransferFunction([2], [1, 0], 0.2)], - ], - ValueError, - ), - ( - [ - [control.TransferFunction([1], [1, 1], 0.1)], - [control.TransferFunction([2], [1, 0], 0)], - ], - ValueError, - ), - # Too few dimensions - ( - [ - control.TransferFunction([1], [1, 1]), - control.TransferFunction([2], [1, 0]), - ], - ValueError, - ), - # Too many dimensions - ( - [ - [[control.TransferFunction([1], [1, 1], 0.1)]], - [[control.TransferFunction([2], [1, 0], 0)]], - ], - ValueError, - ), - # Incompatible dimensions - ( - [ - [ - control.TransferFunction( - [ - [ - [1], - ] - ], - [ - [ - [1, 1], - ] - ], - ), - control.TransferFunction( - [ - [[2], [1]], - [[1], [3]], - ], - [ - [[1, 0], [1, 0]], - [[1, 0], [1, 0]], - ], - ), - ], - ], - ValueError, - ), - ( - [ - [ - control.TransferFunction( - [ - [[2], [1]], - [[1], [3]], - ], - [ - [[1, 0], [1, 0]], - [[1, 0], [1, 0]], - ], - ), - control.TransferFunction( - [ - [ - [1], - ] - ], - [ - [ - [1, 1], - ] - ], - ), - ], - ], - ValueError, - ), - ], - ) - def test_error_combine(self, tf_array, exception): - """Test error cases.""" - with pytest.raises(exception): - dkpy._tf_combine(tf_array) - - -class TestTfEye: - """Test :func:`_tf_eye`.""" - - @pytest.mark.parametrize( - "n, dt, tf_exp", - [ - ( - 1, - None, - control.TransferFunction([1], [1], dt=None), - ), - ( - 2, - 0, - control.TransferFunction( - [ - [[1], [0]], - [[0], [1]], - ], - [ - [[1], [1]], - [[1], [1]], - ], - dt=0, - ), - ), - ( - 3, - 1e-3, - control.TransferFunction( - [ - [[1], [0], [0]], - [[0], [1], [0]], - [[0], [0], [1]], - ], - [ - [[1], [1], [1]], - [[1], [1], [1]], - [[1], [1], [1]], - ], - dt=1e-3, - ), - ), - ], - ) - def test_tf_eye(self, n, dt, tf_exp): - """Test :func:`_tf_eye`.""" - tf = dkpy._tf_eye(n, dt) - assert dkpy._tf_close_coeff(tf, tf_exp) - - -class TestTfZeros: - """Test :func:`_tf_zeros`.""" - - @pytest.mark.parametrize( - "m, n, dt, tf_exp", - [ - ( - 1, - 1, - None, - control.TransferFunction([0], [1], dt=None), - ), - ( - 2, - 3, - None, - control.TransferFunction( - [ - [[0], [0], [0]], - [[0], [0], [0]], - ], - [ - [[1], [1], [1]], - [[1], [1], [1]], - ], - dt=None, - ), - ), - ], - ) - def test_tf_zeros(self, m, n, dt, tf_exp): - """Test :func:`_tf_zeros`.""" - tf = dkpy._tf_zeros(m, n, dt) - assert dkpy._tf_close_coeff(tf, tf_exp) - - -class TestTfOnes: - """Test :func:`_tf_ones`.""" - - @pytest.mark.parametrize( - "m, n, dt, tf_exp", - [ - ( - 1, - 1, - None, - control.TransferFunction([1], [1], dt=None), - ), - ( - 2, - 3, - None, - control.TransferFunction( - [ - [[1], [1], [1]], - [[1], [1], [1]], - ], - [ - [[1], [1], [1]], - [[1], [1], [1]], - ], - dt=None, - ), - ), - ], - ) - def test_tf_ones(self, m, n, dt, tf_exp): - """Test :func:`_tf_ones`.""" - tf = dkpy._tf_ones(m, n, dt) - assert dkpy._tf_close_coeff(tf, tf_exp) - - class TestAutoLmiStrictness: """Test :func:`_auto_lmi_strictness`."""