diff --git a/CITATION.cff b/CITATION.cff index 41f2410..2fa1b17 100644 --- a/CITATION.cff +++ b/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" diff --git a/README.rst b/README.rst index 583d1cd..a58f0ea 100644 --- a/README.rst +++ b/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}, } diff --git a/doc/conf.py b/doc/conf.py index 3a2340b..0eb9fbd 100644 --- a/doc/conf.py +++ b/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 diff --git a/examples/3_example_dk_iter_order_callback.py b/examples/3_example_dk_iter_order_callback.py new file mode 100644 index 0000000..8a5af70 --- /dev/null +++ b/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() diff --git a/pyproject.toml b/pyproject.toml index aae1136..e9c51b1 100644 --- a/pyproject.toml +++ b/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", diff --git a/requirements.txt b/requirements.txt index d7d679e..932fbf3 100644 --- a/requirements.txt +++ b/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 diff --git a/src/dkpy/dk_iteration.py b/src/dkpy/dk_iteration.py index 85eb2da..86dd333 100644 --- a/src/dkpy/dk_iteration.py +++ b/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. diff --git a/tests/test_fit_transfer_functions.py b/tests/test_fit_transfer_functions.py index 0b53e5e..45cc952 100644 --- a/tests/test_fit_transfer_functions.py +++ b/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, ),