Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Improve sys-mor demos #1198

Merged
merged 6 commits into from
Dec 3, 2020
Merged

Improve sys-mor demos #1198

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
4345d92
[demos] improve non-parametric sys-mor demos
pmli Dec 1, 2020
416ee5d
[demos] improve parametric sys-mor demos
pmli Dec 1, 2020
1a3092c
[tests] reduce dimensions of sys-mor demos
pmli Dec 1, 2020
2c28a08
[tests] fix some flake8 issues
pmli Dec 1, 2020
3f575cc
[demos] replace Options with Arguments
pmli Dec 2, 2020
3ed84b3
[tests] run sys-mor demos also without slycot
pmli Dec 2, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
64 changes: 8 additions & 56 deletions src/pymordemos/delay.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)]])

Expand All @@ -30,9 +29,11 @@ def dH(s):

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, '.')
Expand All @@ -41,6 +42,7 @@ def dH(s):
ax.set_ylabel('Im')
plt.show()

# Magnitude plots
w = np.logspace(-1, 3, 200)

fig, ax = plt.subplots()
Expand All @@ -54,56 +56,6 @@ def dH(s):
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)
132 changes: 70 additions & 62 deletions src/pymordemos/heat.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:

Expand All @@ -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(
Expand Down Expand Up @@ -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__':
Expand Down
Loading