diff --git a/CITATION.cff b/CITATION.cff index c7ad650..71a5d66 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -14,5 +14,5 @@ authors: orcid: "https://orcid.org/0000-0002-1987-9268" affiliation: "McGill University" title: "decargroup/dkpy" -version: v0.1.10 +version: v0.1.11 url: "https://github.com/decargroup/dkpy" diff --git a/README.rst b/README.rst index bb19896..2ec16e4 100644 --- a/README.rst +++ b/README.rst @@ -141,6 +141,6 @@ If you use this software in your research, please cite it as below or see url={https://github.com/decargroup/dkpy}, publisher={Zenodo}, author={Steven Dahdah and Timothy Everett Adams and James Richard Forbes}, - version = {{v0.1.10}}, + version = {{v0.1.11}}, year={2025}, } diff --git a/doc/_static/example_1/1_plot_D.png b/doc/_static/example_1/1_plot_D.png new file mode 100644 index 0000000..6bfced4 Binary files /dev/null and b/doc/_static/example_1/1_plot_D.png differ diff --git a/doc/_static/example_1/1_plot_mu.png b/doc/_static/example_1/1_plot_mu.png new file mode 100644 index 0000000..3b6747d Binary files /dev/null and b/doc/_static/example_1/1_plot_mu.png differ diff --git a/doc/_static/example_2/2_plot_D.png b/doc/_static/example_2/2_plot_D.png new file mode 100644 index 0000000..6bfced4 Binary files /dev/null and b/doc/_static/example_2/2_plot_D.png differ diff --git a/doc/_static/example_2/2_plot_mu.png b/doc/_static/example_2/2_plot_mu.png new file mode 100644 index 0000000..3b6747d Binary files /dev/null and b/doc/_static/example_2/2_plot_mu.png differ diff --git a/doc/_static/example_3/3_plot_D.png b/doc/_static/example_3/3_plot_D.png new file mode 100644 index 0000000..126616c Binary files /dev/null and b/doc/_static/example_3/3_plot_D.png differ diff --git a/doc/_static/example_3/3_plot_mu.png b/doc/_static/example_3/3_plot_mu.png new file mode 100644 index 0000000..5689b2c Binary files /dev/null and b/doc/_static/example_3/3_plot_mu.png differ diff --git a/doc/_static/example_4/4_interactive_1.png b/doc/_static/example_4/4_interactive_1.png new file mode 100644 index 0000000..81ffd72 Binary files /dev/null and b/doc/_static/example_4/4_interactive_1.png differ diff --git a/doc/_static/example_4/4_interactive_2.png b/doc/_static/example_4/4_interactive_2.png new file mode 100644 index 0000000..81a6534 Binary files /dev/null and b/doc/_static/example_4/4_interactive_2.png differ diff --git a/doc/_static/example_4/4_interactive_3.png b/doc/_static/example_4/4_interactive_3.png new file mode 100644 index 0000000..f6f1116 Binary files /dev/null and b/doc/_static/example_4/4_interactive_3.png differ diff --git a/doc/_static/example_4/4_interactive_4.png b/doc/_static/example_4/4_interactive_4.png new file mode 100644 index 0000000..0f198f3 Binary files /dev/null and b/doc/_static/example_4/4_interactive_4.png differ diff --git a/doc/_static/example_4/4_plot_D.png b/doc/_static/example_4/4_plot_D.png new file mode 100644 index 0000000..126616c Binary files /dev/null and b/doc/_static/example_4/4_plot_D.png differ diff --git a/doc/_static/example_4/4_plot_mu.png b/doc/_static/example_4/4_plot_mu.png new file mode 100644 index 0000000..5689b2c Binary files /dev/null and b/doc/_static/example_4/4_plot_mu.png differ diff --git a/doc/_static/example_5/5_plot_D.png b/doc/_static/example_5/5_plot_D.png new file mode 100644 index 0000000..126616c Binary files /dev/null and b/doc/_static/example_5/5_plot_D.png differ diff --git a/doc/_static/example_5/5_plot_mu.png b/doc/_static/example_5/5_plot_mu.png new file mode 100644 index 0000000..5689b2c Binary files /dev/null and b/doc/_static/example_5/5_plot_mu.png differ diff --git a/doc/_static/example_6/6_plot_inputs.png b/doc/_static/example_6/6_plot_inputs.png new file mode 100644 index 0000000..55b02f8 Binary files /dev/null and b/doc/_static/example_6/6_plot_inputs.png differ diff --git a/doc/_static/example_6/6_plot_mu.png b/doc/_static/example_6/6_plot_mu.png new file mode 100644 index 0000000..c3cb617 Binary files /dev/null and b/doc/_static/example_6/6_plot_mu.png differ diff --git a/doc/_static/example_6/6_plot_states.png b/doc/_static/example_6/6_plot_states.png new file mode 100644 index 0000000..e4548dd Binary files /dev/null and b/doc/_static/example_6/6_plot_states.png differ diff --git a/doc/_static/interactive_1.png b/doc/_static/interactive_1.png deleted file mode 100644 index c86e585..0000000 Binary files a/doc/_static/interactive_1.png and /dev/null differ diff --git a/doc/_static/interactive_2.png b/doc/_static/interactive_2.png deleted file mode 100644 index 8b7d190..0000000 Binary files a/doc/_static/interactive_2.png and /dev/null differ diff --git a/doc/_static/interactive_3.png b/doc/_static/interactive_3.png deleted file mode 100644 index b687129..0000000 Binary files a/doc/_static/interactive_3.png and /dev/null differ diff --git a/doc/_static/interactive_4.png b/doc/_static/interactive_4.png deleted file mode 100644 index 4a75533..0000000 Binary files a/doc/_static/interactive_4.png and /dev/null differ diff --git a/doc/_static/plot_D.png b/doc/_static/plot_D.png deleted file mode 100644 index 17fb944..0000000 Binary files a/doc/_static/plot_D.png and /dev/null differ diff --git a/doc/_static/plot_mu.png b/doc/_static/plot_mu.png deleted file mode 100644 index 2d2190f..0000000 Binary files a/doc/_static/plot_mu.png and /dev/null differ diff --git a/doc/conf.py b/doc/conf.py index 18484e5..bfc9073 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -9,8 +9,8 @@ project = "dkpy" 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" +version = "0.1.11" +release = "0.1.11" # -- 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 9371a18..55fbadb 100644 --- a/doc/dkpy.rst +++ b/doc/dkpy.rst @@ -18,12 +18,15 @@ with them. The abstract classes are provided below. 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. diff --git a/doc/examples.rst b/doc/examples.rst index 68d128f..7ecd069 100644 --- a/doc/examples.rst +++ b/doc/examples.rst @@ -20,10 +20,10 @@ to 4. Output:: - mu=1.0132789611816406 + mu=1.0368156433105469 -.. image:: _static/plot_mu.png -.. image:: _static/plot_D.png +.. image:: _static/example_1/1_plot_mu.png +.. image:: _static/example_1/1_plot_D.png D-K iteration with list of fit orders ------------------------------------- @@ -37,10 +37,10 @@ integer. Output:: - mu=1.0132789611816406 + mu=1.0368156433105469 -.. image:: _static/plot_mu.png -.. image:: _static/plot_D.png +.. image:: _static/example_2/2_plot_mu.png +.. image:: _static/example_2/2_plot_D.png D-K iteration with automatically selected fit orders ---------------------------------------------------- @@ -53,34 +53,34 @@ with the lowest relative error is selected. Output:: - INFO:DkIterAutoOrder:Iteration: 0, mu: 1.1792325973510742 - INFO:DkIterAutoOrder:Order 0 relative error: 0.5122457769147215 - INFO:DkIterAutoOrder:Order 1 relative error: 0.34431183690576633 - INFO:DkIterAutoOrder:Order 2 relative error: 0.8970100659296376 - INFO:DkIterAutoOrder:Order 3 relative error: 0.030844892155775263 - INFO:DkIterAutoOrder:Order 4 relative error: 0.015896380940858944 + INFO:DkIterAutoOrder:Iteration: 0, mu: 1.1792182922363281 + INFO:DkIterAutoOrder:Order 0 relative error: 0.5066119732511977 + INFO:DkIterAutoOrder:Order 1 relative error: 0.3356881131800033 + INFO:DkIterAutoOrder:Order 2 relative error: 0.643462672042705 + INFO:DkIterAutoOrder:Order 3 relative error: 0.031649076289004166 + INFO:DkIterAutoOrder:Order 4 relative error: 0.012866875729555234 INFO:DkIterAutoOrder:Reached max fit order, selecting order 4 - INFO:DkIterAutoOrder:Iteration: 1, mu: 1.0256481170654297 - INFO:DkIterAutoOrder:Order 0 relative error: 9.350642790785383 - INFO:DkIterAutoOrder:Order 1 relative error: 1.327796980729705 - INFO:DkIterAutoOrder:Order 2 relative error: 7.58063969442474 - INFO:DkIterAutoOrder:Order 3 relative error: 0.13472058625314948 - INFO:DkIterAutoOrder:Order 4 relative error: 0.05262584627773765 + INFO:DkIterAutoOrder:Iteration: 1, mu: 1.0274028778076172 + INFO:DkIterAutoOrder:Order 0 relative error: 9.780743577129419 + INFO:DkIterAutoOrder:Order 1 relative error: 1.4151255816143613 + INFO:DkIterAutoOrder:Order 2 relative error: 8.507402105789147 + INFO:DkIterAutoOrder:Order 3 relative error: 0.15216618301078452 + INFO:DkIterAutoOrder:Order 4 relative error: 0.05185270360448734 INFO:DkIterAutoOrder:Reached max fit order, selecting order 4 - INFO:DkIterAutoOrder:Iteration: 2, mu: 1.0201168060302734 - INFO:DkIterAutoOrder:Order 0 relative error: 28.96943897648902 - INFO:DkIterAutoOrder:Order 1 relative error: 3.9918344794684804 - INFO:DkIterAutoOrder:Order 2 relative error: 7.000978128580445 - INFO:DkIterAutoOrder:Order 3 relative error: 0.22159755671770862 - INFO:DkIterAutoOrder:Order 4 relative error: 0.06377761160336083 + INFO:DkIterAutoOrder:Iteration: 2, mu: 1.0203123092651367 + INFO:DkIterAutoOrder:Order 0 relative error: 30.002715934629506 + INFO:DkIterAutoOrder:Order 1 relative error: 4.172954343879354 + INFO:DkIterAutoOrder:Order 2 relative error: 9.026259323212026 + INFO:DkIterAutoOrder:Order 3 relative error: 0.24098884948353583 + INFO:DkIterAutoOrder:Order 4 relative error: 0.04541136749668787 INFO:DkIterAutoOrder:Reached max fit order, selecting order 4 - INFO:DkIterAutoOrder:Iteration: 3, mu: 1.0132789611816406 + INFO:DkIterAutoOrder:Iteration: 3, mu: 1.0144329071044922 INFO:DkIterAutoOrder:Iteration terminated: reached maximum number of iterations INFO:DkIterAutoOrder:Iteration complete - mu=1.0132789611816406 + mu=1.0144329071044922 -.. image:: _static/plot_mu.png -.. image:: _static/plot_D.png +.. image:: _static/example_3/3_plot_mu.png +.. image:: _static/example_3/3_plot_D.png D-K iteration with interactively selected fit orders ---------------------------------------------------- @@ -92,38 +92,38 @@ singular value plots at each iteration. .. literalinclude:: ../examples/4_example_dk_iter_interactive.py :language: python -.. image:: _static/interactive_1.png +.. image:: _static/example_4/4_interactive_1.png Prompt:: Close plot to continue... Select order ( to end iteration): 4 -.. image:: _static/interactive_2.png +.. image:: _static/example_4/4_interactive_2.png Prompt:: Close plot to continue... Select order ( to end iteration): 4 -.. image:: _static/interactive_3.png +.. image:: _static/example_4/4_interactive_3.png Prompt:: Close plot to continue... Select order ( to end iteration): 4 -.. image:: _static/interactive_4.png +.. image:: _static/example_4/4_interactive_4.png Output:: Close plot to continue... Select order ( to end iteration): Iteration ended. - mu=1.0132789611816406 + mu=1.0144329071044922 -.. image:: _static/plot_mu.png -.. image:: _static/plot_D.png +.. image:: _static/example_4/4_plot_mu.png +.. image:: _static/example_4/4_plot_D.png D-K iteration with a custom fit order selection method ------------------------------------------------------ @@ -136,11 +136,44 @@ after 3 iterations of 4th order fits. Output:: - Iteration 0 with mu of 1.1792325973510742 - Iteration 1 with mu of 1.0256481170654297 - Iteration 2 with mu of 1.0201168060302734 - Iteration 3 with mu of 1.0132789611816406 - mu=1.0132789611816406 + Iteration 0 with mu of 1.1792182922363281 + Iteration 1 with mu of 1.0274028778076172 + Iteration 2 with mu of 1.0203123092651367 + Iteration 3 with mu of 1.0144329071044922 + mu=1.0144329071044922 -.. image:: _static/plot_mu.png -.. image:: _static/plot_D.png +.. image:: _static/example_5/5_plot_mu.png +.. image:: _static/example_5/5_plot_D.png + + +D-K iteration for non-square perturbation and simulation of perturbed systems +----------------------------------------------------------------------------- + +In this example, the perturbation `Δ` is non-square. In other words, the inputs +and outputs of the perturbation are not identical. The orders are specified in +a list. Once the robust controller is synthesized, it is tested on a set of +perturbed systems by generating a set of perturbations `Δ` that have H-infinity +norm less than or equal to 1. + +In this example, a controller is designed for the linearized lateral dynamics +of an aircraft from an example given in Section 14.1 of [M04]_. + +.. literalinclude:: ../examples/6_dk_iteration_non_square_perturbation.py + :language: python + +Output:: + + mu: 0.9824275970458984 + +.. image:: _static/example_6/6_plot_mu.png + +After the controller is synthesized, a set of off-nominal models are generated +by interconnecting various admissible perturbations `Δ` with the nominal model. +Then, the closed-loop systems are formed for each perturbed system. The time +domain response of the systems can be obtained for a step response in the roll +angle reference signal, no external disturbances, and moderate sensor noise. In +particular, the response of the system states and actuator inputs are shown. + +.. image:: _static/example_6/6_plot_states.png + +.. image:: _static/example_6/6_plot_inputs.png diff --git a/doc/references.rst b/doc/references.rst index a6e2ec8..e43f080 100644 --- a/doc/references.rst +++ b/doc/references.rst @@ -10,3 +10,7 @@ References .. [SP06] S. Skogestad and I. Postlethwaite, Multivariable Feedback Control: Analysis and Design. West Sussex, England: Wiley, 2006. ISBN 978-0-47-001168-3. +.. [M04] U. Mackenroth, Robust Control Systems: Theory and Case Studies. + Berlin, Heidelberg: Springer, 2004. ISBN 978-3-642-05891-2. + + diff --git a/examples/1_example_dk_iter_fixed_order.py b/examples/1_example_dk_iter_fixed_order.py index 318dfa4..dc1fb3d 100644 --- a/examples/1_example_dk_iter_fixed_order.py +++ b/examples/1_example_dk_iter_fixed_order.py @@ -12,24 +12,22 @@ def example_dk_iter_fixed_order(): dk_iter = dkpy.DkIterFixedOrder( controller_synthesis=dkpy.HinfSynSlicot(), - structured_singular_value=dkpy.SsvLmiBisection(), + structured_singular_value=dkpy.SsvLmiBisection(n_jobs=1), d_scale_fit=dkpy.DScaleFitSlicot(), n_iterations=3, fit_order=4, ) omega = np.logspace(-3, 3, 61) - # Alternative MATLAB descr. + # Alternative MATLAB block structure description # 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), - ] - ) + uncertainty_structure = [ + 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"], @@ -46,7 +44,7 @@ def example_dk_iter_fixed_order(): ax = None for i, ds in enumerate(iter_results): - _, ax = dkpy.plot_D(ds, ax=ax, plot_kw=dict(label=f"iter{i}")) + fig, ax = dkpy.plot_D(ds, ax=ax, plot_kw=dict(label=f"iter{i}")) plt.show() diff --git a/examples/2_example_dk_iter_list_order.py b/examples/2_example_dk_iter_list_order.py index 96d8484..05dda94 100644 --- a/examples/2_example_dk_iter_list_order.py +++ b/examples/2_example_dk_iter_list_order.py @@ -16,24 +16,22 @@ def example_dk_iter_list_order(): d_scale_fit=dkpy.DScaleFitSlicot(), # fit_orders=[4, 4, 4], # <- an alternative fit_orders=[ - np.diag([4, 4, 0, 0]), - np.diag([4, 4, 0, 0]), - np.diag([4, 4, 0, 0]), + [4, 4, 0], + [4, 4, 0], + [4, 4, 0], ], ) omega = np.logspace(-3, 3, 61) - # Alternative MATLAB descr. + # Alternative MATLAB block structure description # 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), - ] - ) + uncertainty_structure = [ + 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"], @@ -50,7 +48,7 @@ def example_dk_iter_list_order(): ax = None for i, ds in enumerate(iter_results): - _, ax = dkpy.plot_D(ds, ax=ax, plot_kw=dict(label=f"iter{i}")) + fig, ax = dkpy.plot_D(ds, ax=ax, plot_kw=dict(label=f"iter{i}")) plt.show() diff --git a/examples/3_example_dk_iter_auto_order.py b/examples/3_example_dk_iter_auto_order.py index 9ba2b17..3a88c5a 100644 --- a/examples/3_example_dk_iter_auto_order.py +++ b/examples/3_example_dk_iter_auto_order.py @@ -35,6 +35,7 @@ def example_dk_iter_auto_order(): solver="MOSEK", eps=1e-9, ), + n_jobs=1, ), d_scale_fit=dkpy.DScaleFitSlicot(), max_mu=1, @@ -44,17 +45,15 @@ def example_dk_iter_auto_order(): ) omega = np.logspace(-3, 3, 61) - # Alternative MATLAB descr. + # Alternative MATLAB block structure description # 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), - ] - ) + uncertainty_structure = [ + 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"], @@ -71,8 +70,7 @@ def example_dk_iter_auto_order(): ax = None for i, ds in enumerate(iter_results): - _, ax = dkpy.plot_D(ds, ax=ax, plot_kw=dict(label=f"iter{i}")) - + fig, ax = dkpy.plot_D(ds, ax=ax, plot_kw=dict(label=f"iter{i}")) plt.show() diff --git a/examples/4_example_dk_iter_interactive.py b/examples/4_example_dk_iter_interactive.py index 368e9c3..ee5368a 100644 --- a/examples/4_example_dk_iter_interactive.py +++ b/examples/4_example_dk_iter_interactive.py @@ -37,17 +37,15 @@ def example_dk_iter_interactive(): ) omega = np.logspace(-3, 3, 61) - # Alternative MATLAB descr. + # Alternative MATLAB block structure description # 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), - ] - ) + uncertainty_structure = [ + 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"], @@ -64,7 +62,7 @@ def example_dk_iter_interactive(): ax = None for i, ds in enumerate(iter_results): - _, ax = dkpy.plot_D(ds, ax=ax, plot_kw=dict(label=f"iter{i}")) + fig, ax = dkpy.plot_D(ds, ax=ax, plot_kw=dict(label=f"iter{i}")) plt.show() diff --git a/examples/5_example_dk_iteration_custom.py b/examples/5_example_dk_iteration_custom.py index 8f01287..19c0c1a 100644 --- a/examples/5_example_dk_iteration_custom.py +++ b/examples/5_example_dk_iteration_custom.py @@ -31,10 +31,11 @@ def _get_fit_order( iteration, omega, mu_omega, - D_omega, + D_l_omega, + D_r_omega, P, K, - uncertainty_structure, + block_structure, ): print(f"Iteration {self.my_iter_count} with mu of {np.max(mu_omega)}") if self.my_iter_count < 3: @@ -65,22 +66,21 @@ def example_dk_iter_custom(): solver="MOSEK", eps=1e-9, ), + n_jobs=1, ), d_scale_fit=dkpy.DScaleFitSlicot(), ) omega = np.logspace(-3, 3, 61) - # Alternative MATLAB descr. + # Alternative MATLAB block structure description # 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), - ] - ) + uncertainty_structure = [ + 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"], @@ -97,7 +97,7 @@ def example_dk_iter_custom(): ax = None for i, ds in enumerate(iter_results): - _, ax = dkpy.plot_D(ds, ax=ax, plot_kw=dict(label=f"iter{i}")) + fig, ax = dkpy.plot_D(ds, ax=ax, plot_kw=dict(label=f"iter{i}")) plt.show() diff --git a/examples/6_dk_iteration_non_square_perturbation.py b/examples/6_dk_iteration_non_square_perturbation.py new file mode 100644 index 0000000..e5877d2 --- /dev/null +++ b/examples/6_dk_iteration_non_square_perturbation.py @@ -0,0 +1,361 @@ +"""D-K iteration with nonsquare perturbations + +D-K iteration with a fixed number of iterations and fit order and closed-loop +simluation of perturbed models for the linearized lateral dynamics of an +aircraft. + +The quantities of interest in the aircraft dynamics are shown below. + +Airframe States +--------------- +phi: Roll angle +beta: Sideslip angle +p: Roll rate +r: Yaw rate + +Control Inputs +-------------- +zeta: Rudder angle input +xi: Aileron angle input + +Disturbances +------------ +beta_w: Sideslip angle wind disturbance +p_w: Roll rate wind disturbance + +Tracking Reference +------------------ +phi_ref: Roll angle reference +""" + +import numpy as np +import control +from matplotlib import pyplot as plt + +import dkpy + + +def example_dk_iter_list_order_aircraft(): + # Example parameters + eg = dkpy.utilities.example_mackenroth2004_p430() + plant_gen = eg["P"] + + n_u = eg["n_u"] + n_y = eg["n_y"] + n_u_delta = eg["n_u_delta"] + n_y_delta = eg["n_y_delta"] + n_w = eg["n_w"] + n_z = eg["n_z"] + + # Frequency + freq_min = 0.01 + freq_max = 100 + num_freq = 100 + freq = np.logspace(np.log10(freq_min), np.log10(freq_max), num_freq) + omega = 2 * np.pi * freq + + # DK-iteration controller synthesis + dk_iter = dkpy.DkIterListOrder( + controller_synthesis=dkpy.HinfSynLmi( + lmi_strictness=5e-7, # Have to play with tolerances. + solver_params=dict( + solver="MOSEK", + eps=5e-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, + ), + ), + d_scale_fit=dkpy.DScaleFitSlicot(), + fit_orders=[4, 4], + ) + # Alternative MATLAB block structure description + # block_structure = np.array([[n_u_delta, n_y_delta], [n_w, n_z]]) + block_structure = [ + dkpy.ComplexFullBlock(n_y_delta, n_u_delta), + dkpy.ComplexFullBlock(n_w, n_z), + ] + controller, N, mu, iter_results, info = dk_iter.synthesize( + plant_gen, + n_y, + n_u, + omega, + block_structure, + ) + controller.set_inputs(["phi_ref", "phi_meas", "beta_meas", "p_meas", "r_meas"]) + controller.set_outputs(["zeta_c", "xi_c"]) + print("mu:", mu) + + # Plot - DK-iteration results + fig, ax = plt.subplots() + for i, ds in enumerate(iter_results): + dkpy.plot_mu(ds, ax=ax, plot_kw=dict(label=f"iter{i}")) + + # Generate off-nominal models with different perturbations + delta_block_list = [] + # Nominal system uncertainty (no perturbation) + delta_block = control.StateSpace([], [], [], [[0, 0], [0, 0]]) + delta_block.set_inputs(delta_block.ninputs, "y_del") + delta_block.set_outputs(delta_block.noutputs, "u_del") + delta_block_list.append(delta_block) + # Constant gain perturbation + num_offnom_gain = 20 + for gain_unc in np.linspace(-1, 1, num_offnom_gain): + delta = control.StateSpace([], [], [], [gain_unc]) + delta_block = control.append(delta, delta) + delta_block.set_inputs(delta_block.ninputs, "y_del") + delta_block.set_outputs(delta_block.noutputs, "u_del") + delta_block_list.append(delta_block) + # Phase perturbation + a_min = 0.01 + a_max = 10 + num_offnom_phase = 20 + for a in np.linspace(a_min, a_max, num_offnom_phase): + delta_block = control.append( + control.tf([1 / a, -1], [1 / a, 1]), + control.tf([a, -1], [a, 1]), + name="delta_block", + ) + delta_block.set_inputs(delta_block.ninputs, "y_del") + delta_block.set_outputs(delta_block.noutputs, "u_del") + delta_block_list.append(delta_block) + + # Closed-loop simulation system + input_id_sim_list = [ + "phi_ref", + "beta_w", + "p_w", + "n[0]", + "n[1]", + "n[2]", + "n[3]", + ] + output_id_sim_list = [ + "phi", + "beta", + "p", + "r", + "zeta_c", + "xi_c", + "zeta", + "xi", + "rate_xi", + "rate_zeta", + ] + + # Uncertain closed-loop system + airframe = eg["airframe"] + actuator = eg["actuator"] + weight_unc = eg["weight_unc"] + sum_noise = eg["sum_noise"] + sum_uncertainty = eg["sum_uncertainty"] + closed_loop_sim_list = [] + for delta_block in delta_block_list: + closed_loop_sim = control.interconnect( + syslist=[ + airframe, + actuator, + controller, + weight_unc, + delta_block, + sum_noise, + sum_uncertainty, + ], + inplist=input_id_sim_list, + outlist=output_id_sim_list, + name="closed_loop_sim", + ) + closed_loop_sim.set_inputs(input_id_sim_list) + closed_loop_sim.set_outputs(output_id_sim_list) + closed_loop_sim_list.append(closed_loop_sim) + + # Time + time_min = 0 + time_max = 7.5 + dt = 0.01 + num_time = round((time_max - time_min) / dt) + time = np.arange(time_min, time_max, dt) + # Reference signal + phi_ref = 6 * np.ones_like(time) + # Disturbance signals + beta_w = np.zeros_like(time) + p_w = np.zeros_like(time) + # Sensor noise signals (assume sensors have identical noise characteristics) + mean_noise = 0 + std_noise = 0.005 + noise_phi = np.random.normal(mean_noise, std_noise, num_time) + noise_beta = np.random.normal(mean_noise, std_noise, num_time) + noise_p = np.random.normal(mean_noise, std_noise, num_time) + noise_r = np.random.normal(mean_noise, std_noise, num_time) + # Exogenous input signal + inputs_sim = np.array( + [phi_ref, beta_w, p_w, noise_phi, noise_beta, noise_p, noise_r] + ) + + # Simulation response + trd_sim_list = control.forced_response( + sysdata=closed_loop_sim_list, T=time, U=inputs_sim + ) + + # Plot - Reference signal + fig, ax = plt.subplots(layout="constrained") + ax.plot(time, phi_ref) + ax.set_ylabel("$\\phi_r$ (deg)") + ax.set_xlabel("$t$ (s)") + ax.grid(linestyle="--") + + # Plot - Disturbance signals + fig, ax = plt.subplots(2, layout="constrained", sharex=True) + ax[0].plot(time, beta_w) + ax[1].plot(time, p_w) + ax[0].set_ylabel("$\\beta_w$ (deg)") + ax[0].grid(linestyle="--") + ax[1].set_ylabel("$p_w$ (deg)") + ax[1].grid(linestyle="--") + ax[-1].set_xlabel("$t$ (s)") + fig.align_ylabels() + + # Plot - Noise signals + fig, ax = plt.subplots(4, layout="constrained", sharex=True) + ax[0].plot(time, noise_phi) + ax[1].plot(time, noise_beta) + ax[2].plot(time, noise_p) + ax[3].plot(time, noise_r) + ax[0].set_ylabel("$n_{\\phi}$ (deg)") + ax[1].set_ylabel("$n_{\\beta}$ (deg)") + ax[2].set_ylabel("$n_{p}$ (deg/s)") + ax[3].set_ylabel("$n_{r}$ (deg/s)") + ax[-1].set_xlabel("$t$ (s)") + for a in ax: + a.grid(linestyle="--") + fig.align_ylabels() + + # Plot - Output signals + fig, ax = plt.subplots(4, layout="constrained", sharex=True) + # Reference signal + ax[0].plot(time, phi_ref, color="black", linestyle="--", label="Reference") + for idx_sim, trd_sim in enumerate(reversed(trd_sim_list)): + # Output signals + phi = trd_sim.y[0, :] + beta = trd_sim.y[1, :] + p = trd_sim.y[2, :] + r = trd_sim.y[3, :] + # Plot signals (nominal) + if idx_sim == len(trd_sim_list) - 1: + ax[0].plot(time, phi, color="tab:blue", label="Nominal") + ax[1].plot(time, beta, color="tab:blue", label="Nominal") + ax[2].plot(time, p, color="tab:blue", label="Nominal") + ax[3].plot(time, r, color="tab:blue", label="Nominal") + # Plot signals (off-nominal) + else: + ax[0].plot(time, phi, color="tab:orange", alpha=0.25, label="Off-Nominal") + ax[1].plot(time, beta, color="tab:orange", alpha=0.25, label="Off-Nominal") + ax[2].plot(time, p, color="tab:orange", alpha=0.25, label="Off-Nominal") + ax[3].plot(time, r, color="tab:orange", alpha=0.25, label="Off-Nominal") + ax[0].set_ylabel("$\\phi$ (deg)") + ax[1].set_ylabel("$\\beta$ (deg)") + ax[2].set_ylabel("$p$ (deg/s)") + ax[3].set_ylabel("$r$ (deg/s)") + ax[-1].set_xlabel("$t$ (s)") + for a in ax: + a.grid(linestyle="--") + handles, labels = ax[0].get_legend_handles_labels() + legend_dict = dict(zip(labels, handles)) + fig.align_ylabels() + fig.legend( + handles=legend_dict.values(), + labels=legend_dict.keys(), + loc="outside lower center", + ncol=3, + ) + + # Plot - Input signals + fig, ax = plt.subplots(4, layout="constrained", sharex=True) + for idx_sim, trd_sim in enumerate(reversed(trd_sim_list)): + # Control signals + zeta = trd_sim.y[6, :] + xi = trd_sim.y[7, :] + rate_zeta = trd_sim.y[8, :] + rate_xi = trd_sim.y[9, :] + if idx_sim == len(trd_sim_list) - 1: + ax[0].plot( + time, zeta, color="tab:blue", linestyle="-", label="Actuator (Nominal)" + ) + ax[1].plot( + time, + rate_zeta, + color="tab:blue", + linestyle="-", + label="Actuator (Nominal)", + ) + ax[2].plot( + time, xi, color="tab:blue", linestyle="-", label="Actuator (Nominal)" + ) + ax[3].plot( + time, + rate_xi, + color="tab:blue", + linestyle="-", + label="Actuator (Nominal)", + ) + else: + ax[0].plot( + time, + zeta, + color="tab:orange", + alpha=0.25, + linestyle="-", + label="Actuator (Off-Nominal)", + ) + ax[1].plot( + time, + rate_zeta, + color="tab:orange", + alpha=0.25, + linestyle="-", + label="Actuator (Off-Nominal)", + ) + ax[2].plot( + time, + xi, + color="tab:orange", + alpha=0.25, + linestyle="-", + label="Actuator (Off-Nominal)", + ) + ax[3].plot( + time, + rate_xi, + color="tab:orange", + alpha=0.25, + linestyle="-", + label="Actuator (Off-Nominal)", + ) + ax[0].set_ylabel("$\\zeta$ (deg)") + ax[1].set_ylabel("$\\xi$ (deg)") + ax[2].set_ylabel("$\\dot{\\zeta}$ (deg/s)") + ax[3].set_ylabel("$\\dot{\\xi}$ (deg/s)") + ax[-1].set_xlabel("$t$ (s)") + for a in ax: + a.grid(linestyle="--") + handles, labels = ax[0].get_legend_handles_labels() + legend_dict = dict(zip(labels, handles)) + fig.align_ylabels() + fig.legend( + handles=legend_dict.values(), + labels=legend_dict.keys(), + loc="outside lower center", + ncol=2, + ) + plt.show() + + +if __name__ == "__main__": + example_dk_iter_list_order_aircraft() diff --git a/pyproject.toml b/pyproject.toml index 93923f9..e107e49 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "hatchling.build" [project] name = "dkpy" -version = "0.1.10" +version = "0.1.11" dependencies = [ "numpy>=1.21.0", "scipy>=1.7.0", diff --git a/src/dkpy/d_scale_fit.py b/src/dkpy/d_scale_fit.py index cecb99d..9649dda 100644 --- a/src/dkpy/d_scale_fit.py +++ b/src/dkpy/d_scale_fit.py @@ -24,7 +24,8 @@ class DScaleFit(metaclass=abc.ABCMeta): def fit( self, omega: np.ndarray, - D_omega: np.ndarray, + D_l_omega: np.ndarray, + D_r_omega: np.ndarray, order: Union[int, np.ndarray] = 0, block_structure: Optional[ Union[ @@ -40,7 +41,10 @@ def fit( ---------- omega : np.ndarray Angular frequencies (rad/s). - D_omega : np.ndarray + D_l_omega : np.ndarray + Transfer matrix evaluated at each frequency, with frequency as last + dimension. + D_r_omega: np.ndarray, Transfer matrix evaluated at each frequency, with frequency as last dimension. order : Union[int, np.ndarray] @@ -82,18 +86,19 @@ class DScaleFitSlicot(DScaleFit): >>> omega = np.logspace(-3, 3, 61) >>> N = P.lft(K) >>> N_omega = N(1j * omega) - >>> mu_omega, D_omega, info = dkpy.SsvLmiBisection().compute_ssv( + >>> mu_omega, D_l_omega, D_r_omega, info = dkpy.SsvLmiBisection().compute_ssv( ... N_omega, ... block_structure, ... ) - >>> D, D_inv = dkpy.DScaleFitSlicot().fit(omega, D_omega, 2, block_structure) + >>> D, D_inv = dkpy.DScaleFitSlicot().fit(omega, D_l_omega, D_r_omega, 2, block_structure) """ def fit( self, omega: np.ndarray, - D_omega: np.ndarray, - order: Union[int, np.ndarray] = 0, + D_l_omega: np.ndarray, + D_r_omega: np.ndarray, + order: Union[int, List[int], np.ndarray] = 0, block_structure: Optional[ Union[ List[uncertainty_structure.UncertaintyBlock], @@ -102,61 +107,132 @@ def fit( ] ] = None, ) -> Tuple[control.StateSpace, control.StateSpace]: + # Check D-scale dimensions + if D_l_omega.shape[0] != D_l_omega.shape[1]: + raise ValueError( + "The left D-scale must be square " + "(D_l_omega.shape[0] = D_l_omega.shape[1])." + ) + if D_r_omega.shape[0] != D_r_omega.shape[1]: + raise ValueError( + "The right D-scale must be square " + "(D_r_omega.shape[0] = D_r_omega.shape[1])." + ) # Get mask if block_structure is None: - mask = -1 * np.ones((D_omega.shape[0], D_omega.shape[1]), dtype=int) + block_structure = [ + uncertainty_structure.ComplexFullBlock( + D_r_omega.shape[0], D_l_omega.shape[0] + ) + ] else: 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 `uncertainty_structure`." - ) + mask_l, mask_r = _generate_d_scale_mask(block_structure) + if isinstance(order, int): + orders_l = order * np.abs(mask_l) + orders_r = order * np.abs(mask_r) + else: + order = np.array(order) + if order.size != len(block_structure): + raise ValueError( + "Length of `order` must be the same as length of `block_structure`." + ) + orders_l_lst = [ + order[i] + * np.ones( + ( + block_structure[i].n_exogenous_outputs, + block_structure[i].n_exogenous_outputs, + ) + ) + for i in range(len(block_structure)) + ] + orders_r_lst = [ + order[i] + * np.ones( + ( + block_structure[i].n_exogenous_inputs, + block_structure[i].n_exogenous_inputs, + ) + ) + for i in range(len(block_structure)) + ] + orders_l = scipy.linalg.block_diag(*orders_l_lst) + orders_r = scipy.linalg.block_diag(*orders_r_lst) + # Transfer matrix - tf_array = np.zeros((D_omega.shape[0], D_omega.shape[1]), dtype=object) - # Fit SISO transfer functions - for row in range(D_omega.shape[0]): - for col in range(D_omega.shape[1]): - if mask[row, col] == 0: - tf_array[row, col] = control.TransferFunction([0], [1], dt=0) - elif mask[row, col] == 1: - if isinstance(order, np.ndarray) and (orders[row, col] != 0): + tf_l_array = np.zeros((D_l_omega.shape[0], D_l_omega.shape[1]), dtype=object) + tf_r_array = np.zeros((D_r_omega.shape[0], D_r_omega.shape[1]), dtype=object) + # Fit SISO transfer functions of left scale + for row in range(D_l_omega.shape[0]): + for col in range(D_l_omega.shape[1]): + if mask_l[row, col] == 0: + tf_l_array[row, col] = control.TransferFunction([0], [1], dt=0) + elif mask_l[row, col] == 1: + if isinstance(order, np.ndarray) and (orders_l[row, col] != 0): warnings.warn( "Entries of `order` in last uncertainty block " "should be 0 since those transfer functions are " "known to be 1. Ignoring value of " f"`order[{row}, {col}]`." ) - tf_array[row, col] = control.TransferFunction([1], [1], dt=0) + tf_l_array[row, col] = control.TransferFunction([1], [1], dt=0) else: n, A, B, C, D = slycot.sb10yd( discfl=0, # Continuous-time flag=1, # Constrain stable, minimum phase lendat=omega.shape[0], - rfrdat=np.real(D_omega[row, col, :]), - ifrdat=np.imag(D_omega[row, col, :]), + rfrdat=np.real(D_l_omega[row, col, :]), + ifrdat=np.imag(D_l_omega[row, col, :]), omega=omega, - n=orders[row, col], + n=orders_l[row, col], tol=0, # Length of cache array ) sys = control.StateSpace(A, B, C, D, dt=0) - tf_array[row, col] = control.ss2tf(sys) - tf = control.combine_tf(tf_array) - ss = control.tf2ss(tf) - ss_inv = _invert_biproper_ss(ss) - return ss, ss_inv + tf_l_array[row, col] = control.ss2tf(sys) + # Fit SISO transfer functions of right scale + for row in range(D_r_omega.shape[0]): + for col in range(D_r_omega.shape[1]): + if mask_r[row, col] == 0: + tf_r_array[row, col] = control.TransferFunction([0], [1], dt=0) + elif mask_r[row, col] == 1: + if isinstance(order, np.ndarray) and (orders_r[row, col] != 0): + warnings.warn( + "Entries of `order` in last uncertainty block " + "should be 0 since those transfer functions are " + "known to be 1. Ignoring value of " + f"`order[{row}, {col}]`." + ) + tf_r_array[row, col] = control.TransferFunction([1], [1], dt=0) + else: + n, A, B, C, D = slycot.sb10yd( + discfl=0, # Continuous-time + flag=1, # Constrain stable, minimum phase + lendat=omega.shape[0], + rfrdat=np.real(D_r_omega[row, col, :]), + ifrdat=np.imag(D_r_omega[row, col, :]), + omega=omega, + n=orders_r[row, col], + tol=0, # Length of cache array + ) + sys = control.StateSpace(A, B, C, D, dt=0) + tf_r_array[row, col] = control.ss2tf(sys) + + tf_l = control.combine_tf(tf_l_array) + tf_r = control.combine_tf(tf_r_array) + ss_l = control.tf2ss(tf_l) + ss_r = control.tf2ss(tf_r) + ss_r_inv = _invert_biproper_ss(ss_r) + return ss_l, ss_r_inv def _generate_d_scale_mask( block_structure: List[uncertainty_structure.UncertaintyBlock], -) -> np.ndarray: +) -> Tuple[np.ndarray, 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 @@ -169,28 +245,40 @@ def _generate_d_scale_mask( Returns ------- - np.ndarray + Tuple[np.ndarray, np.ndarray] Array of integers indicating zero, one, and unknown elements in the block structure. """ num_blocks = len(block_structure) - X_lst = [] + idx_last_full_block = -1 + for idx, block in enumerate(reversed(block_structure)): + if isinstance(block, uncertainty_structure.ComplexFullBlock): + idx_last_full_block = len(block_structure) - 1 - idx + break + mask_l_lst = [] + mask_r_lst = [] 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.is_diagonal: - raise NotImplementedError("Diagonal perturbations are not yet supported.") - if not block.is_square: - raise NotImplementedError("Nonsquare perturbations are not yet supported.") - # Set last scaling to identity - if i == num_blocks - 1: - X_lst.append(np.eye(block.num_inputs, dtype=int)) - else: - X_lst.append(-1 * np.eye(block.num_inputs, dtype=int)) - X = scipy.linalg.block_diag(*X_lst) - return X + if (i == idx_last_full_block) and (not block.is_diagonal): + # Last scaling is always identity if it is a full perturbation + mask_l_lst.append(np.eye(block.n_exogenous_outputs, dtype=int)) + mask_r_lst.append(np.eye(block.n_exogenous_inputs, dtype=int)) + elif (not block.is_complex) and (not block.is_diagonal): + raise NotImplementedError("Real full perturbations are not supported.") + elif (not block.is_complex) and (block.is_diagonal): + raise NotImplementedError( + "Real diagonal perturbations are not yet supported." + ) + elif (block.is_complex) and (block.is_diagonal): + mask_l_lst.append(-1 * np.tri(block.n_exogenous_outputs, dtype=int).T) + mask_r_lst.append(-1 * np.tri(block.n_exogenous_outputs, dtype=int).T) + elif (block.is_complex) and (not block.is_diagonal): + mask_l_lst.append(-1 * np.eye(block.n_exogenous_outputs, dtype=int)) + mask_r_lst.append(-1 * np.eye(block.n_exogenous_inputs, dtype=int)) + mask_l = scipy.linalg.block_diag(*mask_l_lst) + mask_r = scipy.linalg.block_diag(*mask_r_lst) + return mask_l, mask_r def _invert_biproper_ss(ss: control.StateSpace) -> control.StateSpace: diff --git a/src/dkpy/dk_iteration.py b/src/dkpy/dk_iteration.py index a522ea3..833ac81 100644 --- a/src/dkpy/dk_iteration.py +++ b/src/dkpy/dk_iteration.py @@ -25,7 +25,6 @@ d_scale_fit, structured_singular_value, uncertainty_structure, - utilities, ) @@ -43,10 +42,13 @@ def __init__( self, omega: np.ndarray, mu_omega: np.ndarray, - D_omega: np.ndarray, + D_l_omega: np.ndarray, + D_r_inv_omega: np.ndarray, mu_fit_omega: np.ndarray, - D_fit_omega: np.ndarray, - D_fit: control.StateSpace, + D_l_fit_omega: np.ndarray, + D_r_inv_fit_omega: np.ndarray, + D_l_fit: control.StateSpace, + D_r_inv_fit: control.StateSpace, block_structure: Union[ List[uncertainty_structure.UncertaintyBlock], List[List[int]], np.ndarray ], @@ -59,23 +61,32 @@ def __init__( Angular frequencies to evaluate D-scales (rad/s). mu_omega : np.ndarray Numerically computed structured singular value at each frequency. - D_omega : np.ndarray - Numerically computed D-scale magnitude at each frequency. + D_l_omega : np.ndarray + Numerically computed left D-scale magnitude at each frequency. + D_r_inv_omega : np.ndarray + Numerically computed inverse right D-scale magnitude at each frequency. mu_fit_omega : np.ndarray Fit structured singular value at each frequency. - D_fit_omega : np.ndarray + D_l_fit_omega : np.ndarray + Fit left D-scale magnitude at each frequency. + D_r_inv_fit_omega : np.ndarray Fit D-scale magnitude at each frequency. - D_fit : control.StateSpace - Fit D-scale state-space representation. + D_l_fit_omega : np.ndarray + Fit left D-scale magnitude at each frequency. + D_r_inv_fit : control.StateSpace + Fit right D-scale state-space representation. uncertainty_structure : Union[List[uncertainty_structure.UncertaintyBlock], List[List[int]], np.ndarray] Uncertainty block structure representation. """ self.omega = omega self.mu_omega = mu_omega - self.D_omega = D_omega + self.D_l_omega = D_l_omega + self.D_r_inv_omega = D_r_inv_omega self.mu_fit_omega = mu_fit_omega - self.D_fit_omega = D_fit_omega - self.D_fit = D_fit + self.D_l_fit_omega = D_l_fit_omega + self.D_r_inv_fit_omega = D_r_inv_fit_omega + self.D_l_fit = D_l_fit + self.D_r_inv_fit = D_r_inv_fit self.block_structure = ( uncertainty_structure._convert_block_structure_representation( block_structure @@ -87,11 +98,12 @@ def create_from_fit( cls, omega: np.ndarray, mu_omega: np.ndarray, - D_omega: np.ndarray, + D_l_omega: np.ndarray, + D_r_omega: np.ndarray, P: control.StateSpace, K: control.StateSpace, - D_fit: control.StateSpace, - D_fit_inv: control.StateSpace, + D_l_fit: control.StateSpace, + D_r_inv_fit: control.StateSpace, block_structure: Union[ List[uncertainty_structure.UncertaintyBlock], List[List[int]], np.ndarray ], @@ -104,16 +116,18 @@ def create_from_fit( Angular frequencies to evaluate D-scales (rad/s). mu_omega : np.ndarray Numerically computed structured singular value at each frequency. - D_omega : np.ndarray - Numerically computed D-scale magnitude at each frequency. + D_l_omega : np.ndarray + Numerically computed left D-scale magnitude at each frequency. + D_r_omega : np.ndarray + Numerically computed right D-scale magnitude at each frequency. P : control.StateSpace Generalized plant. K : control.StateSpace Controller. - D_fit : control.StateSpace - Fit D-scale magnitude at each frequency. - D_fit_inv : control.StateSpace - Fit inverse D-scale magnitude at each frequency. + D_l_fit : control.StateSpace + Fit left D-scale magnitude at each frequency. + D_r_inv_fit : control.StateSpace + Fit inverse right D-scale magnitude at each frequency. block_structure : Union[List[uncertainty_structure.UncertaintyBlock], List[List[int]], np.ndarray] Uncertainty block structure representation. @@ -135,26 +149,31 @@ def create_from_fit( >>> omega = np.logspace(-3, 3, 61) >>> N = P.lft(K) >>> N_omega = N(1j * omega) - >>> mu_omega, D_omega, info = dkpy.SsvLmiBisection().compute_ssv( + >>> mu_omega, D_l_omega, D_r_omega, info = dkpy.SsvLmiBisection().compute_ssv( ... N_omega, ... block_structure, ... ) - >>> D, D_inv = dkpy.DScaleFitSlicot().fit(omega, D_omega, 2, block_structure) + >>> D_l, D_r_inv = dkpy.DScaleFitSlicot().fit(omega, D_l_omega, D_r_omega, 2, block_structure) >>> d_scale_fit_info = IterResult.create_from_fit( ... omega, ... mu_omega, - ... D_omega, + ... D_l_omega, + ... D_r_omega, ... P, ... K, - ... D, - ... D_inv, + ... D_l, + ... D_r_inv, ... block_structure, ... ) """ + # Compute `D_r^{-1}(omega)` from `D_r(omega)` + D_r_inv_omega = np.array( + [scipy.linalg.inv(D_r_omega[:, :, i]) for i in range(D_r_omega.shape[2])] + ).transpose(1, 2, 0) # Compute ``mu(omega)`` based on fit D-scales N = P.lft(K) - scaled_cl = (D_fit * N * D_fit_inv)(1j * omega) + scaled_cl = (D_l_fit * N * D_r_inv_fit)(1j * omega) mu_fit_omega = np.array( [ np.max(scipy.linalg.svdvals(scaled_cl[:, :, i])) @@ -162,14 +181,18 @@ def create_from_fit( ] ) # Compute ``D(omega)`` based on fit D-scales - D_fit_omega = D_fit(1j * omega) + D_l_fit_omega = D_l_fit(1j * omega) + D_r_inv_fit_omega = D_r_inv_fit(1j * omega) return cls( omega, mu_omega, - D_omega, + D_l_omega, + D_r_inv_omega, mu_fit_omega, - D_fit_omega, - D_fit, + D_l_fit_omega, + D_r_inv_fit_omega, + D_l_fit, + D_r_inv_fit, block_structure, ) @@ -270,9 +293,11 @@ def synthesize( ) N = P.lft(K) N_omega = N(1j * omega) - mu_omega, D_omega, info = self.structured_singular_value.compute_ssv( - N_omega, - block_structure=block_structure, + mu_omega, D_l_omega, D_r_omega, info = ( + self.structured_singular_value.compute_ssv( + N_omega, + block_structure=block_structure, + ) ) # Start iteration while True: @@ -282,7 +307,8 @@ def synthesize( iteration, omega, mu_omega, - D_omega, + D_l_omega, + D_r_omega, P, K, block_structure, @@ -292,9 +318,10 @@ def synthesize( self._log.info("Iteration complete") break # Fit transfer functions to gridded D-scales - D_fit, D_fit_inv = self.d_scale_fit.fit( + D_l_fit, D_r_fit_inv = self.d_scale_fit.fit( omega, - D_omega, + D_l_omega, + D_r_omega, order=fit_order, block_structure=block_structure, ) @@ -303,33 +330,37 @@ def synthesize( IterResult.create_from_fit( omega, mu_omega, - D_omega, + D_l_omega, + D_r_omega, P, K, - D_fit, - D_fit_inv, + D_l_fit, + D_r_fit_inv, block_structure, ) ) # Augment D-scales with identity transfer functions D_aug, D_aug_inv = _augment_d_scales( - D_fit, - D_fit_inv, + D_l_fit, + D_r_fit_inv, n_y=n_y, n_u=n_u, ) + plant_gen_scale = D_aug * P * D_aug_inv # Synthesize controller K, _, gamma, info = self.controller_synthesis.synthesize( - D_aug * P * D_aug_inv, + plant_gen_scale, 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, + mu_omega, D_l_omega, D_r_omega, info = ( + self.structured_singular_value.compute_ssv( + N_omega, + block_structure=block_structure, + ) ) # Increment iteration iteration += 1 @@ -340,7 +371,8 @@ def _get_fit_order( iteration: int, omega: np.ndarray, mu_omega: np.ndarray, - D_omega: np.ndarray, + D_l_omega: np.ndarray, + D_r_omega: np.ndarray, P: control.StateSpace, K: control.StateSpace, block_structure: List[uncertainty_structure.UncertaintyBlock], @@ -355,8 +387,10 @@ def _get_fit_order( Angular frequencies to evaluate D-scales (rad/s). mu_omega : np.ndarray Numerically computed structured singular value at each frequency. - D_omega : np.ndarray - Numerically computed D-scale magnitude at each frequency. + D_l_omega : np.ndarray + Numerically computed left D-scale magnitude at each frequency. + D_r_omega : np.ndarray + Numerically computed right D-scale magnitude at each frequency. P : control.StateSpace Generalized plant. K : control.StateSpace @@ -437,7 +471,8 @@ def _get_fit_order( iteration: int, omega: np.ndarray, mu_omega: np.ndarray, - D_omega: np.ndarray, + D_l_omega: np.ndarray, + D_r_omega: np.ndarray, P: control.StateSpace, K: control.StateSpace, block_structure: List[uncertainty_structure.UncertaintyBlock], @@ -508,7 +543,8 @@ def _get_fit_order( iteration: int, omega: np.ndarray, mu_omega: np.ndarray, - D_omega: np.ndarray, + D_l_omega: np.ndarray, + D_r_omega: np.ndarray, P: control.StateSpace, K: control.StateSpace, block_structure: List[uncertainty_structure.UncertaintyBlock], @@ -594,7 +630,8 @@ def _get_fit_order( iteration: int, omega: np.ndarray, mu_omega: np.ndarray, - D_omega: np.ndarray, + D_l_omega: np.ndarray, + D_r_omega: np.ndarray, P: control.StateSpace, K: control.StateSpace, block_structure: List[uncertainty_structure.UncertaintyBlock], @@ -610,20 +647,22 @@ def _get_fit_order( fit_order = 0 relative_errors = [] while True: - D_fit, D_fit_inv = self.d_scale_fit.fit( + D_l_fit, D_r_fit_inv = self.d_scale_fit.fit( omega, - D_omega, + D_l_omega, + D_r_omega, order=fit_order, block_structure=block_structure, ) d_scale_fit_info = IterResult.create_from_fit( omega, mu_omega, - D_omega, + D_l_omega, + D_r_omega, P, K, - D_fit, - D_fit_inv, + D_l_fit, + D_r_fit_inv, block_structure, ) error = np.abs(mu_omega - d_scale_fit_info.mu_fit_omega) @@ -637,6 +676,7 @@ def _get_fit_order( if relative_error >= self.max_mu_fit_error: fit_order += 1 else: + best_order = fit_order self._log.info( f"Reached structured singular value target with order {best_order}" ) @@ -701,7 +741,8 @@ def _get_fit_order( iteration: int, omega: np.ndarray, mu_omega: np.ndarray, - D_omega: np.ndarray, + D_l_omega: np.ndarray, + D_r_omega: np.ndarray, P: control.StateSpace, K: control.StateSpace, block_structure: List[uncertainty_structure.UncertaintyBlock], @@ -710,7 +751,8 @@ def _get_fit_order( for fit_order in range(self.max_fit_order + 1): D_fit, D_fit_inv = self.d_scale_fit.fit( omega, - D_omega, + D_l_omega, + D_r_omega, order=fit_order, block_structure=block_structure, ) @@ -718,7 +760,8 @@ def _get_fit_order( IterResult.create_from_fit( omega, mu_omega, - D_omega, + D_l_omega, + D_r_omega, P, K, D_fit, @@ -793,19 +836,20 @@ def plot_mu( >>> omega = np.logspace(-3, 3, 61) >>> N = P.lft(K) >>> N_omega = N(1j * omega) - >>> mu_omega, D_omega, info = dkpy.SsvLmiBisection().compute_ssv( + >>> mu_omega, D_l_omega, D_r_omega, info = dkpy.SsvLmiBisection().compute_ssv( ... N_omega, ... block_structure, ... ) - >>> D, D_inv = dkpy.DScaleFitSlicot().fit(omega, D_omega, 2, block_structure) + >>> D_l, D_r_inv = dkpy.DScaleFitSlicot().fit(omega, D_l_omega, D_r_omega, 2, block_structure) >>> d_scale_fit_info = IterResult.create_from_fit( ... omega, ... mu_omega, - ... D_omega, + ... D_l_omega, + ... D_r_omega, ... P, ... K, - ... D, - ... D_inv, + ... D_l, + ... D_r_inv, ... block_structure, ... ) >>> fig, ax = dkpy.plot_mu(d_scale_fit_info) @@ -843,6 +887,9 @@ def plot_mu( # Set axis labels ax.set_xlabel(r"$\omega$ (rad/s)") ax.set_ylabel(r"$\mu(\omega)$") + ax.set_ylim( + (0.75 * np.min(d_scale_info.mu_omega), 1.25 * np.max(d_scale_info.mu_omega)) + ) ax.grid(linestyle="--") ax.legend(loc="lower left") # Return figure and axes @@ -854,6 +901,7 @@ def plot_D( ax: Optional[np.ndarray] = None, plot_kw: Optional[Dict[str, Any]] = None, hide: Optional[str] = None, + plot_inverse: Optional[bool] = False, ) -> Tuple[plt.Figure, np.ndarray]: """Plot D. @@ -868,6 +916,8 @@ def plot_D( hide : Optional[str] Set to ``'D_omega'`` or ``'D_fit_omega'`` to hide either one of those lines. + plot_inverse: Optional[bool] + Plot the inverse right D-scale. Returns ------- @@ -888,30 +938,31 @@ def plot_D( >>> omega = np.logspace(-3, 3, 61) >>> N = P.lft(K) >>> N_omega = N(1j * omega) - >>> mu_omega, D_omega, info = dkpy.SsvLmiBisection().compute_ssv( + >>> mu_omega, D_l_omega, D_r_omega, info = dkpy.SsvLmiBisection().compute_ssv( ... N_omega, ... block_structure, ... ) - >>> D, D_inv = dkpy.DScaleFitSlicot().fit(omega, D_omega, 2, block_structure) + >>> D_l, D_r_inv = dkpy.DScaleFitSlicot().fit(omega, D_l_omega, D_r_omega, 2, block_structure) >>> d_scale_fit_info = IterResult.create_from_fit( ... omega, ... mu_omega, - ... D_omega, + ... D_l_omega, + ... D_r_omega, ... P, ... K, - ... D, - ... D_inv, + ... D_l, + ... D_r_inv, ... block_structure, ... ) >>> fig, ax = dkpy.plot_D(d_scale_fit_info) """ block_structure = d_scale_info.block_structure - mask = d_scale_fit._generate_d_scale_mask(block_structure) + mask_l, mask_r = d_scale_fit._generate_d_scale_mask(block_structure) # Create figure if not provided if ax is None: fig, ax = plt.subplots( - mask.shape[0], - mask.shape[1], + mask_l.shape[0], + mask_l.shape[1], constrained_layout=True, ) else: @@ -926,37 +977,71 @@ def plot_D( _ = plot_kw.pop("ls", None) _ = plot_kw.pop("linestyle", None) # Plot D - mag_D_omega = np.abs(d_scale_info.D_omega) - mag_D_fit_omega = np.abs(d_scale_info.D_fit_omega) - for i in range(ax.shape[0]): - for j in range(ax.shape[1]): - if mask[i, j] != 0: - dB_omega = 20 * np.log10(mag_D_omega[i, j, :]) - dB_fit_omega = 20 * np.log10(mag_D_fit_omega[i, j, :]) - if hide != "D_omega": - ax[i, j].semilogx( - d_scale_info.omega, - dB_omega, - label=label_D_omega, - ls="--", - **plot_kw, - ) - if hide != "D_fit_omega": - ax[i, j].semilogx( - d_scale_info.omega, - dB_fit_omega, - label=label_D_fit_omega, - **plot_kw, - ) - # Set axis labels - ax[i, j].set_xlabel(r"$\omega$ (rad/s)") - ax[i, j].set_ylabel(rf"$D_{{{i}{j}}}(\omega) (dB)$") - ax[i, j].grid(linestyle="--") - else: - ax[i, j].axis("off") - fig.legend(handles=ax[0, 0].get_lines(), loc="lower left") - # Return figure and axes - return fig, ax + if not plot_inverse: + mag_D_omega = np.abs(d_scale_info.D_l_omega) + mag_D_fit_omega = np.abs(d_scale_info.D_l_fit_omega) + for i in range(ax.shape[0]): + for j in range(ax.shape[1]): + if mask_l[i, j] != 0: + dB_omega = 20 * np.log10(mag_D_omega[i, j, :]) + dB_fit_omega = 20 * np.log10(mag_D_fit_omega[i, j, :]) + if hide != "D_omega": + ax[i, j].semilogx( + d_scale_info.omega, + dB_omega, + label=label_D_omega, + ls="--", + **plot_kw, + ) + if hide != "D_fit_omega": + ax[i, j].semilogx( + d_scale_info.omega, + dB_fit_omega, + label=label_D_fit_omega, + **plot_kw, + ) + # Set axis labels + ax[i, j].set_xlabel(r"$\omega$ (rad/s)") + ax[i, j].set_ylabel(rf"$D_{{{i}{j}}}(\omega) (dB)$") + ax[i, j].grid(linestyle="--") + else: + ax[i, j].axis("off") + fig.legend(handles=ax[0, 0].get_lines(), loc="lower left") + # Return figure and axes + return fig, ax + # Plot D_inv + else: + mag_D_inv_omega = np.abs(d_scale_info.D_r_inv_omega) + mag_D_inv_fit_omega = np.abs(d_scale_info.D_r_inv_fit_omega) + for i in range(ax.shape[0]): + for j in range(ax.shape[1]): + if mask_l[i, j] != 0: + dB_omega = 20 * np.log10(mag_D_inv_omega[i, j, :]) + dB_fit_omega = 20 * np.log10(mag_D_inv_fit_omega[i, j, :]) + if hide != "D_omega": + ax[i, j].semilogx( + d_scale_info.omega, + dB_omega, + label=label_D_omega, + ls="--", + **plot_kw, + ) + if hide != "D_fit_omega": + ax[i, j].semilogx( + d_scale_info.omega, + dB_fit_omega, + label=label_D_fit_omega, + **plot_kw, + ) + # Set axis labels + ax[i, j].set_xlabel(r"$\omega$ (rad/s)") + ax[i, j].set_ylabel(rf"$D_{{{i}{j}}}(\omega)^{{-1}} (dB)$") + ax[i, j].grid(linestyle="--") + else: + ax[i, j].axis("off") + fig.legend(handles=ax[0, 0].get_lines(), loc="lower left") + # Return figure and axes + return fig, ax def _augment_d_scales( diff --git a/src/dkpy/structured_singular_value.py b/src/dkpy/structured_singular_value.py index eba8c7b..c38a0a1 100644 --- a/src/dkpy/structured_singular_value.py +++ b/src/dkpy/structured_singular_value.py @@ -25,10 +25,14 @@ class StructuredSingularValue(metaclass=abc.ABCMeta): def compute_ssv( self, N_omega: np.ndarray, - block_structure: Union[ - List[uncertainty_structure.UncertaintyBlock], List[List[int]], np.ndarray - ], - ) -> Tuple[np.ndarray, np.ndarray, Dict[str, Any]]: + block_structure: Optional[ + Union[ + List[uncertainty_structure.UncertaintyBlock], + List[List[int]], + np.ndarray, + ] + ] = None, + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, Dict[str, Any]]: """Compute structured singular value. Parameters @@ -37,6 +41,8 @@ def compute_ssv( Closed-loop transfer function evaluated at each frequency. block_structure : Union[List[uncertainty_structure.UncertaintyBlock], List[List[int]], np.ndarray] Uncertainty block structure representation. Returns + + Returns ------- Tuple[np.ndarray, np.ndarray, Dict[str, Any]] Structured singular value at each frequency, D-scales at each @@ -66,7 +72,7 @@ class SsvLmiBisection(StructuredSingularValue): >>> omega = np.logspace(-3, 3, 61) >>> N = P.lft(K) >>> N_omega = N(1j * omega) - >>> mu_omega, D_omega, info = dkpy.SsvLmiBisection(n_jobs=None).compute_ssv( + >>> mu_omega, D_l_omega, D_r_omega, info = dkpy.SsvLmiBisection(n_jobs=None).compute_ssv( ... N_omega, ... block_structure, ... ) @@ -138,10 +144,14 @@ def __init__( def compute_ssv( self, N_omega: np.ndarray, - block_structure: Union[ - List[uncertainty_structure.UncertaintyBlock], List[List[int]], np.ndarray - ], - ) -> Tuple[np.ndarray, np.ndarray, Dict[str, Any]]: + block_structure: Optional[ + Union[ + List[uncertainty_structure.UncertaintyBlock], + List[List[int]], + np.ndarray, + ] + ] = None, + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, Dict[str, Any]]: # Solver settings solver_params = ( { @@ -167,13 +177,27 @@ def compute_ssv( else: raise ValueError("`objective` must be `'minimize'` or `'constant'`.") # Convert uncertainty block structure representation - block_structure = uncertainty_structure._convert_block_structure_representation( - block_structure - ) + if block_structure is None: + block_structure = [ + uncertainty_structure.ComplexFullBlock( + N_omega.shape[1], N_omega.shape[0] + ) + ] + else: + block_structure = ( + uncertainty_structure._convert_block_structure_representation( + block_structure + ) + ) def _ssv_at_omega( N_omega: np.ndarray, - ) -> Tuple[Union[float, None], Union[np.ndarray, None], Dict[str, Any]]: + ) -> Tuple[ + Union[float, None], + Union[np.ndarray, 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 @@ -187,10 +211,10 @@ def _ssv_at_omega( Returns ------- Tuple[float, np.ndarray, Dict[str, Any]] - Structured singular value, D-scaling, 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. + Structured singular value, left and right D-scaling, 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. Raises ------ @@ -200,21 +224,23 @@ def _ssv_at_omega( info = {} # Get an optimization variable that shares its block structure with # the D-scalings - X = _generate_ssv_variable(block_structure) + X_l, X_r = _generate_ssv_variable(block_structure) # Set objective function if constant_objective: objective = cvxpy.Minimize(1) else: - # `X` should have a real diagonal because it's Hermitian, but - # CVXPY does not realize this - objective = cvxpy.Minimize(cvxpy.real(cvxpy.trace(X))) + objective = cvxpy.Minimize( + cvxpy.real(cvxpy.trace(X_l)) + cvxpy.real(cvxpy.trace(X_r)) + ) # Set upper bound on structured singular value squared as a parameter gamma_sq = cvxpy.Parameter(1, name="gamma_sq") # Set up the constraints constraints = [ - X.conj().T == X, - X >> lmi_strictness, - N_omega.conj().T @ X @ N_omega - gamma_sq * X << -lmi_strictness, + X_l.conj().T == X_l, + X_r.conj().T == X_r, + X_l >> lmi_strictness, + X_r >> lmi_strictness, + N_omega.conj().T @ X_l @ N_omega - gamma_sq * X_r << -lmi_strictness, ] # Put everything together in the optimization problem problem = cvxpy.Problem(objective, constraints) @@ -250,7 +276,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 + return None, None, None, info # Check if ``gamma`` was increased. if gamma_high > self.initial_guess: warnings.warn( @@ -295,15 +321,16 @@ 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, None, info # Save info info["status"] = "Bisection succeeded." info["gammas"] = gammas info["solver_stats"] = [p.solver_stats for p in problems] info["size_metrics"] = [p.size_metrics for p in problems] info["results"] = results info["iterations"] = n_iterations - D_omega = scipy.linalg.cholesky(X.value) - return (gammas[-1], D_omega, info) + D_l_omega = scipy.linalg.cholesky(X_l.value) + D_r_omega = scipy.linalg.cholesky(X_r.value) + return (gammas[-1], D_l_omega, D_r_omega, info) # Compute structured singular value and D scales for each frequency joblib_results = joblib.Parallel(n_jobs=self.n_jobs)( @@ -313,19 +340,20 @@ def _ssv_at_omega( ] ) # Extract and return results - mu_lst, D_scales_lst, info_lst = tuple(zip(*joblib_results)) + mu_lst, D_l_scales_lst, D_r_scales_lst, info_lst = tuple(zip(*joblib_results)) mu = np.array(mu_lst) - D_scales = np.moveaxis(np.array(D_scales_lst), 0, 2) + D_l_scales = np.array(D_l_scales_lst).transpose(1, 2, 0) + D_r_scales = np.array(D_r_scales_lst).transpose(1, 2, 0) info = {k: [d[k] for d in info_lst] for k in info_lst[0].keys()} info["solver_params"] = solver_params info["lmi_strictness"] = lmi_strictness info["constant_objective"] = constant_objective - return mu, D_scales, info + return mu, D_l_scales, D_r_scales, info def _generate_ssv_variable( block_structure: List[uncertainty_structure.UncertaintyBlock], -) -> cvxpy.Variable: +) -> Tuple[cvxpy.Variable, cvxpy.Variable]: """Get structured singular value optimization variable for the block structure. Parameters @@ -335,50 +363,60 @@ def _generate_ssv_variable( Returns ------- - cvxpy.Variable - CVXPY variable with specified block structure. + Tuple[cvxpy.Variable, cvxpy.Variable] + CVXPY variables with specified block structure. """ num_blocks = len(block_structure) - X_lst = [] + idx_last_full_block = -1 + for idx, block in enumerate(reversed(block_structure)): + if isinstance(block, uncertainty_structure.ComplexFullBlock): + idx_last_full_block = len(block_structure) - 1 - idx + break + X_l_lst = [] + X_r_lst = [] for i in range(num_blocks): - row = [] + row_l = [] + row_r = [] 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 (not block_i.is_complex) and (not block_i.is_diagonal): + if (i == idx_last_full_block) and (not block_i.is_diagonal): + # Last scaling is always identity if it is a full perturbation + row_l.append(np.eye(block_i.n_exogenous_outputs)) + row_r.append(np.eye(block_i.n_exogenous_inputs)) + elif (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): + elif (not block_i.is_complex) and (block_i.is_diagonal): raise NotImplementedError( "Real diagonal perturbations are not yet supported." ) - if (block_i.is_complex) and (block_i.is_diagonal): - raise NotImplementedError( - "Complex diagonal perturbations are not yet supported." - ) - if (block_i.is_complex) and (not block_i.is_square): - raise NotImplementedError( - "Nonsquare perturbations are not yet supported." + elif (block_i.is_complex) and (block_i.is_diagonal): + X_i = cvxpy.Variable( + (block_i.n_exogenous_outputs, block_i.n_exogenous_outputs), + hermitian=True, + name=f"X{i}", ) - if i == num_blocks - 1: - # Last scaling is always identity - row.append(np.eye(block_i.num_inputs)) - else: - # Every other scaling is either a scalar or a scalar - # multiplied by identity - 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_i.num_inputs)) + row_l.append(X_i) + row_r.append(X_i) + elif (block_i.is_complex) and (not block_i.is_diagonal): + x_i = cvxpy.Variable((), complex=False, name=f"x{i}") + row_l.append(x_i * np.eye(block_i.n_exogenous_outputs)) + row_r.append(x_i * np.eye(block_i.n_exogenous_inputs)) else: # If off the block diagonal, insert zeros - row.append(np.zeros((block_i.num_inputs, block_j.num_inputs))) - X_lst.append(row) - X = cvxpy.bmat(X_lst) - return X + row_l.append( + np.zeros((block_i.n_exogenous_outputs, block_j.n_exogenous_outputs)) + ) + row_r.append( + np.zeros((block_i.n_exogenous_inputs, block_j.n_exogenous_inputs)) + ) + X_l_lst.append(row_l) + X_r_lst.append(row_r) + X_l = cvxpy.bmat(X_l_lst) + X_r = cvxpy.bmat(X_r_lst) + return X_l, X_r diff --git a/src/dkpy/uncertainty_structure.py b/src/dkpy/uncertainty_structure.py index 4d925f1..5a08a03 100644 --- a/src/dkpy/uncertainty_structure.py +++ b/src/dkpy/uncertainty_structure.py @@ -10,14 +10,14 @@ class UncertaintyBlock(metaclass=abc.ABCMeta): @property @abc.abstractmethod - def num_inputs(self) -> int: - """Get number of inputs of the uncertainty block.""" + def n_exogenous_inputs(self) -> int: + """Get number of exogenous inputs (w) for the uncertainty block.""" raise NotImplementedError() @property @abc.abstractmethod - def num_outputs(self) -> int: - """Get number of output of the uncertainty block.""" + def n_exogenous_outputs(self) -> int: + """Get number of exogenous outputs (z) for the uncertainty block.""" raise NotImplementedError() @property @@ -48,7 +48,8 @@ def __init__(self, num_channels: int): Parameters ---------- num_channels : int - Number of inputs/outputs to the uncertainty block. + Number of exogenous inputs (w) and exogenous outputs (z) for the + uncertainty block. Raises ------ @@ -61,27 +62,23 @@ def __init__(self, num_channels: int): self._num_channels = num_channels @property - def num_inputs(self) -> int: + def n_exogenous_inputs(self) -> int: return self._num_channels @property - def num_outputs(self) -> int: - """Get number of output of the uncertainty block.""" + def n_exogenous_outputs(self) -> int: 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 @@ -94,7 +91,8 @@ def __init__(self, num_channels: int): Parameters ---------- num_channels : int - Number of inputs/outputs to the uncertainty block. + Number of exogenous inputs (w) and exogenous outputs (z) for the + uncertainty block. Raises ------ @@ -107,82 +105,74 @@ def __init__(self, num_channels: int): self._num_channels = num_channels @property - def num_inputs(self) -> int: + def n_exogenous_inputs(self) -> int: return self._num_channels @property - def num_outputs(self) -> int: - """Get number of output of the uncertainty block.""" + def n_exogenous_outputs(self) -> int: 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): + def __init__(self, n_exogenous_inputs: int, n_exogenous_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. + n_exogenous_inputs : int + Number of exogenous inputs (w) for the uncertainty block. + n_exogenous_outputs : int + Number of exogenous outputs (z) for the uncertainty block. Raises ------ ValueError - If ``num_inputs`` is not greater than zero. + If ``n_exogenous_inputs`` is not greater than zero. ValueError - If ``num_outputs`` is not greater than zero. + If ``n_exogenous_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.") + if n_exogenous_inputs <= 0: + raise ValueError("`num_exogenous_inputs` must be greater than 0.") + if n_exogenous_outputs <= 0: + raise ValueError("`num_exogenous_outputs` must be greater than 0.") - self._num_inputs = num_inputs - self._num_outputs = num_outputs + self._n_exogenous_inputs = n_exogenous_inputs + self._n_exogenous_outputs = n_exogenous_outputs @property - def num_inputs(self) -> int: - return self._num_inputs + def n_exogenous_inputs(self) -> int: + return self._n_exogenous_inputs @property - def num_outputs(self) -> int: - """Get number of output of the uncertainty block.""" - return self._num_outputs + def n_exogenous_outputs(self) -> int: + return self._n_exogenous_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 + return self._n_exogenous_inputs == self._n_exogenous_outputs def _convert_block_structure_representation( @@ -247,7 +237,7 @@ def _convert_block_structure_representation( 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]) + ComplexFullBlock(block_matlab[0], block_matlab[1]) ) else: raise ValueError( diff --git a/src/dkpy/utilities.py b/src/dkpy/utilities.py index 6ab2f6f..2d48e31 100644 --- a/src/dkpy/utilities.py +++ b/src/dkpy/utilities.py @@ -126,6 +126,198 @@ def example_skogestad2006_p325() -> Dict[str, Any]: return out +def example_mackenroth2004_p430(): + """Add generalized plant from [M04]_, Section 14.1 (p. 430).""" + # Airframe model + A_af = np.array( + [ + [-0.0869, 0, 0.039, -1], + [-4.424, -1.184, 0, 0.335], + [0, 1, 0, 0], + [2.148, -0.021, 0, -0.228], + ] + ) + B1_af = np.array([[0.0869, 0], [4.424, 1.184], [0, -1], [-2.148, 0.021]]) + B2_af = np.array([[0.0223, 0], [0.547, 2.12], [0, 0], [-1.169, 0.065]]) + B_af = np.block([[B1_af, B2_af]]) + C_af = np.eye(4) + D_af = np.zeros((4, 4)) + airframe = control.StateSpace(A_af, B_af, C_af, D_af, name="airframe") + airframe.set_inputs(["beta_w", "p_w", "zeta", "xi"]) + airframe.set_outputs(["beta", "p", "phi", "r"]) + + # Actuator model + A_act = np.array( + [[0, 1, 0, 0], [-1600, -56, 0, 0], [0, 0, 0, 1], [0, 0, -1600, -56]] + ) + B_act = np.array([[0, 0], [1600, 0], [0, 0], [0, 1600]]) + C_act = np.eye(4) + D_act = np.zeros((4, 2)) + actuator = control.StateSpace(A_act, B_act, C_act, D_act, name="actuator") + actuator.set_inputs(["zeta_unc", "xi_unc"]) + actuator.set_outputs(["zeta", "rate_zeta", "xi", "rate_xi"]) + + # Performance weight + weight_perf_angle = control.TransferFunction([1, 1.5], [1.5, 0.015]) + weight_perf_rate = control.TransferFunction([1, 1.5], [10, 20]) + weight_perf = control.append( + weight_perf_angle, + weight_perf_angle, + weight_perf_rate, + weight_perf_rate, + name="weight_perf", + ) + weight_perf.set_inputs(["e_phi", "beta", "p", "r"]) + weight_perf.set_outputs(weight_perf.noutputs, "z_p") + + # Actuator weight + weight_actu_angle = control.StateSpace([], [], [], [[1 / 20]]) + weight_actu_rate = control.StateSpace([], [], [], [[1 / 60]]) + weight_actu = control.append( + weight_actu_angle, + weight_actu_angle, + weight_actu_rate, + weight_actu_rate, + name="weight_actu", + ) + weight_actu.set_inputs(["zeta_c", "xi_c", "rate_zeta", "rate_xi"]) + weight_actu.set_outputs(weight_actu.noutputs, "z_u") + + # Disturbance weight + weight_dist_channel = control.StateSpace([], [], [], [[1]]) + weight_dist = control.append( + weight_dist_channel, weight_dist_channel, name="weight_dist" + ) + weight_dist.set_inputs(["beta_w_nor", "p_w_nor"]) + weight_dist.set_outputs(["beta_w", "p_w"]) + + # Reference weight + weight_ref = control.TransferFunction([2.25], [1, 2.1, 2.25]) + weight_ref = control.ss(weight_ref, name="weight_ref") + weight_ref.set_inputs("phi_ref_nor") + weight_ref.set_outputs("phi_ref") + + # Noise weight + weight_noise_channel = 0.0005 * control.TransferFunction([1 / 0.02, 1], [1, 1]) + weight_noise = control.append( + weight_noise_channel, + weight_noise_channel, + weight_noise_channel, + weight_noise_channel, + name="weight_noise", + ) + weight_noise.set_inputs(weight_noise.ninputs, "n_nor") + weight_noise.set_outputs(weight_noise.noutputs, "n") + + # Input multiplicative uncertainty weight + weight_unc_channel = control.TransferFunction([0.25, 0.05], [0.125, 1]) + weight_unc = control.append( + weight_unc_channel, + weight_unc_channel, + name="weight_unc", + ) + weight_unc.set_inputs(["zeta_c", "xi_c"]) + weight_unc.set_outputs(weight_unc.noutputs, "y_del") + + # Uncertainty summation junction + sum_uncertainty = control.StateSpace([], [], [], [[1, 0, 1, 0], [0, 1, 0, 1]]) + sum_uncertainty.set_inputs(["zeta_c", "xi_c", "u_del[0]", "u_del[1]"]) + sum_uncertainty.set_outputs(["zeta_unc", "xi_unc"]) + + # Reference tracking error difference junction + sum_error = control.StateSpace([], [], [], [[-1, 1]]) + sum_error.set_inputs(["phi", "phi_ref"]) + sum_error.set_outputs("e_phi") + + # Noise summation junction + sum_noise = control.StateSpace([], [], [], np.block([[np.eye(4), np.eye(4)]])) + sum_noise.set_inputs(["phi", "beta", "p", "r", "n[0]", "n[1]", "n[2]", "n[3]"]) + sum_noise.set_outputs(["phi_meas", "beta_meas", "p_meas", "r_meas"]) + + # Reference signal passthrough (control.interconnect cannot have the interconnected + # system have an output that is not the output of a sub-system) + passthrough_ref = control.StateSpace([], [], [], np.eye(1)) + passthrough_ref.set_inputs(["phi_ref_exo"]) + passthrough_ref.set_outputs(["phi_ref_nor"]) + + # Generalized plant + input_id_list = [ + "u_del[0]", + "u_del[1]", + "phi_ref_exo", + "n_nor[0]", + "n_nor[1]", + "n_nor[2]", + "n_nor[3]", + "beta_w_nor", + "p_w_nor", + "zeta_c", + "xi_c", + ] + output_id_list = [ + "y_del[0]", + "y_del[1]", + "z_p[0]", + "z_p[1]", + "z_p[2]", + "z_p[3]", + "z_u[0]", + "z_u[1]", + "z_u[2]", + "z_u[3]", + "phi_ref_nor", + "phi_meas", + "beta_meas", + "p_meas", + "r_meas", + ] + plant_gen = control.interconnect( + syslist=[ + airframe, + actuator, + weight_perf, + weight_actu, + weight_dist, + weight_ref, + weight_noise, + weight_unc, + sum_uncertainty, + sum_error, + sum_noise, + passthrough_ref, + ], + inplist=input_id_list, + outlist=output_id_list, + name="plant_gen", + ) + plant_gen.set_inputs(input_id_list) + plant_gen.set_outputs(output_id_list) + plant_gen = plant_gen.minreal() + # Dimensions + n_y = 5 + n_u = 2 + n_u_delta = 2 + n_y_delta = 2 + n_w = 7 + n_z = 8 + # Example output dictionary + out = { + "P": plant_gen, + "airframe": airframe, + "actuator": actuator, + "weight_unc": weight_unc, + "sum_noise": sum_noise, + "sum_uncertainty": sum_uncertainty, + "n_u": n_u, + "n_y": n_y, + "n_u_delta": n_u_delta, + "n_y_delta": n_y_delta, + "n_w": n_w, + "n_z": n_z, + } + return out + + def _tf_close_coeff( tf_a: control.TransferFunction, tf_b: control.TransferFunction, diff --git a/tests/test_d_scale_fit.py b/tests/test_d_scale_fit.py index 2bc3a96..3e64642 100644 --- a/tests/test_d_scale_fit.py +++ b/tests/test_d_scale_fit.py @@ -11,76 +11,56 @@ class TestTfFitSlicot: """Test :class:`DScaleFitSlicot`.""" @pytest.mark.parametrize( - "omega, tf, order, block_structure, atol", + "omega, tf_l, tf_r, order, block_structure, atol", [ - ( - np.logspace(-2, 2, 100), - control.TransferFunction([10], [1]), - 0, - None, - 1e-6, - ), - ( - np.logspace(-2, 2, 100), - control.TransferFunction([1, 1], [1, 10]), - 1, - None, - 1e-2, - ), ( np.logspace(-2, 2, 100), control.TransferFunction( [ [ [1, 1], - [1, 1], + [0], ], [ - [1, 1], - [1, 1], + [0], + [1], ], ], [ [ [1, 10], - [1, 9], + [1], ], [ - [1, 8], - [1, 10], + [1], + [1], ], ], ), - 1, - None, - 1e-2, - ), - ( - np.logspace(-2, 2, 100), control.TransferFunction( [ [ - [1, 2, 1], - [1, 2, 1], + [1, 1], + [0], ], [ - [1, 2, 1], - [1, 2, 1], + [0], + [1], ], ], [ [ - [1, 10, 1], - [1, 9, 2], + [1, 10], + [1], ], [ - [1, 8, 3], - [1, 10, 4], + [1], + [1], ], ], ), - 2, - None, + 1, + [dkpy.ComplexFullBlock(1, 1), dkpy.ComplexFullBlock(1, 1)], 1e-2, ), ( @@ -88,31 +68,25 @@ class TestTfFitSlicot: control.TransferFunction( [ [ - [1, 2, 1], - [1, 2, 1], + [1, 1], + [0], ], [ - [1, 2, 1], - [1, 2, 1], + [0], + [1], ], ], [ [ - [1, 10, 1], - [1, 9, 2], + [1, 10], + [1], ], [ - [1, 8, 3], - [1, 10, 4], + [1], + [1], ], ], ), - 2 * np.ones((2, 2)), - None, - 1e-2, - ), - ( - np.logspace(-2, 2, 100), control.TransferFunction( [ [ @@ -135,7 +109,7 @@ class TestTfFitSlicot: ], ], ), - 1, + [1, 0], [dkpy.ComplexFullBlock(1, 1), dkpy.ComplexFullBlock(1, 1)], 1e-2, ), @@ -144,7 +118,7 @@ class TestTfFitSlicot: control.TransferFunction( [ [ - [1, 1], + [1], [0], ], [ @@ -154,7 +128,7 @@ class TestTfFitSlicot: ], [ [ - [1, 10], + [1], [1], ], [ @@ -163,12 +137,6 @@ class TestTfFitSlicot: ], ], ), - np.diag([1, 0]), - [dkpy.ComplexFullBlock(1, 1), dkpy.ComplexFullBlock(1, 1)], - 1e-2, - ), - ( - np.logspace(-2, 2, 100), control.TransferFunction( [ [ @@ -219,12 +187,6 @@ class TestTfFitSlicot: ], ], ), - 1, - np.array([[1, 1], [1, 1]]), - 1e-2, - ), - ( - np.logspace(-2, 2, 100), control.TransferFunction( [ [ @@ -247,7 +209,7 @@ class TestTfFitSlicot: ], ], ), - np.diag([1, 0]), + 1, np.array([[1, 1], [1, 1]]), 1e-2, ), @@ -256,7 +218,7 @@ class TestTfFitSlicot: control.TransferFunction( [ [ - [1], + [1, 1], [0], ], [ @@ -266,7 +228,29 @@ class TestTfFitSlicot: ], [ [ + [1, 10], [1], + ], + [ + [1], + [1], + ], + ], + ), + control.TransferFunction( + [ + [ + [1, 1], + [0], + ], + [ + [0], + [1], + ], + ], + [ + [ + [1, 10], [1], ], [ @@ -275,25 +259,30 @@ class TestTfFitSlicot: ], ], ), - 1, - np.array([[2, 2]]), + [1, 0], + np.array([[1, 1], [1, 1]]), 1e-2, ), ], ) - def test_tf_fit_slicot(self, omega, tf, order, block_structure, atol): + def test_tf_fit_slicot(self, omega, tf_l, tf_r, order, block_structure, atol): """Test :class:`DScaleFitSlicot`.""" - D_omega = tf(1j * omega) - if D_omega.ndim == 1: - D_omega = D_omega.reshape((1, 1, -1)) - tf_fit, _ = dkpy.DScaleFitSlicot().fit(omega, D_omega, order, block_structure) - D_omega_fit = tf_fit(1j * omega) - if D_omega_fit.ndim == 1: - D_omega_fit = D_omega_fit.reshape((1, 1, -1)) - np.testing.assert_allclose(D_omega, D_omega_fit, atol=atol) + D_l_omega = tf_l(1j * omega) + D_r_omega = tf_r(1j * omega) + if D_l_omega.ndim == 1: + D_l_omega = D_l_omega.reshape((1, 1, -1)) + if D_r_omega.ndim == 1: + D_r_omega = D_r_omega.reshape((1, 1, -1)) + tf_l_fit, _ = dkpy.DScaleFitSlicot().fit( + omega, D_l_omega, D_r_omega, order, block_structure + ) + D_l_omega_fit = tf_l_fit(1j * omega) + if D_l_omega_fit.ndim == 1: + D_l_omega_fit = D_l_omega_fit.reshape((1, 1, -1)) + np.testing.assert_allclose(D_l_omega, D_l_omega_fit, atol=atol) @pytest.mark.parametrize( - "omega, tf, order, block_structure", + "omega, tf_l, tf_r, order, block_structure", [ ( np.logspace(-2, 2, 100), @@ -302,28 +291,27 @@ def test_tf_fit_slicot(self, omega, tf, order, block_structure, atol): [ [1, 1], [1, 1], + [1, 1], ], [ [1, 1], [1, 1], + [1, 1], ], ], [ [ [1, 10], - [1, 10], + [1, 9], + [1, 7], ], [ - [1, 10], + [1, 8], + [1, 6], [1, 10], ], ], ), - 1, - None, - ), - ( - np.logspace(-2, 2, 100), control.TransferFunction( [ [ @@ -355,19 +343,22 @@ def test_tf_fit_slicot(self, omega, tf, order, block_structure, atol): ), ], ) - def test_tf_fit_slicot_error(self, omega, tf, order, block_structure): + def test_tf_fit_slicot_error(self, omega, tf_l, tf_r, order, block_structure): with pytest.raises(ValueError): - D_omega = tf(1j * omega) - if D_omega.ndim == 1: - D_omega = D_omega.reshape((1, 1, -1)) + D_l_omega = tf_l(1j * omega) + D_r_omega = tf_r(1j * omega) + if D_l_omega.ndim == 1: + D_l_omega = D_l_omega.reshape((1, 1, -1)) + if D_r_omega.ndim == 1: + D_r_omega = D_r_omega.reshape((1, 1, -1)) tf_fit, _ = dkpy.DScaleFitSlicot().fit( - omega, D_omega, order, block_structure + omega, D_l_omega, D_r_omega, order, block_structure ) class TestGenerateDScaleMask: @pytest.mark.parametrize( - "block_structure, mask_exp", + "block_structure, mask_l_exp, mask_r_exp", [ ( [dkpy.ComplexFullBlock(1, 1), dkpy.ComplexFullBlock(1, 1)], @@ -378,6 +369,13 @@ class TestGenerateDScaleMask: ], dtype=int, ), + np.array( + [ + [-1, 0], + [0, 1], + ], + dtype=int, + ), ), ( [dkpy.ComplexFullBlock(2, 2), dkpy.ComplexFullBlock(1, 1)], @@ -389,6 +387,14 @@ class TestGenerateDScaleMask: ], dtype=int, ), + np.array( + [ + [-1, 0, 0], + [0, -1, 0], + [0, 0, 1], + ], + dtype=int, + ), ), ( [dkpy.ComplexFullBlock(1, 1), dkpy.ComplexFullBlock(2, 2)], @@ -400,13 +406,98 @@ class TestGenerateDScaleMask: ], dtype=int, ), + np.array( + [ + [-1, 0, 0], + [0, 1, 0], + [0, 0, 1], + ], + dtype=int, + ), + ), + ( + [dkpy.ComplexDiagonalBlock(2), dkpy.ComplexFullBlock(2, 3)], + np.array( + [ + [-1, -1, 0, 0, 0], + [0, -1, 0, 0, 0], + [0, 0, 1, 0, 0], + [0, 0, 0, 1, 0], + [0, 0, 0, 0, 1], + ], + dtype=int, + ), + np.array( + [ + [-1, -1, 0, 0], + [0, -1, 0, 0], + [0, 0, 1, 0], + [0, 0, 0, 1], + ], + dtype=int, + ), + ), + ( + [dkpy.ComplexFullBlock(2, 3), dkpy.ComplexDiagonalBlock(2)], + np.array( + [ + [1, 0, 0, 0, 0], + [0, 1, 0, 0, 0], + [0, 0, 1, 0, 0], + [0, 0, 0, -1, -1], + [0, 0, 0, 0, -1], + ], + dtype=int, + ), + np.array( + [ + [1, 0, 0, 0], + [0, 1, 0, 0], + [0, 0, -1, -1], + [0, 0, 0, -1], + ], + dtype=int, + ), + ), + ( + [ + dkpy.ComplexDiagonalBlock(3), + dkpy.ComplexFullBlock(2, 3), + dkpy.ComplexDiagonalBlock(2), + ], + np.array( + [ + [-1, -1, -1, 0, 0, 0, 0, 0], + [0, -1, -1, 0, 0, 0, 0, 0], + [0, 0, -1, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 0, 0, -1, -1], + [0, 0, 0, 0, 0, 0, 0, -1], + ], + dtype=int, + ), + np.array( + [ + [-1, -1, -1, 0, 0, 0, 0], + [0, -1, -1, 0, 0, 0, 0], + [0, 0, -1, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 0, -1, -1], + [0, 0, 0, 0, 0, 0, -1], + ], + dtype=int, + ), ), ], ) - def test_generate_d_scale_mask(self, block_structure, mask_exp): + def test_generate_d_scale_mask(self, block_structure, mask_l_exp, mask_r_exp): """Test :func:`_mask_from_block_strucure`.""" - mask = dkpy.d_scale_fit._generate_d_scale_mask(block_structure) - np.testing.assert_allclose(mask_exp, mask) + mask_l, mask_r = dkpy.d_scale_fit._generate_d_scale_mask(block_structure) + np.testing.assert_allclose(mask_l_exp, mask_l) + np.testing.assert_allclose(mask_r_exp, mask_r) class TestInvertBiproperSs: diff --git a/tests/test_structured_singular_value.py b/tests/test_structured_singular_value.py index 6d34c2e..58a88b1 100644 --- a/tests/test_structured_singular_value.py +++ b/tests/test_structured_singular_value.py @@ -13,17 +13,31 @@ class TestGenerateSsvVariable: [ ( [dkpy.ComplexFullBlock(1, 1), dkpy.ComplexFullBlock(2, 2)], - cvxpy.bmat( - [ + ( + cvxpy.bmat( [ - cvxpy.Variable((1, 1), complex=True, name="x0"), - np.zeros((1, 2)), - ], + [ + cvxpy.Variable((), complex=True, name="x0") * np.eye(1), + np.zeros((1, 2)), + ], + [ + np.zeros((2, 1)), + np.eye(2), + ], + ] + ), + cvxpy.bmat( [ - np.zeros((2, 1)), - np.eye(2), - ], - ] + [ + cvxpy.Variable((), complex=True, name="x0") * np.eye(1), + np.zeros((1, 2)), + ], + [ + np.zeros((2, 1)), + np.eye(2), + ], + ] + ), ), ), ( @@ -32,24 +46,200 @@ class TestGenerateSsvVariable: dkpy.ComplexFullBlock(2, 2), dkpy.ComplexFullBlock(1, 1), ], - cvxpy.bmat( - [ + ( + cvxpy.bmat( [ - cvxpy.Variable((1, 1), complex=True, name="x0"), - np.zeros((1, 2)), - np.zeros((1, 1)), - ], + [ + cvxpy.Variable((), complex=True, name="x0") * np.eye(1), + np.zeros((1, 2)), + np.zeros((1, 1)), + ], + [ + np.zeros((2, 1)), + cvxpy.Variable((), complex=True, name="x1") * np.eye(2), + np.zeros((2, 1)), + ], + [ + np.zeros((1, 1)), + np.zeros((1, 2)), + np.eye(1), + ], + ] + ), + cvxpy.bmat( [ - np.zeros((2, 1)), - cvxpy.Variable(1, complex=True, name="x1") * np.eye(2), - np.zeros((2, 1)), + [ + cvxpy.Variable((), complex=True, name="x0") * np.eye(1), + np.zeros((1, 2)), + np.zeros((1, 1)), + ], + [ + np.zeros((2, 1)), + cvxpy.Variable((), complex=True, name="x1") * np.eye(2), + np.zeros((2, 1)), + ], + [ + np.zeros((1, 1)), + np.zeros((1, 2)), + np.eye(1), + ], + ] + ), + ), + ), + ( + [ + dkpy.ComplexDiagonalBlock(2), + dkpy.ComplexDiagonalBlock(1), + dkpy.ComplexFullBlock(2, 3), + dkpy.ComplexFullBlock(3, 3), + dkpy.ComplexDiagonalBlock(3), + dkpy.ComplexFullBlock(2, 1), + ], + ( + cvxpy.bmat( + [ + [ + cvxpy.Variable((2, 2), hermitian=True, name="X0"), + np.zeros((2, 1)), + np.zeros((2, 3)), + np.zeros((2, 3)), + np.zeros((2, 3)), + np.zeros((2, 1)), + ], + [ + np.zeros((1, 2)), + cvxpy.Variable((1, 1), hermitian=True, name="X1"), + np.zeros((1, 3)), + np.zeros((1, 3)), + np.zeros((1, 3)), + np.zeros((1, 1)), + ], + [ + np.zeros((3, 2)), + np.zeros((3, 1)), + cvxpy.Variable((), complex=True, name="x2") * np.eye(3), + np.zeros((3, 3)), + np.zeros((3, 3)), + np.zeros((3, 1)), + ], + [ + np.zeros((3, 2)), + np.zeros((3, 1)), + np.zeros((3, 3)), + cvxpy.Variable((), complex=True, name="x3") * np.eye(3), + np.zeros((3, 3)), + np.zeros((3, 1)), + ], + [ + np.zeros((3, 2)), + np.zeros((3, 1)), + np.zeros((3, 3)), + np.zeros((3, 3)), + cvxpy.Variable((3, 3), hermitian=True, name="X4"), + np.zeros((3, 1)), + ], + [ + np.zeros((1, 2)), + np.zeros((1, 1)), + np.zeros((1, 3)), + np.zeros((1, 3)), + np.zeros((1, 3)), + np.eye(1), + ], ], + ), + cvxpy.bmat( [ - np.zeros((1, 1)), - np.zeros((1, 2)), - np.eye(1), + [ + cvxpy.Variable((2, 2), hermitian=True, name="X0"), + np.zeros((2, 1)), + np.zeros((2, 2)), + np.zeros((2, 3)), + np.zeros((2, 3)), + np.zeros((2, 2)), + ], + [ + np.zeros((1, 2)), + cvxpy.Variable((1, 1), hermitian=True, name="X1"), + np.zeros((1, 2)), + np.zeros((1, 3)), + np.zeros((1, 3)), + np.zeros((1, 2)), + ], + [ + np.zeros((2, 2)), + np.zeros((2, 1)), + cvxpy.Variable((), name="x2") * np.eye(2), + np.zeros((2, 3)), + np.zeros((2, 3)), + np.zeros((2, 2)), + ], + [ + np.zeros((3, 2)), + np.zeros((3, 1)), + np.zeros((3, 2)), + cvxpy.Variable((), complex=True, name="x3") * np.eye(3), + np.zeros((3, 3)), + np.zeros((3, 2)), + ], + [ + np.zeros((3, 2)), + np.zeros((3, 1)), + np.zeros((3, 2)), + np.zeros((3, 3)), + cvxpy.Variable((3, 3), hermitian=True, name="X4"), + np.zeros((3, 2)), + ], + [ + np.zeros((2, 2)), + np.zeros((2, 1)), + np.zeros((2, 2)), + np.zeros((2, 3)), + np.zeros((2, 3)), + np.eye(2), + ], ], - ] + ), + ), + ), + ( + [ + dkpy.ComplexDiagonalBlock(2), + dkpy.ComplexFullBlock(3, 3), + dkpy.ComplexDiagonalBlock(3), + ], + ( + cvxpy.bmat( + [ + [ + cvxpy.Variable((2, 2), hermitian=True, name="X0"), + np.zeros((2, 3)), + np.zeros((2, 3)), + ], + [np.zeros((3, 2)), np.eye(3), np.zeros((3, 3))], + [ + np.zeros((3, 2)), + np.zeros((3, 3)), + cvxpy.Variable((3, 3), hermitian=True, name="X2"), + ], + ] + ), + cvxpy.bmat( + [ + [ + cvxpy.Variable((2, 2), hermitian=True, name="X0"), + np.zeros((2, 3)), + np.zeros((2, 3)), + ], + [np.zeros((3, 2)), np.eye(3), np.zeros((3, 3))], + [ + np.zeros((3, 2)), + np.zeros((3, 3)), + cvxpy.Variable((3, 3), hermitian=True, name="X2"), + ], + ] + ), ), ), ], @@ -59,7 +249,64 @@ def test_generate_ssv_variable(self, block_structure, variable_exp): variable = dkpy.structured_singular_value._generate_ssv_variable( block_structure ) - assert variable.ndim == variable_exp.ndim - assert variable.shape == variable_exp.shape - # This is the only way I can think of to compare expressions - assert variable.name() == variable_exp.name() + assert variable[0].ndim == variable_exp[0].ndim + assert variable[1].ndim == variable_exp[1].ndim + assert variable[0].shape == variable_exp[0].shape + assert variable[1].shape == variable_exp[1].shape + assert variable[0].name() == variable_exp[0].name() + assert variable[1].name() == variable_exp[1].name() + + +class TestSsvLmiBisection: + @pytest.mark.parametrize( + "N_omega, block_structure, mu_exp", + [ + ( + np.array( + [ + [1j, 2, 2j, 0, -1, -1 + 3j, 2 + 3j], + [3 + 1j, 2 - 1j, -1 + 1j, 2 + 1j, -1 + 1j, 1, -1 + 1j], + [3 + 1j, 1j, 2 + 2j, -1 + 2j, 3 - 1j, 3j, -1 + 1j], + [-1 + 1j, -1 - 1j, 1j, 0, 1 - 1j, 2 - 1j, 2 + 2j], + [3, 1j, 1 + 1j, 3j, 1 + 1j, 3j, -1j], + [1, 3 + 2j, 2 + 2j, 3j, 1 + 2j, 2 + 1j, -1 + 2j], + [2 + 1j, -1 - 1j, -1, 3 + 3j, 2 + 3j, 2j, 1 - 1j], + ] + ), + [ + dkpy.ComplexDiagonalBlock(2), + dkpy.ComplexDiagonalBlock(1), + dkpy.ComplexFullBlock(2, 3), + dkpy.ComplexFullBlock(2, 1), + ], + 10.5523, + ), + ( + np.array( + [ + [1j, 2, 2j, 0], + [3 + 1j, 2 - 1j, -1 + 1j, 2 + 1j], + [3 + 1j, 1j, 2 + 2j, -1 + 2j], + [3 - 1j, 3j, -1 + 1j, 2 + 2j], + [-1 + 1j, -1 - 1j, 1j, 0], + [1 - 1j, 2 - 1j, 2 + 2j, 1], + [3, 1j, 1 + 1j, 3j], + [1 + 2j, 2 + 1j, -1 + 2j, 1j], + [2 + 1j, -1 - 1j, -1, 3 + 3j], + [2 + 3j, 2j, 1 - 1j, 0], + ] + ), + [ + dkpy.ComplexFullBlock(2, 4), + dkpy.ComplexFullBlock(2, 6), + ], + 10.2438, + ), + ], + ) + def test_compute_ssv(self, N_omega, block_structure, mu_exp): + N_omega = np.expand_dims(N_omega, N_omega.ndim) + mu, D_l_scales, D_r_scales, info = dkpy.SsvLmiBisection(n_jobs=1).compute_ssv( + N_omega, block_structure + ) + assert np.isclose(mu, mu_exp) diff --git a/tests/test_uncertainty_structure.py b/tests/test_uncertainty_structure.py index 45b576a..e590f02 100644 --- a/tests/test_uncertainty_structure.py +++ b/tests/test_uncertainty_structure.py @@ -1,7 +1,6 @@ """Test :mod:`uncertainty_structure`.""" import numpy as np -import cvxpy import pytest import dkpy @@ -14,8 +13,8 @@ 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.n_exogenous_inputs == 5 + assert block.n_exogenous_outputs == 5 assert block.is_diagonal assert not block.is_complex @@ -26,8 +25,8 @@ class TestComplexDiagonalBlock: 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.n_exogenous_inputs == 5 + assert block.n_exogenous_outputs == 5 assert block.is_diagonal assert block.is_complex @@ -38,8 +37,8 @@ class TestComplexFullBlock: 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 block.n_exogenous_inputs == 5 + assert block.n_exogenous_outputs == 10 assert not block.is_diagonal assert block.is_complex @@ -48,7 +47,7 @@ class TestUncertaintyBlockStructure: """Test :class:`UncertaintyBlockStructure""" @pytest.mark.parametrize( - "block_structure, num_inputs_list, num_outputs_list," + "block_structure, n_exogenous_inputs_list, n_exogenous_outputs_list," " is_diagonal_list, is_complex_list", [ ( @@ -60,15 +59,15 @@ class TestUncertaintyBlockStructure: ), ( np.array([[-3, 0], [6, 0], [1, 2]]), - [3, 6, 2], [3, 6, 1], + [3, 6, 2], [True, True, False], [False, True, True], ), ( np.array([[3, 6], [-1, 0], [5, 1], [10, 0]]), - [6, 1, 1, 10], [3, 1, 5, 10], + [6, 1, 1, 10], [False, True, False, True], [True, False, True, True], ), @@ -81,15 +80,15 @@ class TestUncertaintyBlockStructure: ), ( [[-3, 0], [6, 0], [1, 2]], - [3, 6, 2], [3, 6, 1], + [3, 6, 2], [True, True, False], [False, True, True], ), ( [[3, 6], [-1, 0], [5, 1], [10, 0]], - [6, 1, 1, 10], [3, 1, 5, 10], + [6, 1, 1, 10], [False, True, False, True], [True, False, True, True], ), @@ -108,22 +107,22 @@ class TestUncertaintyBlockStructure: [ dkpy.RealDiagonalBlock(3), dkpy.ComplexDiagonalBlock(6), - dkpy.ComplexFullBlock(2, 1), + dkpy.ComplexFullBlock(1, 2), ], - [3, 6, 2], [3, 6, 1], + [3, 6, 2], [True, True, False], [False, True, True], ), ( [ - dkpy.ComplexFullBlock(6, 3), + dkpy.ComplexFullBlock(3, 6), dkpy.RealDiagonalBlock(1), - dkpy.ComplexFullBlock(1, 5), + dkpy.ComplexFullBlock(5, 1), dkpy.ComplexDiagonalBlock(10), ], - [6, 1, 1, 10], [3, 1, 5, 10], + [6, 1, 1, 10], [False, True, False, True], [True, False, True, True], ), @@ -132,8 +131,8 @@ class TestUncertaintyBlockStructure: def test_convert_block_structure_representation( self, block_structure, - num_inputs_list, - num_outputs_list, + n_exogenous_inputs_list, + n_exogenous_outputs_list, is_diagonal_list, is_complex_list, ): @@ -145,7 +144,7 @@ def test_convert_block_structure_representation( 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.n_exogenous_inputs == n_exogenous_inputs_list[idx_block] + assert block.n_exogenous_outputs == n_exogenous_outputs_list[idx_block] assert block.is_diagonal == is_diagonal_list[idx_block] assert block.is_complex == is_complex_list[idx_block]