Merge pull request #19 from decargroup/feature/10-automatically-selec…

…t-order-by-minimizing-peak-mu-error

Automatically select order by minimizing peak mu error

examples/2_example_dk_iter_list_order.py

@@ -7,7 +7,7 @@
 import dkpy
 
 
-def example_dk_iter_fixed_order():
+def example_dk_iter_list_order():
     """D-K iteration with fixed number of iterations and fit order."""
     # Plant
     G0 = np.array(
@@ -117,4 +117,4 @@ def example_dk_iter_fixed_order():
 
 
 if __name__ == "__main__":
-    example_dk_iter_fixed_order()
+    example_dk_iter_list_order()

examples/3_example_dk_iter_order_callback.py → examples/3_example_dk_iter_custom.py

@@ -53,7 +53,7 @@ def _get_fit_order(
             return int(selected_order_str)
 
 
-def example_dk_iter_fixed_order():
+def example_dk_iter_custom():
     """D-K iteration with fixed number of iterations and fit order."""
     # Plant
     G0 = np.array(
@@ -157,4 +157,4 @@ def example_dk_iter_fixed_order():
 
 
 if __name__ == "__main__":
-    example_dk_iter_fixed_order()
+    example_dk_iter_custom()

examples/4_example_dk_iter_auto_order.py

@@ -0,0 +1,118 @@
+"""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_auto_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
+
+    dk_iter = dkpy.DkIterAutoOrder(
+        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(),
+        max_mu=1,
+        max_mu_fit_error=1e-2,
+        max_iterations=3,
+        max_fit_order=4,
+    )
+
+    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):
+        dkpy.plot_mu(ds, ax=ax, plot_kw=dict(label=f"iter{i}"))
+
+    ax = None
+    for i, ds in enumerate(d_scale_fit_info):
+        _, ax = dkpy.plot_D(ds, ax=ax, plot_kw=dict(label=f"iter{i}"))
+
+    plt.show()
+
+
+if __name__ == "__main__":
+    example_dk_iter_auto_order()

src/dkpy/dk_iteration.py

@@ -5,6 +5,7 @@
     "DkIteration",
     "DkIterFixedOrder",
     "DkIterListOrder",
+    "DkIterAutoOrder",
 ]
 
 import abc
@@ -397,6 +398,97 @@ def _get_fit_order(
             return None
 
 
+class DkIterAutoOrder(DkIteration):
+    """D-K iteration with a fixed list of fit orders."""
+
+    def __init__(
+        self,
+        controller_synthesis: controller_synthesis.ControllerSynthesis,
+        structured_singular_value: structured_singular_value.StructuredSingularValue,
+        transfer_function_fit: fit_transfer_functions.TransferFunctionFit,
+        max_mu: float = 1,
+        max_mu_fit_error: float = 1e-2,
+        max_iterations: Optional[int] = None,
+        max_fit_order: Optional[int] = 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.
+        max_mu : float
+            Maximum acceptable structured singular value.
+        max_mu_fit_error : float
+            Maximum relative fit error for structured singular value.
+        max_iterations : Optional[int]
+            Maximum number of iterations. If ``None``, there is no upper limit.
+        max_fit_order : Optional[int]
+            Maximum fit order. If ``None``, there is no upper limit.
+        """
+        self.controller_synthesis = controller_synthesis
+        self.structured_singular_value = structured_singular_value
+        self.transfer_function_fit = transfer_function_fit
+        self.max_mu = max_mu
+        self.max_mu_fit_error = max_mu_fit_error
+        self.max_iterations = max_iterations
+        self.max_fit_order = max_fit_order
+
+    def _get_fit_order(
+        self,
+        iteration: int,
+        omega: np.ndarray,
+        mu_omega: np.ndarray,
+        D_omega: np.ndarray,
+        P: control.StateSpace,
+        K: control.StateSpace,
+        block_structure: np.ndarray,
+    ) -> Optional[Union[int, np.ndarray]]:
+        print(f"{iteration}")  # TODO
+        # Check termination conditions
+        if (self.max_iterations is not None) and (iteration >= self.max_iterations):
+            # TODO LOG
+            return None
+        if np.max(mu_omega) < self.max_mu:
+            # TODO LOG
+            return None
+        # Determine fit order
+        fit_order = 0
+        relative_errors = []
+        while True:
+            D_fit, D_fit_inv = self.transfer_function_fit.fit(
+                omega,
+                D_omega,
+                order=fit_order,
+                block_structure=block_structure,
+            )
+            d_scale_fit_info = DScaleFitInfo.create_from_fit(
+                omega,
+                mu_omega,
+                D_omega,
+                P,
+                K,
+                D_fit,
+                D_fit_inv,
+                block_structure,
+            )
+            error = np.abs(mu_omega - d_scale_fit_info.mu_fit_omega)
+            relative_error = np.max(error / np.max(np.abs(mu_omega)))
+            print(f"    {fit_order}: {relative_error}")  # TODO
+            relative_errors.append(relative_error)
+            if (self.max_fit_order is not None) and (fit_order >= self.max_fit_order):
+                # TODO LOG
+                return int(np.argmin(relative_errors))
+            if relative_error >= self.max_mu_fit_error:
+                fit_order += 1
+            else:
+                return fit_order
+
+
 def _augment_d_scales(
     D: Union[control.TransferFunction, control.StateSpace],
     D_inv: Union[control.TransferFunction, control.StateSpace],