Unverified Commit aa3011ac authored by Petar Mlinarić's avatar Petar Mlinarić Committed by GitHub
Browse files

Merge pull request #1198 from pymor/sys-mor-demos

Improve sys-mor demos
parents 6c81dc07 3ed84b36
Pipeline #69712 passed with stages
in 33 minutes and 18 seconds
......@@ -6,22 +6,21 @@
import matplotlib.pyplot as plt
import numpy as np
import scipy.linalg as spla
from typer import Option, run
from typer import Argument, run
from pymor.models.iosys import TransferFunction
from pymor.reductors.interpolation import TFBHIReductor
from pymor.reductors.h2 import TFIRKAReductor
def main(
tau: float = Option(0.1, help='The time delay.'),
r: int = Option(10, help='Order of the TF-IRKA ROM.'),
tau: float = Argument(0.1, help='The time delay.'),
r: int = Argument(10, help='Order of the TF-IRKA ROM.'),
):
"""Delay demo
"""Delay demo.
Cascade of delay and integrator
"""
# Transfer function
def H(s):
return np.array([[np.exp(-s) / (tau * s + 1)]])
......@@ -30,9 +29,11 @@ def main(
tf = TransferFunction(1, 1, H, dH)
# Transfer function IRKA (TF-IRKA)
tf_irka_reductor = TFIRKAReductor(tf)
rom = tf_irka_reductor.reduce(r, maxit=1000)
# Final interpolation points
sigma_list = tf_irka_reductor.sigma_list
fig, ax = plt.subplots()
ax.plot(sigma_list[-1].real, sigma_list[-1].imag, '.')
......@@ -41,6 +42,7 @@ def main(
ax.set_ylabel('Im')
plt.show()
# Magnitude plots
w = np.logspace(-1, 3, 200)
fig, ax = plt.subplots()
......@@ -54,56 +56,6 @@ def main(
ax.set_title('Magnitude plots of the error system')
plt.show()
# step response
E = rom.E.matrix
A = rom.A.matrix
B = rom.B.matrix
C = rom.C.matrix
nt = 1000
t = np.linspace(0, 4, nt)
x_old = np.zeros(rom.order)
y = np.zeros(nt)
for i in range(1, nt):
h = t[i] - t[i - 1]
x_new = spla.solve(E - h / 2 * A, (E + h / 2 * A).dot(x_old) + h * B[:, 0])
x_old = x_new
y[i] = C.dot(x_new)[0]
step_response = np.piecewise(t, [t < 1, t >= 1], [0, 1]) * (1 - np.exp(-(t - 1) / tau))
fig, ax = plt.subplots()
ax.plot(t, step_response, '-', t, y, '--')
ax.set_title('Step responses of the full and reduced model')
ax.set_xlabel('$t$')
plt.show()
# match steady state (add interpolation point at 0)
sigma_ss = list(sigma_list[-1]) + [0]
b_ss = np.ones((r+1, tf.dim_input))
c_ss = np.ones((r+1, tf.dim_output))
interp_reductor = TFBHIReductor(tf)
rom_ss = interp_reductor.reduce(sigma_ss, b_ss, c_ss)
# step response
E_ss = rom_ss.E.matrix
A_ss = rom_ss.A.matrix
B_ss = rom_ss.B.matrix
C_ss = rom_ss.C.matrix
x_ss_old = np.zeros(rom_ss.order)
y_ss = np.zeros(nt)
for i in range(1, nt):
h = t[i] - t[i - 1]
x_ss_new = spla.solve(E_ss - h / 2 * A_ss, (E_ss + h / 2 * A_ss).dot(x_ss_old) + h * B_ss[:, 0])
x_ss_old = x_ss_new
y_ss[i] = C_ss.dot(x_ss_new)[0]
fig, ax = plt.subplots()
ax.plot(t, step_response, '-', t, y_ss, '--')
ax.set_title('Step responses of the full and reduced model 2')
ax.set_xlabel('$t$')
plt.show()
if __name__ == '__main__':
run(main)
......@@ -5,21 +5,75 @@
import numpy as np
import matplotlib.pyplot as plt
from typer import Option, run
from typer import Argument, run
from pymor.basic import (InstationaryProblem, StationaryProblem, RectDomain, ConstantFunction, ExpressionFunction,
discretize_instationary_cg, BTReductor, IRKAReductor)
from pymor.analyticalproblems.domaindescriptions import RectDomain
from pymor.analyticalproblems.elliptic import StationaryProblem
from pymor.analyticalproblems.functions import ConstantFunction, ExpressionFunction
from pymor.analyticalproblems.instationary import InstationaryProblem
from pymor.core.config import config
from pymor.core.logger import set_log_levels
from pymor.discretizers.builtin import discretize_instationary_cg
from pymor.reductors.bt import BTReductor, LQGBTReductor, BRBTReductor
from pymor.reductors.h2 import IRKAReductor, TSIAReductor, OneSidedIRKAReductor
def run_mor_method(lti, w, reductor, reductor_short_name, r, **reduce_kwargs):
"""Run a model order reduction method.
Parameters
----------
lti
The full-order |LTIModel|.
w
Array of frequencies.
reductor
The reductor object.
reductor_short_name
A short name for the reductor.
r
The order of the reduced-order model.
reduce_kwargs
Optional keyword arguments for the reduce method.
"""
# Reduction
rom = reductor.reduce(r, **reduce_kwargs)
err = lti - rom
# Poles of the reduced-order model
poles_rom = rom.poles()
fig, ax = plt.subplots()
ax.plot(poles_rom.real, poles_rom.imag, '.')
ax.set_title(f"{reductor_short_name} reduced model's poles")
plt.show()
# Errors
print(f'{reductor_short_name} relative H_2-error: {err.h2_norm() / lti.h2_norm():e}')
if config.HAVE_SLYCOT:
print(f'{reductor_short_name} relative H_inf-error: {err.hinf_norm() / lti.hinf_norm():e}')
else:
print('Skipped H_inf-norm calculation due to missing slycot.')
print(f'{reductor_short_name} relative Hankel-error: {err.hankel_norm() / lti.hankel_norm():e}')
# Magnitude plot of the full and reduced model
fig, ax = plt.subplots()
lti.mag_plot(w, ax=ax)
rom.mag_plot(w, ax=ax, linestyle='dashed')
ax.set_title(f'Magnitude plot of the full and {reductor_short_name} reduced model')
plt.show()
import logging
logging.getLogger('pymor.algorithms.gram_schmidt.gram_schmidt').setLevel(logging.ERROR)
# Magnitude plot of the error system
fig, ax = plt.subplots()
err.mag_plot(w, ax=ax)
ax.set_title(f'Magnitude plot of the {reductor_short_name} error system')
plt.show()
def main(
diameter: float = Option(0.1, help='Diameter option for the domain discretizer.'),
r: int = Option(5, help='Order of the ROMs.'),
diameter: float = Argument(0.1, help='Diameter option for the domain discretizer.'),
r: int = Argument(5, help='Order of the ROMs.'),
):
r"""2D heat equation demo
r"""2D heat equation demo.
Discretization of the PDE:
......@@ -38,6 +92,7 @@ def main(
where :math:`u(t)` is the input and :math:`y(t)` is the output.
"""
set_log_levels({'pymor.algorithms.gram_schmidt.gram_schmidt': 'WARNING'})
p = InstationaryProblem(
StationaryProblem(
......@@ -89,60 +144,13 @@ def main(
print('Skipped H_inf-norm calculation due to missing slycot.')
print(f'FOM Hankel-norm: {lti.hankel_norm():e}')
# Balanced Truncation
reductor = BTReductor(lti)
rom_bt = reductor.reduce(r, tol=1e-5)
err_bt = lti - rom_bt
print(f'BT relative H_2-error: {err_bt.h2_norm() / lti.h2_norm():e}')
if config.HAVE_SLYCOT:
print(f'BT relative H_inf-error: {err_bt.hinf_norm() / lti.hinf_norm():e}')
else:
print('Skipped H_inf-norm calculation due to missing slycot.')
print(f'BT relative Hankel-error: {err_bt.hankel_norm() / lti.hankel_norm():e}')
# Magnitude plot of the full and BT reduced model
fig, ax = plt.subplots()
lti.mag_plot(w, ax=ax)
rom_bt.mag_plot(w, ax=ax, linestyle='dashed')
ax.set_title('Magnitude plot of the full and BT reduced model')
plt.show()
# Magnitude plot of the BT error system
fig, ax = plt.subplots()
err_bt.mag_plot(w, ax=ax)
ax.set_title('Magnitude plot of the BT error system')
plt.show()
# Iterative Rational Krylov Algorithm
irka_reductor = IRKAReductor(lti)
rom_irka = irka_reductor.reduce(r, compute_errors=True)
# Shift distances
fig, ax = plt.subplots()
ax.semilogy(irka_reductor.conv_crit, '.-')
ax.set_title('Distances between shifts in IRKA iterations')
plt.show()
err_irka = lti - rom_irka
print(f'IRKA relative H_2-error: {err_irka.h2_norm() / lti.h2_norm():e}')
if config.HAVE_SLYCOT:
print(f'IRKA relative H_inf-error: {err_irka.hinf_norm() / lti.hinf_norm():e}')
else:
print('Skipped H_inf-norm calculation due to missing slycot.')
print(f'IRKA relative Hankel-error: {err_irka.hankel_norm() / lti.hankel_norm():e}')
# Magnitude plot of the full and IRKA reduced model
fig, ax = plt.subplots()
lti.mag_plot(w, ax=ax)
rom_irka.mag_plot(w, ax=ax, linestyle='dashed')
ax.set_title('Magnitude plot of the full and IRKA reduced model')
plt.show()
# Magnitude plot of the IRKA error system
fig, ax = plt.subplots()
err_irka.mag_plot(w, ax=ax)
ax.set_title('Magnitude plot of the IRKA error system')
plt.show()
# Model order reduction
run_mor_method(lti, w, BTReductor(lti), 'BT', r, tol=1e-5)
run_mor_method(lti, w, LQGBTReductor(lti), 'LQGBT', r, tol=1e-5)
run_mor_method(lti, w, BRBTReductor(lti), 'BRBT', r, tol=1e-5)
run_mor_method(lti, w, IRKAReductor(lti), 'IRKA', r)
run_mor_method(lti, w, TSIAReductor(lti), 'TSIA', r)
run_mor_method(lti, w, OneSidedIRKAReductor(lti, 'V'), 'OS-IRKA', r)
if __name__ == '__main__':
......
#!/usr/bin/env python
# coding: utf-8
# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2020 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
# In[ ]:
import numpy as np
import scipy.linalg as spla
import matplotlib as mpl
import matplotlib.pyplot as plt
from typer import Option, run
from pymor.basic import *
def main(r: int = Option(10, help='Order of the TF-IRKA ROM.')):
# # Model
import numpy as np
from typer import Argument, run
# In[ ]:
from pymor.models.iosys import TransferFunction
from pymor.reductors.h2 import TFIRKAReductor
def main(r: int = Argument(10, help='Order of the TF-IRKA ROM.')):
"""Parametric delay demo."""
# Model
def H(s, mu):
tau = mu['tau'][0]
return np.array([[np.exp(-s) / (tau * s + 1)]])
......@@ -29,101 +23,44 @@ def main(r: int = Option(10, help='Order of the TF-IRKA ROM.')):
tau = mu['tau'][0]
return np.array([[-(tau * s + tau + 1) * np.exp(-s) / (tau * s + 1) ** 2]])
# In[ ]:
fom = TransferFunction(1, 1,
H, dH,
parameters={'tau': 1})
# # Magnitude plot
# In[ ]:
mu_list_short = [0.01, 0.1, 1]
# In[ ]:
# Magnitude plot
mu_list = [0.01, 0.1, 1]
w = np.logspace(-2, 4, 100)
fig, ax = plt.subplots()
for mu in mu_list_short:
for mu in mu_list:
fom.mag_plot(w, ax=ax, mu=mu, label=fr'$\tau = {mu}$')
ax.legend()
plt.show()
# In[ ]:
w_list = np.logspace(-2, 4, 100)
mu_list = np.logspace(-2, 0, 50)
fom_w_mu = np.zeros((len(w_list), len(mu_list)))
for i, mu in enumerate(mu_list):
fom_w_mu[:, i] = spla.norm(fom.freq_resp(w_list, mu=mu), axis=(1, 2))
# In[ ]:
fig, ax = plt.subplots()
out = ax.contourf(w_list, mu_list, fom_w_mu.T,
norm=mpl.colors.LogNorm(),
levels=np.logspace(np.log10(fom_w_mu.min()), np.log10(fom_w_mu.max()), 100))
ax.set_xlabel(r'Frequency $\omega$')
ax.set_ylabel(r'Parameter $\mu$')
ax.set_xscale('log')
ax.set_yscale('log')
fig.colorbar(out, ticks=np.logspace(-4, 1, 6))
plt.show()
# # TF-IRKA
# In[ ]:
# TF-IRKA
roms_tf_irka = []
for mu in mu_list_short:
for mu in mu_list:
tf_irka = TFIRKAReductor(fom, mu=mu)
rom = tf_irka.reduce(r, conv_crit='h2', maxit=1000, num_prev=5)
roms_tf_irka.append(rom)
# In[ ]:
fig, ax = plt.subplots()
for mu, rom in zip(mu_list_short, roms_tf_irka):
for mu, rom in zip(mu_list, roms_tf_irka):
poles = rom.poles()
ax.plot(poles.real, poles.imag, '.', label=fr'$\tau = {mu}$')
ax.set_title("Poles of TF-IRKA's ROMs")
ax.legend()
plt.show()
# In[ ]:
fig, ax = plt.subplots()
for mu, rom in zip(mu_list_short, roms_tf_irka):
for mu, rom in zip(mu_list, roms_tf_irka):
rom.mag_plot(w, ax=ax, label=fr'$\tau = {mu}$')
ax.set_title("Magnitude plot of TF-IRKA's ROMs")
ax.legend()
plt.show()
# In[ ]:
fig, ax = plt.subplots()
for mu, rom in zip(mu_list_short, roms_tf_irka):
for mu, rom in zip(mu_list, roms_tf_irka):
(fom - rom).mag_plot(w, ax=ax, mu=mu, label=fr'$\tau = {mu}$')
ax.set_title("Magnitude plot of error systems")
ax.legend()
......
#!/usr/bin/env python
# coding: utf-8
# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2020 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
# In[ ]:
# In[ ]:
import numpy as np
import scipy.linalg as spla
import matplotlib.pyplot as plt
import matplotlib as mpl
from typer import Option, run
from typer import Argument, run
from pymor.basic import *
from pymor.analyticalproblems.domaindescriptions import LineDomain
from pymor.analyticalproblems.elliptic import StationaryProblem
from pymor.analyticalproblems.functions import ConstantFunction, ExpressionFunction, LincombFunction
from pymor.analyticalproblems.instationary import InstationaryProblem
from pymor.core.config import config
from pymor.core.logger import set_log_levels
from pymor.discretizers.builtin import discretize_instationary_cg
from pymor.parameters.functionals import ProjectionParameterFunctional
from pymor.reductors.bt import BTReductor, LQGBTReductor, BRBTReductor
from pymor.reductors.h2 import IRKAReductor, TSIAReductor, OneSidedIRKAReductor
def run_mor_method_param(fom, r, w, mus, reductor_cls, reductor_short_name, **reductor_kwargs):
"""Plot reductor errors for different parameter values.
Parameters
----------
fom
The full-order |LTIModel|.
r
The order of the reduced-order model.
w
Array of frequencies.
mus
An array of parameter values.
reductor_cls
The reductor class.
reductor_short_name
A short name for the reductor.
reductor_kwargs
Optional keyword arguments for the reductor class.
"""
# Reduction
roms = []
for mu in mus:
rom = reductor_cls(fom, mu=mu, **reductor_kwargs).reduce(r)
roms.append(rom)
# Poles
fig, ax = plt.subplots()
for rom in roms:
poles_rom = rom.poles()
ax.plot(poles_rom.real, poles_rom.imag, '.', label=fr'$\mu = {mu}$')
ax.set_title(f"{reductor_short_name} reduced model's poles")
plt.show()
# Magnitude plots
fig, ax = plt.subplots()
for mu, rom in zip(mus, roms):
rom.mag_plot(w, ax=ax, label=fr'$\mu = {mu}$')
ax.set_title('Magnitude plot of {reductor_short_name} reduced models')
ax.legend()
plt.show()
def main(
diameter: float = Option(0.01, help='Diameter option for the domain discretizer.'),
r: int = Option(5, help='Order of the ROMs.'),
):
set_log_levels({'pymor.algorithms.gram_schmidt.gram_schmidt': 'WARNING'})
set_defaults({'pymor.discretizers.builtin.gui.jupyter.get_visualizer.backend': 'not pythreejs'})
fig, ax = plt.subplots()
for mu, rom in zip(mus, roms):
(fom - rom).mag_plot(w, ax=ax, mu=mu, label=fr'$\mu = {mu}$')
ax.set_title('Magnitude plot of the {reductor_short_name} error system')
ax.legend()
plt.show()
# # Model
# Errors
for mu, rom in zip(mus, roms):
err = fom - rom
print(f'mu = {mu}')
print(f' {reductor_short_name} relative H_2-error:'
f' {err.h2_norm(mu=mu) / fom.h2_norm(mu=mu):e}')
if config.HAVE_SLYCOT:
print(f' {reductor_short_name} relative H_inf-error:'
f' {err.hinf_norm(mu=mu) / fom.hinf_norm(mu=mu):e}')
print(f' {reductor_short_name} relative Hankel-error:'
f' {err.hankel_norm(mu=mu) / fom.hankel_norm(mu=mu):e}')
# In[ ]:
def main(
diameter: float = Argument(0.01, help='Diameter option for the domain discretizer.'),
r: int = Argument(5, help='Order of the ROMs.'),
):
"""Parametric 1D heat equation example."""
set_log_levels({'pymor.algorithms.gram_schmidt.gram_schmidt': 'WARNING'})
# Model
p = InstationaryProblem(
StationaryProblem(
domain=LineDomain([0.,1.], left='robin', right='robin'),
domain=LineDomain([0., 1.], left='robin', right='robin'),
diffusion=LincombFunction([ExpressionFunction('(x[...,0] <= 0.5) * 1.', 1),
ExpressionFunction('(0.5 < x[...,0]) * 1.', 1)],
[1,
......@@ -53,101 +105,37 @@ def main(
fom, _ = discretize_instationary_cg(p, diameter=diameter, nt=100)
# In[ ]:
fom.visualize(fom.solve(mu=0.1))
# In[ ]:
fom.visualize(fom.solve(mu=1))
# In[ ]:
fom.visualize(fom.solve(mu=10))
# In[ ]:
lti = fom.to_lti()
# # System analysis
# In[ ]:
print(f'order of the model = {lti.order}')
print(f'number of inputs = {lti.dim_input}')
print(f'number of outputs = {lti.dim_output}')
# In[ ]:
mu_list = [0.1, 1, 10]
w_list = np.logspace(-1, 3, 100)
fig, ax = plt.subplots(len(mu_list), 1, sharex=True, sharey=True)
for i, mu in enumerate(mu_list):
# System poles
fig, ax = plt.subplots()
for mu in mu_list:
poles = lti.poles(mu=mu)
ax[i].plot(poles.real, poles.imag, '.')
ax[i].set_xscale('symlog')
ax[i].set_title(fr'$\mu = {mu}$')
fig.suptitle('System poles')
fig.subplots_adjust(hspace=0.5)
ax.plot(poles.real, poles.imag, '.', label=fr'$\mu = {mu}$')
ax.set_title('System poles')
ax.legend()
plt.show()
# In[ ]:
mu_list = [0.1, 1, 10]
# Magnitude plots
fig, ax = plt.subplots()
w = np.logspace(-1, 3, 100)
for mu in mu_list:
lti.mag_plot(w, ax=ax, mu=mu, label=fr'$\mu = {mu}$')
lti.mag_plot(w_list, ax=ax, mu=mu, label=fr'$\mu = {mu}$')
ax.set_title('Magnitude plot of the full model')
ax.legend()
plt.show()
# In[ ]: