Merge pull request #13 from decargroup/dev

Fix `matplotlib` interaction with `joblib`

CITATION.cff

@@ -10,5 +10,5 @@ authors:
     orcid: "https://orcid.org/0000-0002-1987-9268"
     affiliation: "McGill University"
 title: "decargroup/dkpy"
-version: v0.1.4
+version: v0.1.5
 url: "https://github.com/decargroup/dkpy"

README.rst

@@ -69,6 +69,6 @@ If you use this software in your research, please cite it as below or see
         title={{decargroup/dkpy}},
         url={https://github.com/decargroup/dkpy},
         author={Steven Dahdah and James Richard Forbes},
-        version = {{v0.1.4}},
+        version = {{v0.1.5}},
         year={2024},
     }

doc/conf.py

@@ -9,8 +9,8 @@
 project = "dkpy"
 copyright = "2024, Steven Dahdah and James Richard Forbes"
 author = "Steven Dahdah and James Richard Forbes"
-version = "0.1.4"
-release = "0.1.4"
+version = "0.1.5"
+release = "0.1.5"
 
 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration

examples/3_example_dk_iter_order_callback.py

@@ -0,0 +1,154 @@
+"""D-K iteration with fixed number of iterations and fit order."""
+
+import control
+import numpy as np
+from matplotlib import pyplot as plt
+
+import dkpy
+
+
+def example_dk_iter_fixed_order():
+    """D-K iteration with fixed number of iterations and fit order."""
+    # Plant
+    G0 = np.array(
+        [
+            [87.8, -86.4],
+            [108.2, -109.6],
+        ]
+    )
+    G = control.append(
+        control.TransferFunction([1], [75, 1]),
+        control.TransferFunction([1], [75, 1]),
+    ) * control.TransferFunction(
+        G0.reshape(2, 2, 1),
+        np.ones((2, 2, 1)),
+    )
+    # Weights
+    Wp = 0.5 * control.append(
+        control.TransferFunction([10, 1], [10, 1e-5]),
+        control.TransferFunction([10, 1], [10, 1e-5]),
+    )
+    Wi = control.append(
+        control.TransferFunction([1, 0.2], [0.5, 1]),
+        control.TransferFunction([1, 0.2], [0.5, 1]),
+    )
+    G.name = "G"
+    Wp.name = "Wp"
+    Wi.name = "Wi"
+    sum_w = control.summing_junction(
+        inputs=["u_w", "u_G"],
+        dimension=2,
+        name="sum_w",
+    )
+    sum_del = control.summing_junction(
+        inputs=["u_del", "u_u"],
+        dimension=2,
+        name="sum_del",
+    )
+    split = control.summing_junction(
+        inputs=["u"],
+        dimension=2,
+        name="split",
+    )
+    P = control.interconnect(
+        syslist=[G, Wp, Wi, sum_w, sum_del, split],
+        connections=[
+            ["G.u", "sum_del.y"],
+            ["sum_del.u_u", "split.y"],
+            ["sum_w.u_G", "G.y"],
+            ["Wp.u", "sum_w.y"],
+            ["Wi.u", "split.y"],
+        ],
+        inplist=["sum_del.u_del", "sum_w.u_w", "split.u"],
+        outlist=["Wi.y", "Wp.y", "-sum_w.y"],
+    )
+    # Dimensions
+    n_y = 2
+    n_u = 2
+
+    def callback(
+        dk_iteration,
+        iteration,
+        omega,
+        mu_omega,
+        D_omega,
+        P,
+        K,
+        block_structure,
+    ):
+        d_info = []
+        for fit_order in range(5):
+            D_fit, D_fit_inv = dk_iteration.transfer_function_fit.fit(
+                omega,
+                D_omega,
+                order=fit_order,
+                block_structure=block_structure,
+            )
+            d_info.append(
+                dkpy.DScaleFitInfo.create_from_fit(
+                    omega,
+                    mu_omega,
+                    D_omega,
+                    P,
+                    K,
+                    D_fit,
+                    D_fit_inv,
+                    block_structure,
+                )
+            )
+        fig, ax = plt.subplots()
+        d_info[0].plot_mu(ax=ax, plot_kw=dict(label="true"), hide="mu_fit_omega")
+        for i, ds in enumerate(d_info):
+            ds.plot_mu(ax=ax, plot_kw=dict(label=f"order={i}"), hide="mu_omega")
+        plt.show()
+        selected_order = int(input("Selected order: "))
+        done = input("Done? (y/N): ") == "y"
+        return (selected_order, done)
+
+    dk_iter = dkpy.DkIterOrderCallback(
+        controller_synthesis=dkpy.HinfSynLmi(
+            lmi_strictness=1e-7,
+            solver_params=dict(
+                solver="MOSEK",
+                eps=1e-8,
+            ),
+        ),
+        structured_singular_value=dkpy.SsvLmiBisection(
+            bisection_atol=1e-5,
+            bisection_rtol=1e-5,
+            max_iterations=1000,
+            lmi_strictness=1e-7,
+            solver_params=dict(
+                solver="MOSEK",
+                eps=1e-9,
+            ),
+        ),
+        transfer_function_fit=dkpy.TfFitSlicot(),
+        fit_order_callback=callback,
+    )
+
+    omega = np.logspace(-3, 3, 61)
+    block_structure = np.array([[1, 1], [1, 1], [2, 2]])
+    K, N, mu, d_scale_fit_info, info = dk_iter.synthesize(
+        P,
+        n_y,
+        n_u,
+        omega,
+        block_structure,
+    )
+
+    print(mu)
+
+    fig, ax = plt.subplots()
+    for i, ds in enumerate(d_scale_fit_info):
+        ds.plot_mu(ax=ax, plot_kw=dict(label=f"iter{i}"))
+
+    ax = None
+    for i, ds in enumerate(d_scale_fit_info):
+        _, ax = ds.plot_D(ax=ax, plot_kw=dict(label=f"iter{i}"))
+
+    plt.show()
+
+
+if __name__ == "__main__":
+    example_dk_iter_fixed_order()

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "hatchling.build"
 
 [project]
 name = "dkpy"
-version = "0.1.4"
+version = "0.1.5"
 dependencies = [
     "numpy>=1.21.0",
     "scipy>=1.7.0",

requirements.txt

@@ -1,14 +1,22 @@
+# Dependencies
 numpy
 scipy
 control
+slycot
 cvxpy
 joblib
-slycot
 matplotlib
+
+# Extra solvers
+mosek>=9.2.49
+
+# Development
 pytest
 pytest-regressions
 sphinx
 sphinx-rtd-theme
 ruff==0.8.0
 pre-commit
-mosek>=9.2.49
+
+# Examples
+PyQt6

src/dkpy/dk_iteration.py

@@ -5,10 +5,11 @@
     "DkIteration",
     "DkIterFixedOrder",
     "DkIterListOrder",
+    "DkIterOrderCallback",
 ]
 
 import abc
-from typing import Any, Dict, List, Tuple, Union, Optional
+from typing import Any, Callable, Dict, List, Tuple, Union, Optional
 
 import control
 import numpy as np
@@ -534,6 +535,120 @@ def synthesize(
         return (K, N, np.max(mu_omega), d_scale_fit_info, info)
 
 
+class DkIterOrderCallback(DkIteration):
+    """D-K iteration with a callback to select fit orders."""
+
+    def __init__(
+        self,
+        controller_synthesis: controller_synthesis.ControllerSynthesis,
+        structured_singular_value: structured_singular_value.StructuredSingularValue,
+        transfer_function_fit: fit_transfer_functions.TransferFunctionFit,
+        fit_order_callback: Optional[Callable] = None,
+    ):
+        """Instantiate :class:`DkIterListOrder`.
+
+        Parameters
+        ----------
+        controller_synthesis : dkpy.ControllerSynthesis
+            A controller synthesis object.
+        structured_singular_value : dkpy.StructuredSingularValue
+            A structured singular value computation object.
+        transfer_function_fit : dkpy.TransferFunctionFit
+            A transfer function fit object.
+        fit_order_callback : Optional[Callable]
+            Callback to select D-scale fit orders.
+        """
+        self.controller_synthesis = controller_synthesis
+        self.structured_singular_value = structured_singular_value
+        self.transfer_function_fit = transfer_function_fit
+        self.fit_order_callback = fit_order_callback
+
+    def synthesize(
+        self,
+        P: control.StateSpace,
+        n_y: int,
+        n_u: int,
+        omega: np.ndarray,
+        block_structure: np.ndarray,
+    ) -> Tuple[
+        control.StateSpace,
+        control.StateSpace,
+        float,
+        List[DScaleFitInfo],
+        Dict[str, Any],
+    ]:
+        # Solution information
+        info = {}
+        d_scale_fit_info = []
+        # Set up initial D-scales
+        D = _get_initial_d_scales(block_structure)
+        D_inv = _get_initial_d_scales(block_structure)
+        D_aug, D_aug_inv = _augment_d_scales(D, D_inv, n_y=n_y, n_u=n_u)
+        iteration = 0
+        done = False
+        while not done:
+            # Synthesize controller
+            K, _, gamma, info = self.controller_synthesis.synthesize(
+                D_aug * P * D_aug_inv,
+                n_y,
+                n_u,
+            )
+            N = P.lft(K)
+            # Compute structured singular values on grid
+            N_omega = N(1j * omega)
+            mu_omega, D_omega, info = self.structured_singular_value.compute_ssv(
+                N_omega,
+                block_structure=block_structure,
+            )
+            fit_order, done = self.fit_order_callback(
+                self,
+                iteration,
+                omega,
+                mu_omega,
+                D_omega,
+                P,
+                K,
+                block_structure,
+            )
+            # Fit transfer functions to gridded D-scales
+            D_fit, D_fit_inv = self.transfer_function_fit.fit(
+                omega,
+                D_omega,
+                order=fit_order,
+                block_structure=block_structure,
+            )
+            # Add D-scale fit info
+            d_scale_fit_info.append(
+                DScaleFitInfo.create_from_fit(
+                    omega,
+                    mu_omega,
+                    D_omega,
+                    P,
+                    K,
+                    D_fit,
+                    D_fit_inv,
+                    block_structure,
+                )
+            )
+            # Augment D-scales with identity transfer functions
+            D_aug, D_aug_inv = _augment_d_scales(
+                D_fit,
+                D_fit_inv,
+                n_y=n_y,
+                n_u=n_u,
+            )
+            # Increment iteration
+            iteration += 1
+        # Synthesize controller one last time
+        K, _, gamma, info = self.controller_synthesis.synthesize(
+            D_aug * P * D_aug_inv,
+            n_y,
+            n_u,
+        )
+        N = P.lft(K)
+        return (K, N, np.max(mu_omega), d_scale_fit_info, info)
+
+
 def _get_initial_d_scales(block_structure: np.ndarray) -> control.StateSpace:
     """Generate initial identity D-scales based on block structure.
 

tests/test_fit_transfer_functions.py

@@ -163,7 +163,7 @@ class TestTfFitSlicot:
                         ],
                     ],
                 ),
-                np.diag([1, 1]),
+                np.diag([1, 0]),
                 np.array([[1, 1], [1, 1]]),
                 1e-2,
             ),