From 584e20112a3d92110793c46e074325a178b88d42 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Petar=20Mlinari=C4=87?= Date: Fri, 3 Nov 2023 11:54:21 -0400 Subject: [PATCH] Add scripts --- runme_cd_r6.py | 70 ++++++++++++++++++++++++++++++++++++++++++ runme_gab3_r1_c1.py | 66 ++++++++++++++++++++++++++++++++++++++++ runme_gab3_r2_c0.py | 69 ++++++++++++++++++++++++++++++++++++++++++ runme_gab3_r2_c2.py | 74 +++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 279 insertions(+) create mode 100644 runme_cd_r6.py create mode 100644 runme_gab3_r1_c1.py create mode 100644 runme_gab3_r2_c0.py create mode 100644 runme_gab3_r2_c2.py diff --git a/runme_cd_r6.py b/runme_cd_r6.py new file mode 100644 index 0000000..ad46ccc --- /dev/null +++ b/runme_cd_r6.py @@ -0,0 +1,70 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.15.2 +# kernelspec: +# display_name: Python 3 (ipykernel) +# language: python +# name: python3 +# --- + +# %% +import numpy as np +from pymor.models.iosys import LTIModel + +from models import slicot_model +from reductors import IRKACauchyIndexReductor, IRKAWithLineSearchReductor +from tools import plot_combined, plot_fom, plot_rom, settings + +settings() + +# %% [markdown] +# # FOM + +# %% +fom = slicot_model('CDplayer') + +# %% +print(fom) + +# %% +w = (1e-2, 1e4) +_ = plot_fom(fom, w) + +# %% [markdown] +# # IRKA + +# %% +r = 6 +irka = IRKACauchyIndexReductor(fom) +A0 = np.diag(np.arange(-1, -r - 1, -1)) +B0 = np.ones((r, fom.dim_input)) +C0 = np.ones((fom.dim_output, r)) +rom0 = LTIModel.from_matrices(A0, B0, C0) +rom = irka.reduce(rom0, compute_errors=True, conv_crit='h2') + +# %% +_ = plot_rom(fom, rom, w, irka, 'IRKA') + +# %% +print((fom - rom).h2_norm() / fom.h2_norm()) + +# %% [markdown] +# # IRKA with line search + +# %% +irka_rgd = IRKAWithLineSearchReductor(fom) +rom_rgd = irka_rgd.reduce(rom0, compute_errors=True, conv_crit='h2') + +# %% +_ = plot_rom(fom, rom_rgd, w, irka_rgd, 'RGD-IRKA') + +# %% +print((fom - rom_rgd).h2_norm() / fom.h2_norm()) + +# %% +_ = plot_combined(irka, irka_rgd, 'cd_r6.pdf', h2_error_scale='log') diff --git a/runme_gab3_r1_c1.py b/runme_gab3_r1_c1.py new file mode 100644 index 0000000..df7986d --- /dev/null +++ b/runme_gab3_r1_c1.py @@ -0,0 +1,66 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.15.2 +# kernelspec: +# display_name: Python 3 (ipykernel) +# language: python +# name: python3 +# --- + +# %% +import numpy as np +from pymor.models.iosys import LTIModel + +from models import gab08 +from reductors import IRKACauchyIndexReductor, IRKAWithLineSearchReductor +from tools import plot_combined, plot_fom, plot_rom, settings + +settings() + +# %% [markdown] +# # FOM + +# %% +fom = gab08() + +# %% +print(fom) + +# %% +w = (1e-2, 1e2) +_ = plot_fom(fom, w) + +# %% [markdown] +# # IRKA + +# %% +irka = IRKACauchyIndexReductor(fom) +rom0 = LTIModel.from_matrices(np.array([[-0.27]]), np.eye(1), np.eye(1)) +rom = irka.reduce(rom0, compute_errors=True, conv_crit='h2') + +# %% +_ = plot_rom(fom, rom, w, irka, 'IRKA') + +# %% +print((fom - rom).h2_norm() / fom.h2_norm()) + +# %% [markdown] +# # IRKA with line search + +# %% +irka_rgd = IRKAWithLineSearchReductor(fom) +rom_rgd = irka_rgd.reduce(rom0, compute_errors=True, conv_crit='h2') + +# %% +_ = plot_rom(fom, rom_rgd, w, irka_rgd, 'RGD-IRKA') + +# %% +print((fom - rom_rgd).h2_norm() / fom.h2_norm()) + +# %% +_ = plot_combined(irka, irka_rgd, 'gab3_r1_c1.pdf') diff --git a/runme_gab3_r2_c0.py b/runme_gab3_r2_c0.py new file mode 100644 index 0000000..7d3a569 --- /dev/null +++ b/runme_gab3_r2_c0.py @@ -0,0 +1,69 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.15.2 +# kernelspec: +# display_name: Python 3 (ipykernel) +# language: python +# name: python3 +# --- + +# %% +import numpy as np +from pymor.models.iosys import LTIModel + +from models import gab08 +from reductors import IRKACauchyIndexReductor, IRKAWithLineSearchReductor +from tools import plot_combined, plot_fom, plot_rom, settings + +settings() + +# %% [markdown] +# # FOM + +# %% +fom = gab08() + +# %% +print(fom) + +# %% +w = (1e-2, 1e2) +_ = plot_fom(fom, w) + +# %% [markdown] +# # IRKA + +# %% +irka = IRKACauchyIndexReductor(fom) +A0 = np.array([[-1, 1], [-1, -1]]) +B0 = np.ones((2, 1)) +C0 = np.ones((1, 2)) +rom0 = LTIModel.from_matrices(A0, B0, C0) +rom = irka.reduce(rom0, compute_errors=True, conv_crit='h2') + +# %% +_ = plot_rom(fom, rom, w, irka, 'IRKA') + +# %% +print((fom - rom).h2_norm() / fom.h2_norm()) + +# %% [markdown] +# # IRKA with line search + +# %% +irka_rgd = IRKAWithLineSearchReductor(fom) +rom_rgd = irka_rgd.reduce(rom0, compute_errors=True, conv_crit='h2') + +# %% +_ = plot_rom(fom, rom_rgd, w, irka_rgd, 'RGD-IRKA') + +# %% +print((fom - rom_rgd).h2_norm() / fom.h2_norm()) + +# %% +_ = plot_combined(irka, irka_rgd, 'gab3_r2_c0.pdf', show_ci=False) diff --git a/runme_gab3_r2_c2.py b/runme_gab3_r2_c2.py new file mode 100644 index 0000000..5e0fa80 --- /dev/null +++ b/runme_gab3_r2_c2.py @@ -0,0 +1,74 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.15.2 +# kernelspec: +# display_name: Python 3 (ipykernel) +# language: python +# name: python3 +# --- + +# %% +import numpy as np +from pymor.models.iosys import LTIModel, _lti_to_poles_b_c + +from models import gab08 +from reductors import IRKACauchyIndexReductor, IRKAWithLineSearchReductor +from tools import plot_combined, plot_fom, plot_rom, settings + +settings() + +# %% [markdown] +# # FOM + +# %% +fom = gab08() + +# %% +print(fom) + +# %% +w = (1e-2, 1e2) +_ = plot_fom(fom, w) + +# %% [markdown] +# # IRKA + +# %% +irka = IRKACauchyIndexReductor(fom) +A0 = np.array([[-1, 0], [0, -2]]) +B0 = np.ones((2, 1)) +C0 = np.ones((1, 2)) +rom0 = LTIModel.from_matrices(A0, B0, C0) +rom = irka.reduce(rom0, compute_errors=True, conv_crit='h2') + +# %% +_ = plot_rom(fom, rom, w, irka, 'IRKA') + +# %% +print((fom - rom).h2_norm() / fom.h2_norm()) + +# %% [markdown] +# # IRKA with line search + +# %% +irka_rgd = IRKAWithLineSearchReductor(fom) +rom_rgd = irka_rgd.reduce(rom0, compute_errors=True, conv_crit='h2') + +# %% +_ = plot_rom(fom, rom_rgd, w, irka_rgd, 'RGD-IRKA') + +# %% +print((fom - rom_rgd).h2_norm() / fom.h2_norm()) + +# %% +_ = plot_combined(irka, irka_rgd, 'gab3_r2_c2.pdf') + +# %% +p, b, c = _lti_to_poles_b_c(rom_rgd) +print('RGD-IRKA poles:', p) +print('RGD-IRKA residues:', [bi.dot(ci) for bi, ci in zip(b, c)])