diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..3dba2f5 --- /dev/null +++ b/.gitignore @@ -0,0 +1,436 @@ +### General ### +# Jetbrains +.idea + +# VS Code +.vscode + +# macOS +.DS_store + +### Python ### +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +pip-wheel-metadata/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +.python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + +### LaTeX ### +## Core latex/pdflatex auxiliary files: +*.aux +*.lof +*.log +*.lot +*.fls +*.out +*.toc +*.fmt +*.fot +*.cb +*.cb2 +.*.lb + +## Intermediate documents: +*.dvi +*.xdv +*-converted-to.* +# these rules might exclude image files for figures etc. +# *.ps +# *.eps +# *.pdf + +## Generated if empty string is given at "Please type another file name for output:" +MorNU21.pdf +MorNU22.pdf + +## Bibliography auxiliary files (bibtex/biblatex/biber): +*.bbl +*.bcf +*.blg +*-blx.aux +*-blx.bib +*.run.xml + +## Build tool auxiliary files: +*.fdb_latexmk +*.synctex +*.synctex(busy) +*.synctex.gz +*.synctex.gz(busy) +*.pdfsync + +## Build tool directories for auxiliary files +# latexrun +latex.out/ + +## Auxiliary and intermediate files from other packages: +# algorithms +*.alg +*.loa + +# achemso +acs-*.bib + +# amsthm +*.thm + +# beamer +*.nav +*.pre +*.snm +*.vrb + +# changes +*.soc + +# comment +*.cut + +# cprotect +*.cpt + +# elsarticle (documentclass of Elsevier journals) +*.spl + +# endnotes +*.ent + +*.lox + +# feynmf/feynmp +*.mf +*.mp +*.t[1-9] +*.t[1-9][0-9] +*.tfm + +#(r)(e)ledmac/(r)(e)ledpar +*.end +*.?end +*.[1-9] +*.[1-9][0-9] +*.[1-9][0-9][0-9] +*.[1-9]R +*.[1-9][0-9]R +*.[1-9][0-9][0-9]R +*.eledsec[1-9] +*.eledsec[1-9]R +*.eledsec[1-9][0-9] +*.eledsec[1-9][0-9]R +*.eledsec[1-9][0-9][0-9] +*.eledsec[1-9][0-9][0-9]R + +# glossaries +*.acn +*.acr +*.glg +*.glo +*.gls +*.glsdefs +*.lzo +*.lzs + +# uncomment this for glossaries-extra (will ignore makeindex's style files!) +# *.ist + +# gnuplottex +*-gnuplottex-* + +# gregoriotex +*.gaux +*.glog +*.gtex + +# htlatex +*.4ct +*.4tc +*.idv +*.lg +*.trc +*.xref + +# hyperref +*.brf + +# knitr +*-concordance.tex +# Uncomment the next line if you use knitr and want to ignore its generated tikz files +# *.tikz +*-tikzDictionary + +# listings +*.lol + +# luatexja-ruby +*.ltjruby + +# makeidx +*.idx +*.ilg +*.ind + +# minitoc +*.maf +*.mlf +*.mlt +*.mtc[0-9]* +*.slf[0-9]* +*.slt[0-9]* +*.stc[0-9]* + +# minted +_minted* +*.pyg + +# morewrites +*.mw + +# newpax +*.newpax + +# nomencl +*.nlg +*.nlo +*.nls + +# pax +*.pax + +# pdfpcnotes +*.pdfpc + +# sagetex +*.sagetex.sage +*.sagetex.py +*.sagetex.scmd + +# scrwfile +*.wrt + +# sympy +*.sout +*.sympy +sympy-plots-for-*.tex/ + +# pdfcomment +*.upa +*.upb + +# pythontex +*.pytxcode +pythontex-files-*/ + +# tcolorbox +*.listing + +# thmtools +*.loe + +# TikZ & PGF +*.dpth +*.md5 +*.auxlock + +# todonotes +*.tdo + +# vhistory +*.hst +*.ver + +*.lod + +# xcolor +*.xcp + +# xmpincl +*.xmpi + +# xindy +*.xdy + +# xypic precompiled matrices and outlines +*.xyc +*.xyd + +# endfloat +*.ttt +*.fff + +# Latexian +TSWLatexianTemp* + +## Editors: +# WinEdt +*.bak +*.sav + +# Texpad +.texpadtmp + +# LyX +*.lyx~ + +# Kile +*.backup + +# gummi +.*.swp + +# KBibTeX +*~[0-9]* + +# TeXnicCenter +*.tps + +# auto folder when using emacs and auctex +./auto/* +*.el + +# expex forward references with \gathertags +*-tags.tex + +# standalone packages +*.sta + +# Makeindex log files +*.lpz + +# xwatermark package +*.xwm + +# REVTeX puts footnotes in the bibliography by default, unless the nofootinbib +# option is specified. Footnotes are the stored in a file with suffix Notes.bib. +# Uncomment the next line to have this generated file ignored. +#*Notes.bib + +### LaTeX Patch ### +# LIPIcs / OASIcs +*.vtc + +# glossaries +*.glstex diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..147b9fb --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2021 Jonas Nicodemus + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.md b/README.md new file mode 100644 index 0000000..ecc8286 --- /dev/null +++ b/README.md @@ -0,0 +1,88 @@ + + + + + + +# Port-Hamiltonian Dynamic Mode Decomposition + +We present a novel physics-informed system identification method to construct a passive linear time-invariant system. In more detail, for a given quadratic energy functional, measurements of the input, state, and output of a system in the time domain, we find a realization that approximates the data well while guaranteeing that the energy functional satisfies a dissipation inequality. To this end, we use the framework of port-Hamiltonian (pH) systems and modify the dynamic mode decomposition to be feasible for continuous-time pH systems. We propose an iterative numerical method to solve the corresponding least-squares minimization problem. We construct an effective initialization of the algorithm by studying the least-squares problem in a weighted norm, for which we present the analytical minimum-norm solution. The efficiency of the proposed method is demonstrated with several numerical examples. + + +
+

Table of Contents

 +
    +
  1. + Citing +
    2. +
  2. + Installation +
    3. +
  3. Usage
    4. +
  4. License
    5. +
  5. Contact
    6. +
+
+ + + +## Installation +A python environment is required with at least **Python 3.10**. + +Install dependencies via `pip`: + ```sh + pip install -r requirements.txt + ``` + + +## Usage + +The executable script `main.py` executes the `pHDMD` algorithm for the current configuration, defined in `config.py`. Both files are located in `src`. + + +## Documentation + +Documentation is available [online][docs-url] +or you can build it yourself from inside the `docs` directory +by executing: + + make html + +This will generate HTML documentation in `docs/build/html`. +It is required to have the `sphinx` dependencies installed. This can be done by + + pip install -r requirements.txt + +within the `docs` directory. + + + +## License + +Distributed under the MIT License. See `LICENSE` for more information. + + +## Contact +Jonas Nicodemus - jonas.nicodemus@simtech.uni-stuttgart.de + +Benjamin Unger - benjamin.unger@simtech.uni-stuttgart.de\ +Riccardo Morandin - morandin@math.tu-berlin.de + +Project Link: [https://github.com/Jonas-Nicodemus/phdmd][project-url] + +[license-shield]: https://img.shields.io/github/license/Jonas-Nicodemus/phdmd.svg?style=for-the-badge +[license-url]: https://github.com/Jonas-Nicodemus/phdmd/blob/main/LICENSE +[linkedin-shield]: https://img.shields.io/badge/-LinkedIn-black.svg?style=for-the-badge&logo=linkedin&colorB=555 +[linkedin-url]: https://linkedin.com/in/jonas-nicodemus-a34931209/ +[doi-shield]: https://img.shields.io/badge/DOI-10.5281%20%2F%20zenodo.5520662-blue.svg?style=for-the-badge +[doi-url]: https://zenodo.org/badge/latestdoi/409099116 +[arxiv-shield]: https://img.shields.io/badge/arXiv-2109.10793-b31b1b.svg?style=for-the-badge +[arxiv-url]: https://arxiv.org/abs/2109.10793 +[project-url]:https://github.com/Jonas-Nicodemus/phdmd +[docs-url]:https://jonas-nicodemus.github.io/phdmd/ \ No newline at end of file diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 0000000..d0c3cbf --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = source +BUILDDIR = build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/make.bat b/docs/make.bat new file mode 100644 index 0000000..dc1312a --- /dev/null +++ b/docs/make.bat @@ -0,0 +1,35 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=source +set BUILDDIR=build + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.https://www.sphinx-doc.org/ + exit /b 1 +) + +if "%1" == "" goto help + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/docs/requirements.txt b/docs/requirements.txt new file mode 100644 index 0000000..56e0564 --- /dev/null +++ b/docs/requirements.txt @@ -0,0 +1,4 @@ +Sphinx>=4.5.0 +sphinx-rtd-theme>=1.0.0 +sphinxcontrib-bibtex>=2.4.2 +sphinxcontrib-katex>=0.8.6 \ No newline at end of file diff --git a/docs/source/algorithm.rst b/docs/source/algorithm.rst new file mode 100644 index 0000000..be010af --- /dev/null +++ b/docs/source/algorithm.rst @@ -0,0 +1,37 @@ +algorithm package +================= + +Submodules +---------- + +algorithm.phdmd module +---------------------- + +.. automodule:: algorithm.phdmd + :members: + :undoc-members: + :show-inheritance: + +algorithm.skew\_procrustes module +--------------------------------- + +.. automodule:: algorithm.skew_procrustes + :members: + :undoc-members: + :show-inheritance: + +algorithm.spsd\_procrustes module +--------------------------------- + +.. automodule:: algorithm.spsd_procrustes + :members: + :undoc-members: + :show-inheritance: + +Module contents +--------------- + +.. automodule:: algorithm + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/bibliography.rst b/docs/source/bibliography.rst new file mode 100644 index 0000000..73b294f --- /dev/null +++ b/docs/source/bibliography.rst @@ -0,0 +1,8 @@ +************ +Bibliography +************ + +.. + see https://sphinxcontrib-bibtex.readthedocs.io/en/latest/usage.html for options + +.. bibliography:: diff --git a/docs/source/conf.py b/docs/source/conf.py new file mode 100644 index 0000000..1b4a7bc --- /dev/null +++ b/docs/source/conf.py @@ -0,0 +1,130 @@ +# Configuration file for the Sphinx documentation builder. +# +# This file only contains a selection of the most common options. For a full +# list see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +# -- Path setup -------------------------------------------------------------- + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +import os +import sys +sys.path.insert(0, os.path.abspath('../../src')) + + +# -- Project information ----------------------------------------------------- + +project = 'pHDMD' +copyright = '2022, Jonas Nicodemus' +author = 'Jonas Nicodemus' + +# The full version, including alpha/beta/rc tags +release = '1.0.0' + + +# -- General configuration --------------------------------------------------- + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = ['sphinx.ext.autodoc', + 'sphinx.ext.napoleon', + 'sphinxcontrib.katex', + 'sphinxcontrib.bibtex', +] + +import sphinxcontrib.katex as katex + +latex_macros = r""" + \def \i {\mathrm{i}} + \def \e #1{\mathrm{e}^{#1}} + \def \vec #1{\mathbf{#1}} + \def \x {\vec{x}} + \def \d {\operatorname{d}\!} + \def \dirac #1{\operatorname{\delta}\left(#1\right)} + \def \scalarprod #1#2{\left\langle#1,#2\right\rangle} + \def \timeStep {\delta_t} + \def \nrSnapshots {M} + \def \nrSnapshotsGen {\tilde{\nrSnapshots}} + \def \state {x} + \def \stateVar {\state} + \def \stateDim {n} + \def \stateDimRed {r} + \def \inpVar {u} + \def \inputVar {\inpVar} + \def \inpVarDim {m} + \def \inputDim {\inpVarDim} + \def \outVar {y} + \def \outputVar {\outVar} + \def \dmdV {X_+} + \def \dmdW {X_{-}} + \def \dmdY {Y} + \def \dmdU {U} + \def \dataZ {\mathcal{Z}} + \def \dataT {\mathcal{T}} + \def \dmdJJ {\widetilde{\mathcal{J}}} + \def \dmdRR {\widetilde{\mathcal{R}}} +""" + +# Translate LaTeX macros to KaTeX and add to options for HTML builder +katex_macros = katex.latex_defs_to_katex_macros(latex_macros) +katex_options = 'macros: {' + katex_macros + '}' + +# Add LaTeX macros for LATEX builder +latex_elements = {'preamble': latex_macros} + +bibtex_bibfiles = ['refs.bib'] + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['_templates'] + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This pattern also affects html_static_path and html_extra_path. +exclude_patterns = [] + + +# -- Options for HTML output ------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# +html_theme = 'sphinx_rtd_theme' + +# # Material theme options (see theme.conf for more information) +# html_theme_options = { +# +# # Set the name of the project to appear in the navigation. +# 'nav_title': 'pHDMD', +# +# # Set you GA account ID to enable tracking +# 'google_analytics_account': 'UA-XXXXX', +# +# # Specify a base_url used to generate sitemap.xml. If not +# # specified, then no sitemap will be built. +# 'base_url': 'https://project.github.io/project', +# +# # Set the color and the accent color +# 'color_primary': 'blue', +# 'color_accent': 'light-blue', +# +# # Set the repo location to get a badge with stats +# # TODO: change when release +# 'repo_url': 'https://github.com/Jonas-Nicodemus/crispy-pancake', +# 'repo_name': 'crispy-pancake', +# +# # Visible levels of the global TOC; -1 means unlimited +# 'globaltoc_depth': 3, +# # If False, expand all TOC entries +# 'globaltoc_collapse': False, +# # If True, show hidden TOC entries +# 'globaltoc_includehidden': False, +# } + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ['_static'] \ No newline at end of file diff --git a/docs/source/contact.rst b/docs/source/contact.rst new file mode 100644 index 0000000..7a460e6 --- /dev/null +++ b/docs/source/contact.rst @@ -0,0 +1,10 @@ +Contact +======= + +Jonas Nicodemus - jonas.nicodemus@simtech.uni-stuttgart.de + +Benjamin Unger - benjamin.unger@simtech.uni-stuttgart.de + +Riccardo Morandin - morandin@math.tu-berlin.de + +Project link: https://github.com/Jonas-Nicodemus/phdmd \ No newline at end of file diff --git a/docs/source/index.rst b/docs/source/index.rst new file mode 100644 index 0000000..43258af --- /dev/null +++ b/docs/source/index.rst @@ -0,0 +1,46 @@ +.. pHDMD documentation master file, created by + sphinx-quickstart on Thu Apr 21 11:21:24 2022. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. + +Port-Hamiltonian Dynamic Mode Decomposition +=========================================== + +We present a novel physics-informed system identification method to construct +a passive linear time-invariant system. In more detail, for a given quadratic +energy functional, measurements of the input, state, and output of a system +in the time domain, we find a realization that approximates the data well +while guaranteeing that the energy functional satisfies a dissipation inequality. +To this end, we use the framework of port-Hamiltonian (pH) systems and modify +the dynamic mode decomposition to be feasible for continuous-time pH systems. +We propose an iterative numerical method to solve the corresponding least-squares +minimization problem. We construct an effective initialization of the algorithm +by studying the least-squares problem in a weighted norm, for which we present +the analytical minimum-norm solution. The efficiency of the proposed method is +demonstrated with several numerical examples. + +Citing +========= +If you use this project for academic work, please consider citing our +`Publication `_: + + R. Morandin, J. Nicodemus, and B. Unger + Port-Hamiltonian Dynamic Mode Decomposition + ArXiv e-print TODO, 2022. + + + +.. toctree:: + :maxdepth: 2 + :caption: Contents: + + modules + bibliography + contact + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` diff --git a/docs/source/linalg.rst b/docs/source/linalg.rst new file mode 100644 index 0000000..b531521 --- /dev/null +++ b/docs/source/linalg.rst @@ -0,0 +1,37 @@ +linalg package +============== + +Submodules +---------- + +linalg.definiteness module +-------------------------- + +.. automodule:: linalg.definiteness + :members: + :undoc-members: + :show-inheritance: + +linalg.svd module +----------------- + +.. automodule:: linalg.svd + :members: + :undoc-members: + :show-inheritance: + +linalg.symmetric module +----------------------- + +.. automodule:: linalg.symmetric + :members: + :undoc-members: + :show-inheritance: + +Module contents +--------------- + +.. automodule:: linalg + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/model.rst b/docs/source/model.rst new file mode 100644 index 0000000..9d6a5c5 --- /dev/null +++ b/docs/source/model.rst @@ -0,0 +1,29 @@ +model package +============= + +Submodules +---------- + +model.msd module +---------------- + +.. automodule:: model.msd + :members: + :undoc-members: + :show-inheritance: + +model.poro module +----------------- + +.. automodule:: model.poro + :members: + :undoc-members: + :show-inheritance: + +Module contents +--------------- + +.. automodule:: model + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/modules.rst b/docs/source/modules.rst new file mode 100644 index 0000000..58d11c3 --- /dev/null +++ b/docs/source/modules.rst @@ -0,0 +1,11 @@ +API Reference +============= + +.. toctree:: + :maxdepth: 4 + + algorithm + linalg + model + system + utils diff --git a/docs/source/refs.bib b/docs/source/refs.bib new file mode 100755 index 0000000..5b49a9d --- /dev/null +++ b/docs/source/refs.bib @@ -0,0 +1,23 @@ +@STRING{MathCompModDynSys = {Math. Comput. Model. Dyn. Sys.}} %Mathematical and Computer Modelling of Dynamical Systems} + +@article{AltMU21, + author = {{Altmann}, R. and {Mehrmann}, V. and {Unger}, B.}, + title = {Port-{H}amiltonian formulations of poroelastic network models}, + journal = MathCompModDynSys, + year = {2021}, + doi = {10.1080/13873954.2021.1975137}, + volume = {27}, + number = {1}, + pages = {429--452}, +} + +@Article{GugPBV12, + author = {Gugercin, S. and Polyuga, Rostyslav V. and Beattie, C. and {van der Schaft}, A.}, + title = {Structure-preserving tangential interpolation for model reduction of port-{H}amiltonian systems}, + journal = {Automatica}, + volume = {48}, + number = {9}, + pages = {1963--1974}, + year = {2012}, + doi = {10.1016/j.automatica.2012.05.052}, +} \ No newline at end of file diff --git a/docs/source/system.rst b/docs/source/system.rst new file mode 100644 index 0000000..2432f41 --- /dev/null +++ b/docs/source/system.rst @@ -0,0 +1,21 @@ +system package +============== + +Submodules +---------- + +system.phlti module +------------------- + +.. automodule:: system.phlti + :members: + :undoc-members: + :show-inheritance: + +Module contents +--------------- + +.. automodule:: system + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/utils.rst b/docs/source/utils.rst new file mode 100644 index 0000000..1033bb7 --- /dev/null +++ b/docs/source/utils.rst @@ -0,0 +1,29 @@ +utils package +============= + +Submodules +---------- + +utils.norm module +----------------- + +.. automodule:: utils.norm + :members: + :undoc-members: + :show-inheritance: + +utils.plotting module +--------------------- + +.. automodule:: utils.plotting + :members: + :undoc-members: + :show-inheritance: + +Module contents +--------------- + +.. automodule:: utils + :members: + :undoc-members: + :show-inheritance: diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..b18fa54 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,5 @@ +numpy >= 1.22.3 +scipy >= 1.8.0 +matplotlib >= 3.5.1 +pymor >= 2021.2.1 +slycot >= 0.4.0 \ No newline at end of file diff --git a/src/algorithm/__init__.py b/src/algorithm/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/algorithm/phdmd.py b/src/algorithm/phdmd.py new file mode 100644 index 0000000..b3933e5 --- /dev/null +++ b/src/algorithm/phdmd.py @@ -0,0 +1,224 @@ +import logging + +import numpy as np + +from algorithm.skew_procrustes import skew_procrustes +from linalg.definiteness import project_spsd +from linalg.svd import svd +from linalg.symmetric import skew + + +def phdmd(X, Y, U, delta_t, E, max_iter=20, delta=1e-10): + r""" + The pHDMD algorithm identifies a port-Hamiltonian system from state, output and input measurements. + + Define + + .. math:: + \begin{align*} + \dmdW &\vcentcolon= \tfrac{1}{\timeStep} + \begin{bmatrix} + \state_1-\state_0 & \ldots & \state_{\nrSnapshots} - \state_{\nrSnapshots-1} + \end{bmatrix}\in\R^{\stateDim\times\nrSnapshots},\\ + \dmdV &\vcentcolon= \tfrac{1}{2}\begin{bmatrix} + \state_1+\state_0 & \ldots & \state_{\nrSnapshots} + \state_{\nrSnapshots-1} + \end{bmatrix}\in\R^{\stateDim\times\nrSnapshots},\\ + \dmdU &\vcentcolon= \tfrac{1}{2}\begin{bmatrix} + \inpVar_1 + \inpVar_0 & \ldots & \inpVar_{\nrSnapshots} + \inpVar_{\nrSnapshots-1} + \end{bmatrix} \in\R^{\inpVarDim\times\nrSnapshots},\\ + \dmdY &\vcentcolon= \tfrac{1}{2}\begin{bmatrix} + \outVar_1 + \outVar_0 & \ldots & \outVar_{\nrSnapshots} + \outVar_{\nrSnapshots-1} + \end{bmatrix}\in\R^{\inpVarDim\times\nrSnapshots}, + \end{align*} + + and + + .. math:: + \dataZ \vcentcolon= \begin{bmatrix}H \dmdW \\ -\dmdY\end{bmatrix} \quad + \mathrm{and} \quad \dataT \vcentcolon= \begin{bmatrix}\dmdV \\ \dmdU\end{bmatrix}. + + Solve the minimization problem + + .. math:: + \min_{\mathcal{J},\mathcal{R}} \|Z - (\mathcal{J} - \mathcal{R}) T\|_F, + \quad \mathrm{s.t.} \quad \mathcal{J}=-\mathcal{J}^T, \mathcal{R}\in\mathbb{S}^{n}_{\succeq}. + + Parameters + ---------- + X : numpy.ndarray + Sequence of states. + Y : numpy.ndarray + Sequence of outputs. + U : numpy.ndarray + Sequence of inputs. + delta_t : float + Time step size. + E : numpy.ndarray + Hamiltonian matrix. + max_iter : int, optional + Maximum number of iterations. Default 20. + delta : float, optional + Convergence criteria. Default 1e-10. + + Returns + ------- + J : numpy.ndarray + Conservation of energy matrix. + R : numpy.ndarray + Dissipation matrix. + """ + W = (X[:, 1:] - X[:, :-1]) / delta_t + V = 1 / 2 * (X[:, 1:] + X[:, :-1]) + U = 1 / 2 * (U[:, 1:] + U[:, :-1]) + Y = 1 / 2 * (Y[:, 1:] + Y[:, :-1]) + + T = np.concatenate((V, U)) + Z = np.concatenate((E @ W, -Y)) + + # Fantastic initialization + J, R = init_phdmd(T, Z) + + J, R, e = phdmd_FGM(T, Z, J, R, max_iter, delta) + + return J, R + + +def init_phdmd(T, Z, tol=1e-12): + r""" + Returns an initialization for the pHDMD algorithm by solving the related weighted minimization problem + + .. math:: + \min_{\mathcal{J},\mathcal{R}} \|T^TZ - T^T(\mathcal{J} - \mathcal{R}) T\|_F, + \quad \mathrm{s.t.} \quad \mathcal{J}=-\mathcal{J}^T, \mathcal{R}\in\mathbb{S}^{n}_{\succeq}. + + Parameters + ---------- + T : numpy.ndarray + Stacked data matrix + Z : numpy.ndarray + Stacked data matrix + tol : float, optional + Zero tolerance. + + Returns + ------- + J : numpy.ndarray + Conservation of energy matrix. + R : numpy.ndarray + Dissipation matrix. + """ + + n, m = T.shape + + U, s, V = svd(T) + r = np.argmax(s / s[0] < tol) + r = r if r > 0 else len(s) + s = s[:r] + + S = np.diag(s) + S_inv = np.diag(1 / s) + + if r < n: + logging.warning(f'Rank(T) < n + m ({r} < {n})') + + Z_1 = U.T @ Z @ V + + J_11 = skew(S @ Z_1[:r, :r]) + R_11 = project_spsd(-S @ Z_1[:r, :r]) + + R = U[:, :r] @ S_inv @ R_11 @ S_inv @ U[:, :r].T + J = U[:, :r] @ S_inv @ J_11 @ S_inv @ U[:, :r].T + + if r < n: + # compensate rank deficit + J_21 = np.linalg.lstsq(np.diag(s), Z_1[r:, :r].T, rcond=None)[0].T + J_cmp = np.zeros((n, n)) + J_cmp[r:, :r] = J_21 + J_cmp[:r, r:] = -J_21.T + J = J + U @ J_cmp @ U.T + + return J, R + + +def phdmd_FGM(T, Z, J0, R0, max_iter=20, delta=1e-10): + r""" + Iterative algorithm to solve the pHDMD problem via a fast-gradient method. + + Parameters + ---------- + T : numpy.ndarray + Stacked data matrix. + Z : numpy.ndarray + Stacked data matrix. + J0 : numpy.ndarray + Initial matrix for `J`. + R0 : numpy.ndarray + Initial matrix for `R`. + max_iter : int, optional + Maximum number of iterations. Default 20. + delta : float, optional + Convergence criteria. Default 1e-10. + + Returns + ------- + J : numpy.ndarray + Conservation of energy matrix. + R : numpy.ndarray + Dissipation matrix. + e : numpy.ndarray + Value of the cost functional at each iteration. + """ + R = R0 + J = J0 + + # Precomputations + TTt = T @ T.T + w, _ = np.linalg.eigh(TTt) + L = max(w) # Lipschitz constant + mu = min(w) + q = mu / L + + beta = np.zeros(max_iter) + alpha = np.zeros(max_iter) + e = np.zeros(max_iter) + + # Parameters and initialization + alpha_0 = 0.1 # Parameter of the FGM in (0,1) - can be tuned. + + Q = R + alpha[0] = alpha_0 + e[0] = np.linalg.norm(Z - (J - R) @ T, 'fro') / np.linalg.norm(Z) + logging.info(f'|Z - (J^(0) - R^(0) T|_F / |Z|_F={e[0]:.2e}') + + for i in range(max_iter): + # Previous iterate + Rp = R + Jp = J + + Z_1 = Z + R @ T + # Solution of the skew-symmetric Procrustes + J, _ = skew_procrustes(T, Z_1) + + Z_2 = J @ T - Z + # Projected gradient step from Y + G = Q @ TTt - Z_2 @ T.T + R = project_spsd(Q - G / L) + + # FGM Coefficients + alpha[i + 1] = (np.sqrt((alpha[i] ** 2 - q) ** 2 + 4 * alpha[i] ** 2) + (q - alpha[i] ** 2)) / 2 + beta[i + 1] = alpha[i] * (1 - alpha[i]) / (alpha[i] ** 2 + alpha[i + 1]) + + # Linear combination of iterates + Q = R + beta[i] * (R - Rp) + + e[i + 1] = np.linalg.norm(Z - (J - R) @ T, 'fro') / np.linalg.norm(Z) + logging.info(f'|Z - (J^({i + 1}) - R^({i + 1}) T|_F / |Z|_F={e[i + 1]:.2e}') + + eps = np.linalg.norm(Jp - J, 'fro') / (np.linalg.norm(J, 'fro')) + \ + np.linalg.norm(Rp - R, 'fro') / (np.linalg.norm(R, 'fro')) + if eps < delta or e[i] - e[i + 1] < delta: + e = e[:i+2] + logging.info(f'Converged after {i + 1} iterations.') + break + + return J, R, e diff --git a/src/algorithm/skew_procrustes.py b/src/algorithm/skew_procrustes.py new file mode 100644 index 0000000..5c22260 --- /dev/null +++ b/src/algorithm/skew_procrustes.py @@ -0,0 +1,81 @@ +import logging + +import numpy as np + +from linalg.svd import svd + + +def skew_procrustes(X, Y, trunc_tol=1e-12, zero_tol=1e-14): + r""" + Solves the skew-symmetric Procrustes problem + + .. math:: + \min_{A} \|Y - AX\|_F \quad s.t. \quad A = -A^T. + + Parameters + ---------- + X : numpy.ndarray + Matrix. + Y : numpy.ndarray + Matrix. + Returns + ------- + A : numpy.ndarray + The minimizer of the skew-symmetric Procrustes problem. + e : float + Error of the skew-symmetric Procrustes problem. + """ + + n = X.shape[0] + U, s, V = svd(X) + r = np.argmax(s / s[0] < trunc_tol) + r = r if r > 0 else len(s) + + s = s[:r] + + Z_transform = U.T @ Y @ V + Z1 = Z_transform[:r, :r] + Z3 = Z_transform[r:, :r] + + # Assemble Phi + [ii, jj] = np.mgrid[:r, :r] + tmp = np.square(s[ii]) + np.square(s[jj]) + Phi = 1. / tmp + + A1 = Phi * (Z1 @ np.diag(s) - np.diag(s) @ Z1.T) + + A2 = -np.linalg.solve(np.diag(s), Z3.T) + A4 = np.zeros((n - r, n - r)) + + A = U @ np.array(np.concatenate( + (np.hstack((A1, A2)), np.hstack((-A2.T, A4))) + )) @ U.T + + max_A = np.amax(np.abs(A)) + indices = np.where(np.abs(A) / max_A < zero_tol) + A[indices] = 0 + + e = np.linalg.norm(A @ X - Y, 'fro') + + return A, e + + +if __name__ == "__main__": + logging.basicConfig() + logging.getLogger().setLevel(logging.INFO) + + n = 6 + M = 500 + A = np.random.random((n, n)) + A = (A - A.T) / 2 + + X = np.random.random((n, M)) + r = np.linalg.matrix_rank(X) + + Y = A @ X + + X_res, res = skew_procrustes(X, Y) + + logging.info(f'Error of the skew symmetric Procrustes ||A @ X - Y||_F:\t{res:.2e}') + + assert res <= 1e-2 diff --git a/src/algorithm/spsd_procrustes.py b/src/algorithm/spsd_procrustes.py new file mode 100644 index 0000000..67508f9 --- /dev/null +++ b/src/algorithm/spsd_procrustes.py @@ -0,0 +1,129 @@ +import logging + +import matplotlib.pyplot as plt +import numpy as np + +from linalg.definiteness import project_spsd + + +def spsd_procrustes(X, Y, A0=None, max_iter=1000): + """ + Solves the symmetric positive semi-definite Procrustes problem via the fast-gradient method + + .. math:: + \min_{A} \| Y - A X \|_F \quad s.t. \quad A\in\mathbb{S}^{n}_{\succeq}. + + Parameters + ---------- + X : numpy.ndarray + Matrix. + Y : numpy.ndarray + Matrix. + A0 : numpy.ndarraym, optional + Initial matrix for `A`. Default `None`. + max_iter : int, optional + Maximum number of iterations. Default 1000. + + Returns + ------- + A : numpy.ndarray + Symmetric positive-definite matrix that (most) closely maps `X` to `Y`. + e : float + Remaining error of the symmetric positive-definite Procrustes problem. + """ + if A0 is None: + A0, _ = init_spsd_procrustes(X, Y) + + A = A0 + + # Precomputations + XXt = X @ X.T + x, _ = np.linalg.eigh(XXt) + Lx = max(x) # Lipschitz constant + mux = min(x) + qx = mux / Lx + YXt = Y @ X.T + + # Parameters and initalization + alpha0 = 0.1 # Parameter of the FGM in (0,1) -it can be tuned. + + beta = np.zeros(max_iter) + alpha = np.zeros(max_iter + 1) + e = np.zeros(max_iter + 1) + + B = A + alpha[0] = alpha0 + e[0] = np.linalg.norm(A @ X - Y, 'fro') + + for i in range(max_iter): + # Previous iterate + Ap = A + + # Projected gradient step from Y + A = project_spsd(B - (B @ XXt - YXt) / Lx) + + # FGM Coefficients + alpha[i + 1] = (np.sqrt((alpha[i] ** 2 - qx) ** 2 + 4 * alpha[i] ** 2) + (qx - alpha[i] ** 2)) / 2 + beta[i] = alpha[i] * (1 - alpha[i]) / (alpha[i] ** 2 + alpha[i + 1]) + + # Linear combination of iterates + B = A + beta[i] * (A - Ap) + + e[i + 1] = np.linalg.norm(A @ X - Y, 'fro') + + return A, e + + +def init_spsd_procrustes(X, Y): + """ + Returns an initialization of the symmetric positive semi-definite Procrustes problem. + + Parameters + ---------- + X : numpy.ndarray + Matrix. + Y : numpy.ndarray + Matrix. + + Returns + ------- + A0 : numpy.ndarray + Initial matrix for `A` for the Symmetric positive-definite Procrustes problem. + e0 : float + Error for the initial matrix `A` of the symmetric positive-definite Procrustes problem. + """ + n = Y.shape[0] + + A0 = np.zeros((n, n)) + for i in range(n): + A0[i, i] = max(0, X[i, :] @ Y[i, :].T) / np.linalg.norm(X[i, :]) ** 2 + 1e-6 + + e0 = np.linalg.norm(A0 @ X - Y, 'fro') + + return A0, e0 + + +if __name__ == "__main__": + logging.basicConfig() + logging.getLogger().setLevel(logging.INFO) + + n = 100 + m = 1200 + + A = np.random.randn(n, n) + A = 1 / 2 * (A + A.T) + A = A + n * np.eye(n) + + X = np.random.randn(n, m) + Y = A @ X + + A0 = None + + A_proc, e = spsd_procrustes(X, Y, A0=A0, max_iter=200) + + e_rel = e / np.linalg.norm(Y, 'fro') + plt.semilogy(e_rel) + plt.show() + + logging.info(f'||Y - AX||_F = {e[-1]:.2e}') + logging.info(f'||Y - AX||_F / ||Y||_F = {e_rel[-1]:.2e}') diff --git a/src/config.py b/src/config.py new file mode 100644 index 0000000..0d294eb --- /dev/null +++ b/src/config.py @@ -0,0 +1,90 @@ +import os + +import numpy as np + +from scipy.signal import sawtooth + +from model.msd import msd +from model.poro import poro + +save_results = False # If true all figures will be saved as pgf plots +figures_path = os.path.join('../figures') + +match 1: # 1 - SISO MSD with few data, 2 - MIMO MSD with few data, 3 - MIMO MSD with many data, 4 - PORO + case 1: # SISO MSD with few data + model = 'msd' # msd, poro + n = 6 + m = 1 + ph_matrices = lambda: msd(n, m) + u = lambda t: np.exp(-0.5 * t) * np.sin(t ** 2) + exp_id = model + f'_n_{n}_m_{m}' + + x0 = np.zeros(n) + N = 100 + T_end = 4 + T = np.linspace(0, T_end, N) + delta = T[1] - T[0] + + T_test = np.linspace(0, 10, 1000) + freq = 0.5 # Hz + u_test = lambda t: sawtooth(2 * np.pi * freq * t) + + case 2: # MIMO MSD with few data + model = 'msd' # msd, poro + n = 100 + m = 2 + ph_matrices = lambda: msd(n, m) + factor = 100 + u = lambda t: np.array([np.exp(-0.5 * t) * np.sin(t ** 2), + np.exp(-0.5 * t) * np.cos(t ** 2)]) + exp_id = model + f'_n_{n}_m_{m}_short' + + x0 = np.zeros(n) + N = 100 + T_end = 4 + T = np.linspace(0, T_end, N) + delta = T[1] - T[0] + + T_test = np.linspace(0, 10, 1000) + freq = 0.5 # Hz + u_test = lambda t: np.array([sawtooth(2 * np.pi * freq * t), - sawtooth(2 * np.pi * freq * t)]) + + case 3: # MIMO MSD with many data + model = 'msd' # msd, poro + n = 100 + m = 2 + ph_matrices = lambda: msd(n, m) + factor = 100 + u = lambda t: np.array([np.exp(-0.5 / factor * t) * np.sin(1 / factor * t ** 2), + np.exp(-0.5 / factor * t) * np.cos(1 / factor * t ** 2)]) + exp_id = model + f'_n_{n}_m_{m}_long' + + x0 = np.zeros(n) + N = 100 * factor + T_end = 4 * factor + T = np.linspace(0, T_end, N) + delta = T[1] - T[0] + + T_test = np.linspace(0, 10, 1000) + freq = 0.5 # Hz + u_test = lambda t: np.array([sawtooth(2 * np.pi * freq * t), - sawtooth(2 * np.pi * freq * t)]) + + case 4: # PORO + model = 'poro' # msd, poro + n = 980 + m = 2 + ph_matrices = lambda: poro(n) + factor = 100 + u = lambda t: np.array([np.exp(-0.5 / factor * t) * np.sin(1 / factor * t ** 2), + np.exp(-0.5 / factor * t) * np.cos(1 / factor * t ** 2)]) + exp_id = model + f'_n_{n}_m_{m}' + + x0 = np.zeros(n) + N = 100 * factor + T_end = 4 * factor + T = np.linspace(0, T_end, N) + delta = T[1] - T[0] + + T_test = np.linspace(0, 10, 10000) + freq = 0.5 # Hz + u_test = lambda t: np.array([sawtooth(2 * np.pi * freq * t), - sawtooth(2 * np.pi * freq * t)]) diff --git a/src/linalg/__init__.py b/src/linalg/__init__.py new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/src/linalg/__init__.py @@ -0,0 +1 @@ + diff --git a/src/linalg/definiteness.py b/src/linalg/definiteness.py new file mode 100644 index 0000000..30ca2e7 --- /dev/null +++ b/src/linalg/definiteness.py @@ -0,0 +1,51 @@ +import numpy as np + +from linalg.symmetric import is_sym, sym + + +def is_spsd(A): + """ + Checks if a matrix is symmetric positive semi-definite. + + Parameters + ---------- + A : numpy.ndarray + The matrix to check for symmetric positive semi-definiteness. + + Returns + ------- + is_spsd : bool + `True` if `A` is symmetric positive semi-definiteness, else `False`. + """ + if is_sym(A): + try: + np.linalg.cholesky(A) + return True + except np.linalg.LinAlgError: + A += 1e-5 * np.eye(A.shape[0]) + try: + np.linalg.cholesky(A) + return True + except np.linalg.LinAlgError: + return False + else: + return False + + +def project_spsd(A): + """ + Projects a matrix on the set of symmetric positive semi-definite matrices. + + Parameters + ---------- + A : array_like + The matrix to project on the set of symmetric positive semi-definite matrices. + + Returns + ------- + A_spsd : numpy.ndarray + The projected symmetric positive semi-definite matrix. + """ + w, U = np.linalg.eigh(sym(A)) + + return U @ np.diag(w.clip(min=0)) @ U.T diff --git a/src/linalg/svd.py b/src/linalg/svd.py new file mode 100644 index 0000000..7c93812 --- /dev/null +++ b/src/linalg/svd.py @@ -0,0 +1,102 @@ +import numpy as np + + +def svd(A, hermitian=False, variant=None, tol=1e-12, t=None): + """ + Wrapper around `numpy.svd` to perform the four variants of the singular value decomposition + (https://en.wikipedia.org/wiki/Singular_value_decomposition). + + Parameters + ---------- + A : numpy.ndarray + A real or complex matrix to perform the svd on. + hermitian : bool, optional + If True, `A` is assumed to be Hermitian (symmetric if real-valued), + enabling a more efficient method for finding singular values. + Defaults to False. + variant : str, optional + Variant of the svd. Should be one of + + - 'full' + - 'thin' + - 'skinny' + - 'trunc' + tol : float, optional + In the case of 'skinny' the zero tolerance and in the case of 'trunc' the truncation tolerance. + t : int, optional + In the case of 'trunc' the truncation index. + + Returns + ------- + U : numpy.ndarray + Left singular vectors. + s : numpy.ndarray + Singular values. + V : numpy.ndarray + Right singular vectors. + """ + U, s, Vh = np.linalg.svd(A, hermitian=hermitian) + V = Vh.conj().T + + match variant: + case 'full': + pass + case 'thin': + n, m = A.shape + k = min(n, m) + U = U[:, :k] + V = V[:, :k] + case 'skinny': + r = np.argmax(s / s[0] < tol) + r = r if r > 0 else len(s) + U = U[:, :r] + V = V[:, :r] + s = s[:r] + case 'trunc': + t = np.argmax(s / s[0] < tol) if t is None else t + t = t if t > 0 else len(s) + U = U[:, :t] + V = V[:, :t] + s = s[:t] + case None: + pass + + return U, s, V + + +if __name__ == "__main__": + n = 10 + m = 100 + A = np.random.randint(10, size=(n, m)) + A[-1] = 0 + r = np.linalg.matrix_rank(A) + assert r == 9 + + U, s, V = svd(A) + assert U.shape == (n, n) + assert V.shape == (m, m) + assert len(s) == n + + U, s, V = svd(A, variant='thin') + assert U.shape == (n, n) + assert V.shape == (m, n) + assert len(s) == n + + U, s, V = svd(A, variant='skinny') + assert U.shape == (n, r) + assert V.shape == (m, r) + assert len(s) == r + + t = 5 + U, s, V = svd(A, variant='trunc', t=t) + assert U.shape == (n, t) + assert V.shape == (m, t) + assert len(s) == t + + + + + + + + diff --git a/src/linalg/symmetric.py b/src/linalg/symmetric.py new file mode 100644 index 0000000..78da3ad --- /dev/null +++ b/src/linalg/symmetric.py @@ -0,0 +1,105 @@ +import numpy as np + + +def toeplitz(A): + r""" + Toeplitz decomposition. + + Returns the symmetric and skew-symmetric part, + + .. math:: + + S = \tfrac{1}{2} (A + A^T),\quad N = \tfrac{1}{2} (A - A^T) + + of a square matrix :math:`A`. + + Parameters + ---------- + A : numpy.ndarray + Square matrix to split in symmetric and skew-symmetric part. + + Returns + ------- + S : numpy.ndarray + Symmetric part of `A`. + N : numpy.ndarray + Skew-symmetric part of `A`. + """ + return sym(A), skew(A) + + +def sym(A): + """ + Symmetric part of a square matrix `A`. + + Parameters + ---------- + A : numpy.ndarray + Square matrix. + + Returns + ------- + S : numpy.ndarray + Symmetric part of `A`. + """ + return 1 / 2 * (A + A.T) + + +def skew(A): + """ + Skew-symmetric part of a square matrix `A`. + + Parameters + ---------- + A : numpy.ndarray + Square matrix. + + Returns + ------- + N : numpy.ndarray + Skew-symmetric part of `A`. + """ + return 1 / 2 * (A - A.T) + + +def is_sym(A, rtol=1e-05, atol=1e-08): + """ + Checks if a square matrix `A` is symmetric. + + Parameters + ---------- + A : numpy.ndarray + Square matrix. + + Returns + ------- + is_sym : bool + `True` if `A` is symmetric, else `False`. + """ + return np.allclose(A, A.T, rtol=rtol, atol=atol) + + +def is_skew(A, rtol=1e-05, atol=1e-08): + """ + Checks if a square matrix `A` is skew-symmetric. + + Parameters + ---------- + A : numpy.ndarray + Square matrix. + + Returns + ------- + is_sym : bool + `True` if `A` is skew-symmetric, else `False`. + """ + return np.allclose(A, -A.T, rtol=rtol, atol=atol) + + +if __name__ == "__main__": + n = 5 + A = np.random.random((n, n)) + + A_sym, A_skew = toeplitz(A) + + assert is_sym(A_sym) and is_skew(A_skew) diff --git a/src/main.py b/src/main.py new file mode 100644 index 0000000..70542ce --- /dev/null +++ b/src/main.py @@ -0,0 +1,106 @@ +import logging +import os +import config + +import numpy as np + +from algorithm.phdmd import phdmd +from system.phlti import PHLTI +from utils.norm import lnorm +from utils.plotting import trajectories_plot, abs_error_plot, bode_plot, magnitude_plot, poles_plot + + +def main(): + logging.basicConfig() + logging.getLogger().setLevel(logging.INFO) + + if config.save_results: + import matplotlib + matplotlib.use("pgf") + matplotlib.rcParams.update({ + "pgf.texsystem": "pdflatex", + 'font.family': 'serif', + 'text.usetex': True, + 'pgf.rcfonts': False, + }) + + if not os.path.exists(config.figures_path): + os.makedirs(config.figures_path) + + # Initialize the original ph system + E, J, R, G, P, D, Q = config.ph_matrices() + + ph = PHLTI(E, J, R, G, P, D, Q) + lti = ph.to_lti(matrices=False) + + # Generate training data + U_train, X_train, Y_train = ph.sim(config.u, config.T, config.x0) + + trajectories_plot(config.T, U_train, label='u', train=True) + trajectories_plot(config.T, Y_train, train=True) + + # Perform pHDMD + J_dmd, R_dmd = phdmd(X_train, Y_train, U_train, config.delta, + E=ph.E, max_iter=100) + ph_dmd = PHLTI(ph.E, J_dmd, R_dmd) + + # Testing + U_train_dmd, X_train_dmd, Y_train_dmd = ph_dmd.sim(config.u, config.T, config.x0) + U_test, X_test, Y_test = ph.sim(config.u_test, config.T_test, config.x0) + U_dmd, X_dmd, Y_dmd = ph_dmd.sim(config.u_test, config.T_test, config.x0) + + # Error lti system + lti_dmd = ph_dmd.to_lti(matrices=False) + lti_error = lti - lti_dmd + + # Trajectories + trajectories_plot(config.T_test, U_test, label='u') + trajectories_plot(config.T_test, Y_test, Y_dmd) + + # Absolute Error of the trajectories + abs_error_plot(config.T_test, Y_test, Y_dmd) + + # Only first part of the trajectories + trajectories_plot(config.T_test, Y_test, Y_dmd, zoom=100) + abs_error_plot(config.T_test, Y_test, Y_dmd, zoom=100) + + # Bode + w = np.logspace(-1, 3, 100) + bode_plot(w, lti, lti_dmd) + + # Magnitude error system + magnitude_plot(w, lti_error) + + # Poles + poles_plot(lti, lti_dmd) + + # Norms + l2_train = lnorm(Y_train - Y_train_dmd, p=2, tau=config.delta) + l2_rel_train = l2_train / lnorm(Y_train, p=2, tau=config.delta) + logging.info(f'Relative L2 error (Training data): {l2_rel_train:.2e}') + + linf_train = lnorm(Y_train - Y_train_dmd, p=np.inf, tau=config.delta) + linf_rel_train = linf_train / lnorm(Y_train, p=np.inf, tau=config.delta) + logging.info(f'Relative Linf error (Training data): {linf_rel_train:.2e}') + + l2 = lnorm(Y_test - Y_dmd, p=2, tau=config.delta) + l2_rel = l2 / lnorm(Y_test, p=2, tau=config.delta) + logging.info(f'Relative L2 error: {l2_rel:.2e}') + + linf_train = lnorm(Y_test - Y_dmd, p=np.inf, tau=config.delta) + linf_rel = linf_train / lnorm(Y_test, p=np.inf, tau=config.delta) + logging.info(f'Relative Linf error: {linf_rel:.2e}') + + # H-norm calculation fails for the poro benchmark system + if config.model != 'poro': + h2 = lti_error.h2_norm() + h2_rel = h2 / lti.h2_norm() + logging.info(f'Relative H2 error: {h2_rel:.2e}') + + hinf = lti_error.hinf_norm() + hinf_rel = hinf / lti.hinf_norm() + logging.info(f'Relative Hinf error: {hinf_rel:.2e}') + + +if __name__ == "__main__": + main() diff --git a/src/model/__init__.py b/src/model/__init__.py new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/src/model/__init__.py @@ -0,0 +1 @@ + diff --git a/src/model/msd.py b/src/model/msd.py new file mode 100644 index 0000000..f654da2 --- /dev/null +++ b/src/model/msd.py @@ -0,0 +1,84 @@ +import numpy as np + +from scipy.linalg import block_diag + + +def msd(n=6, m=2, m_i=4, k_i=4, c_i=1): + """ + Returns the mass-spring-damper benchmark system (cf. :cite:`GugPBV12`), as port-Hamiltonian system. + + Parameters + ---------- + n : int, optional + Dimension of the state. + m : int, optional + Dimension of the input resp. output. (only 1 or 2 are implemented) + m_i : float, optional + Weight of the masses. + k_i : float, optional + Stiffness of the springs. + c_i : float, optional + Amount of damping. + + Returns + ------- + E : numpy.ndarray + Hamiltonian matrix. + J : numpy.ndarray + Skew-symmetric matrix. + R : numpy.ndarray + Dissipative matrix. + G : numpy.ndarray + Input matrix. + P : numpy.ndarray + Input matrix. + D : numpy.ndarray + Feed trough matrix. + Q : numpy.ndarray + Hamiltonian matrix. + """ + assert (n % 2) == 0 + + n = int(n / 2) + + if m == 2: + B = np.array([[0, 1, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0]]).T + C = np.array([[0, 1 / m_i, 0, 0, 0, 0], [0, 0, 0, 1 / m_i, 0, 0]]) + elif m == 1: + B = np.array([[0, 1, 0, 0, 0, 0]]).T + C = np.array([[0, 1 / m_i, 0, 0, 0, 0]]) + else: + assert False + + A = np.array( + [[0, 1 / m_i, 0, 0, 0, 0], [-k_i, -c_i / m_i, k_i, 0, 0, 0], + [0, 0, 0, 1 / m_i, 0, 0], [k_i, 0, -2 * k_i, -c_i / m_i, k_i, 0], + [0, 0, 0, 0, 0, 1 / m_i], [0, 0, k_i, 0, -2 * k_i, -c_i / m_i]]) + + J_i = np.array([[0, 1], [-1, 0]]) + J = np.kron(np.eye(3), J_i) + R_i = np.array([[0, 0], [0, c_i]]) + R = np.kron(np.eye(3), R_i) + + G = B + for i in range(4, n + 1): + G = np.vstack((G, np.zeros((2, m)))) + C = np.hstack((C, np.zeros((m, 2)))) + J = block_diag(J, J_i) + R = block_diag(R, R_i) + + A = block_diag(A, np.zeros((2, 2))) + A[2 * i - 2, 2 * i - 2] = 0 + A[2 * i - 1, 2 * i - 1] = -c_i / m_i + A[2 * i - 3, 2 * i - 2] = k_i + A[2 * i - 2, 2 * i - 1] = 1 / m_i + A[2 * i - 2, 2 * i - 3] = 0 + A[2 * i - 1, 2 * i - 2] = -2 * k_i + A[2 * i - 1, 2 * i - 4] = k_i + + Q = np.linalg.solve(J - R, A) + P = np.zeros(G.shape) + D = np.zeros((m, m)) + E = np.eye(2 * n) + + return E, J, R, G, P, D, Q \ No newline at end of file diff --git a/src/model/poro.py b/src/model/poro.py new file mode 100644 index 0000000..4f7c778 --- /dev/null +++ b/src/model/poro.py @@ -0,0 +1,95 @@ +import os +import pkg_resources + +import numpy as np + +from scipy.io import loadmat + + +def poro(n=320): + """ + Returns a port-Hamiltonian model of linear poroelasticity in a + bounded Lipschitz domain as described in :cite:`AltMU21`. + + Parameters + ---------- + n : int, optional + System dimension (can only be either: 320, 980, or 1805). Default = 980. + + Returns + ------- + Returns + ------- + E : numpy.ndarray + Hamiltonian matrix. + J : numpy.ndarray + Skew-symmetric matrix. + R : numpy.ndarray + Dissipative matrix. + G : numpy.ndarray + Input matrix. + P : numpy.ndarray + Input matrix. + D : numpy.ndarray + Feed trough matrix. + Q : numpy.ndarray + Hamiltonian matrix. + """ + # load matrices + path = os.path.join(f'resources/poro-n{n}.mat') + path = pkg_resources.resource_stream(__name__, path) + + data = loadmat(path) + A = data['A'] + + # parameter settings + rho = 1e-3 + alpha = 0.79 + Minv = 7.80e3 + kappaNu = 633.33 + Rshift = 1e-3 + + # update discretization matrices to reflect the parameter settings + Y = rho * data['Y'] + D = alpha * data['D'] + M = Minv * data['M'] + K = kappaNu * data['K'] + Bp = data['Bp'].T + Bf = data['Bf'].T + + # construct first-order port-Hamiltonian descriptor system + # representation, similarly as in + # R.Altmann, V.Mehrmann, B.Unger: Port - Hamiltonian + # formulations of poroelastic network models, arXiv preprint + # arXiv: 2012.01949, 2020. + n = A.shape[0] + m = M.shape[0] + E = np.block([ + [Y, np.zeros((n, n + m))], + [np.zeros((n, n)), A, np.zeros((n, m))], + [np.zeros((m, n + n)), M] + ]) + + J = np.block([ + [np.zeros((n, n)), -A, D.T], + [A, np.zeros((n, n + m))], + [-D, np.zeros((m, n + m))] + ]) + + R = np.block([ + [np.zeros((n, 2 * n + m))], + [np.zeros((n, 2 * n + m))], + [np.zeros((m, 2 * n)), K] + ]) + Rshift * np.eye(2 * n + m) + + G = np.block([ + [Bf, np.zeros((n, 1))], + [np.zeros((n, 2))], + [np.zeros((m, 1)), Bp] + ]) + + P = np.zeros(G.shape) + D = np.zeros((G.shape[1], G.shape[1])) + Q = np.eye(J.shape[0]) + + return E, J, R, G, P, D, Q \ No newline at end of file diff --git a/src/model/resources/poro-n1805.mat b/src/model/resources/poro-n1805.mat new file mode 100644 index 0000000..5385275 Binary files /dev/null and b/src/model/resources/poro-n1805.mat differ diff --git a/src/model/resources/poro-n320.mat b/src/model/resources/poro-n320.mat new file mode 100644 index 0000000..9c4bb05 Binary files /dev/null and b/src/model/resources/poro-n320.mat differ diff --git a/src/model/resources/poro-n980.mat b/src/model/resources/poro-n980.mat new file mode 100644 index 0000000..ed8359b Binary files /dev/null and b/src/model/resources/poro-n980.mat differ diff --git a/src/system/__init__.py b/src/system/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/system/phlti.py b/src/system/phlti.py new file mode 100644 index 0000000..a43092c --- /dev/null +++ b/src/system/phlti.py @@ -0,0 +1,189 @@ +import numpy as np + +from pymor.models.iosys import LTIModel + +from linalg.symmetric import sym, skew + + +class PHLTI(object): + """ + Class representing a linear time invariant port-Hamiltonian system. + """ + def __init__(self, E, J, R, G=None, P=None, D=None, Q=None): + """ + Constructor. + + Parameters + ---------- + E : numpy.ndarray + Hamiltonian matrix. + J : numpy.ndarray + Skew-symmetric matrix. + R : numpy.ndarray + Dissipative matrix. + G : numpy.ndarray, optional + Input matrix. If `None` it is assumed that `J` and `R` contain already all other system matrices. + P : numpy.ndarray + Input matrix. If `None` it is assumed that `J` and `R` contain already all other system matrices. + D : numpy.ndarray + Feed trough matrix. If `None` it is assumed that `J` and `R` contain already all other system matrices. + Q : numpy.ndarray + Hamiltonian matrix. If `None` it is assumed that `Q` ist the identity i.e., `Q` already on the lhs, + otherwise `Q` will be brought to the lhs by multiplying with `Q.T` from left. + """ + if G is None: + n = E.shape[0] + P = R[:n, n:] + S = R[n:, n:] + + G = J[:n, n:] + N = J[n:, n:] + + R = R[:n, :n] + J = J[:n, :n] + D = S + N + + if Q is None: + # Q already on lhs + self.E = E + self.J = J + self.R = R + self.G = G + self.P = P + self.D = D + else: + # bring Q on lhs + self.E = Q.T @ E + self.J = Q.T @ J @ Q + self.R = Q.T @ R @ Q + self.G = Q.T @ G + self.P = Q.T @ P + self.D = D + + self.S = sym(self.D) + self.N = -skew(self.D) + + self.n, self.m = G.shape + self.p = self.m + + def sim(self, U, T=None, x0=None): + """ + Simulate the system for a given input signal and returns the resulting discrete input, state and output. + + Parameters + ---------- + U : numpy.ndarray, callable + Input signal as `numpy.ndarray` or `callable` object. + T : numpy.ndarray, optional + Time steps at which the input is defined and the output und state is returned. + If `None` the input signal `U` is assumed to be an `numpy.ndarray` and the time steps are calculated from it + with a step-size of 1. + x0 : numpy.ndarray, optional + Initial condition of the state. If `None` the initial condition is set to zero for all states. + + Returns + ------- + U : numpy.ndarray + Input values at the time steps `T`. + X : numpy.ndarray + States at the time steps `T`. + Y : numpy.ndarray + Output values at the time steps `T`. + """ + if T is None: + assert isinstance(U, np.ndarray) + T = np.linspace(0, len(U)) + if not isinstance(U, np.ndarray): + U = U(T) + if U.ndim < 2: + U = U[np.newaxis, :] + + if x0 is None: + x0 = np.zeros(self.n) + + return self.implicit_midpoint(U, T, x0) + + def implicit_midpoint(self, U, T, x0): + """ + Simulate the system via the implicit midpoint rule. + + Parameters + ---------- + U : numpy.ndarray + Input signal. + T : numpy.ndarray + Time steps at which the input is defined and the output und state is returned. + x0 : numpy.ndarray + Initial condition of the state. + + Returns + ------- + U : numpy.ndarray + Input values at the time steps `T`. + X : numpy.ndarray + States at the time steps `T`. + Y : numpy.ndarray + Output values at the time steps `T`. + """ + delta = T[1] - T[0] + + M = (self.E - delta / 2 * (self.J - self.R)) + A = (self.E + delta / 2 * (self.J - self.R)) + + X = np.zeros((self.n, len(T))) + X[:, 0] = x0 + + for i in range(len(T) - 1): + U_midpoint = 1 / 2 * (U[:, i] + U[:, i + 1]) + X[:, i + 1] = np.linalg.solve(M, A @ X[:, i] + delta * (self.G - self.P) @ U_midpoint) + + Y = (self.G + self.P).T @ X + self.D @ U + + return U, X, Y + + def to_lti(self, matrices=True): + """ + Converts the port-Hamiltonian linear time-invariant system in a standard linear time-invariant system. + + Parameters + ---------- + matrices : bool, optional + If `True` the lti matrices are returned, else a pymor `LTIModel`. Default `True`. + + Returns + ------- + A : numpy.ndarray + `A` matrix of the lti system. + B : numpy.ndarray + `B` matrix of the lti system. + C : numpy.ndarray + `C` matrix of the lti system. + D : numpy.ndarray + `D` matrix of the lti system. + E : numpy.ndarray + `E` matrix of the lti system. + lti : `LTIModel` + Instance of pyMOR `LTIModel`. + """ + A = self.J - self.R + B = self.G - self.P + C = (self.G + self.P).T + D = self.D + E = self.E + + if matrices: + return A, B, C, D, E + else: + return LTIModel.from_matrices(A, B, C, D, E=E) + + +if __name__ == "__main__": + from model.msd import msd + E, J, R, G, P, D, Q = msd() + n = J.shape[0] + x0 = np.zeros(n) + u = lambda t: np.array([2 * np.sin(t ** 2), -np.cos(t ** 2)]) + T = np.linspace(0, 6, 500) + ph = PHLTI(E, J, R, G, P, D, Q) + + U, X, Y = ph.sim(u, T, x0) diff --git a/src/utils/__init__.py b/src/utils/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/utils/norm.py b/src/utils/norm.py new file mode 100644 index 0000000..3f0d31c --- /dev/null +++ b/src/utils/norm.py @@ -0,0 +1,41 @@ +import numpy as np + + +def lnorm(X, p=2, tau=1.0): + r""" + Approximates the :math:`\mathcal{L}^p` norm for discrete values. + + Parameters + ---------- + X : numpy.ndarray + Function values. + p : int, optional + Order. Default 2. + tau : float, optional + Stepsize. Default 1.0. + + Returns + ------- + norm : float + Norm. + """ + if p == np.inf: + return X.max() + + w = np.ones(X.shape[1]) * tau / 2 + w[0] = tau + w[-1] = tau + x = np.zeros(X.shape[1]) + for x_i in X: + x += x_i ** p + + return np.inner(w, x) ** (1 / p) + + +if __name__ == '__main__': + x = np.vstack((np.arange(1, 11), np.arange(10, 0, -1))) + + l2 = lnorm(x, p=2) + print(f'|f(x)|_l2={l2:.2e}') + linf = lnorm(x, p=np.inf) + print(f'|f(x)|_linf={linf:.2e}') diff --git a/src/utils/plotting.py b/src/utils/plotting.py new file mode 100644 index 0000000..720b2d0 --- /dev/null +++ b/src/utils/plotting.py @@ -0,0 +1,313 @@ +import os + +import matplotlib.pyplot as plt +import numpy as np + +import config + + +def fig_size(width_pt, fraction=1, ratio=(5 ** .5 - 1) / 2, subplots=(1, 1)): + """ + Returns the width and heights in inches for a matplotlib figure. + + Parameters + ---------- + width_pt : float + Document width in points, in latex can be determined with `\showthe\textwidth`. + fraction : float, optional + The fraction of the width with which the figure will occupy. Default 1. + ratio : float, optional + Ratio of the figure. Default is the golden ratio. + subplots : tuple, optional + The shape of subplots. + + Returns + ------- + fig_width_in : float + Width of the figure in inches. + fig_height_in : float + Height of the figure in inches. + """ + # Width of figure (in pts) + fig_width_pt = width_pt * fraction + # Convert from pt to inches + inches_per_pt = 1 / 72.27 + + # Golden ratio to set aesthetic figure height + golden_ratio = ratio + + # Figure width in inches + fig_width_in = fig_width_pt * inches_per_pt + # Figure height in inches + fig_height_in = fig_width_in * golden_ratio * (subplots[0] / subplots[1]) + + return fig_width_in, fig_height_in + + +def new_fig(width_pt=420, fraction=1 / 2, ratio=(5 ** .5 - 1) / 2, subplots=(1, 1)): + """ + Creates new instance of a `matplotlib.pyplot.figure` by using the `fig_size` function. + + Parameters + ---------- + width_pt : float, optional + Document width in points, in latex can be determined with `\showthe\textwidth`. + Default is 420. + fraction : float, optional + The fraction of the width with which the figure will occupy. Default 1. + ratio : float, optional + Ratio of the figure. Default is the golden ratio. + subplots : tuple, optional + The shape of subplots. + + Returns + ------- + fig : matplotlib.pyplot.figure + The figure. + """ + + fig = plt.figure(figsize=fig_size(width_pt, fraction, ratio, subplots)) + return fig + + +def trajectories_plot(T, Y, Y_dmd=None, zoom=None, label='y', train=False): + """ + Creates plot of a trajectory (input, state, output, ...). + + Parameters + ---------- + T : numpy.ndarray + Time steps. + Y : numpy.ndarray + Trajectory value at the time steps. + Y_dmd : numpy.ndarray, optional + Approximated trajectory value at the time steps. Default `None`. + zoom : int, optional + Number of samples to plot. Default `None`. + label : str, optional + Label of the trajectory. Default 'y'. + train : bool, optional + Flag if trajectory is related to the training. + """ + fraction = 1 / 2 if config.save_results else 1 + + for i, y in enumerate(Y): + label_ = f'{label}_{i + 1}' if Y.shape[0] > 1 else f'{label}' + title = f'Training ${label_}$' if train else f'Testing ${label_}$' + match label: + case 'y': + ylabel = 'Output' + case 'u': + ylabel = 'Input' + case _: + ylabel = '' + + fig = new_fig(fraction=fraction) + ax = fig.add_subplot(1, 1, 1) + ax.set_title(title) + if zoom is not None: + T = T[:zoom] + ax.set(xlim=[np.min(T), np.max(T)]) + ax.set(xlabel='Time (s)', ylabel=ylabel) + if zoom is not None: + ax.plot(T, Y[i, :zoom], label=f'${label_}$') + else: + ax.plot(T, Y[i], label=f'${label_}$') + + if Y_dmd is not None: + dmd_label = r'\widetilde{' + label + '}' + dmd_ylabel = f'{dmd_label}_{i + 1}' if Y.shape[0] > 1 else f'{dmd_label}' + if zoom is not None: + ax.plot(T, Y_dmd[i,:zoom], ls='--', label=f'${dmd_ylabel}$') + else: + ax.plot(T, Y_dmd[i], ls='--', label=f'${dmd_ylabel}$') + + ax.legend(loc='best') + + if config.save_results: + assert config.figures_path is not None and config.exp_id is not None + ax.set_title('') + fig.tight_layout() + filename = f'{config.exp_id}_{label_}_train' if train else f'{config.exp_id}_{label_}' + if zoom is not None: + filename += '_zoom' + filename += '.pgf' + fig.savefig(os.path.join(config.figures_path, filename)) + else: + plt.show() + + +def abs_error_plot(T, Y, Y_dmd, label='y', zoom=None): + """ + Creates plot of an absolute error between two trajectories (state, output, ...). + + Parameters + ---------- + T : numpy.ndarray + Time steps. + Y : numpy.ndarray + Trajectory value at the time steps. + Y_dmd : numpy.ndarray, optional + Approximated trajectory value at the time steps. Default `None`. + label : str, optional + Label of the trajectory. Default 'y'. + zoom : int, optional + Number of samples to plot. Default `None`. + """ + + assert Y.shape == Y_dmd.shape + assert len(T) == Y.shape[1] + + fraction = 1 / 2 if config.save_results else 1 + + for i, y in enumerate(Y): + if zoom is not None: + T = T[:zoom] + + fig = new_fig(fraction=fraction) + ax = fig.add_subplot(1, 1, 1) + ax.set(xlim=[np.min(T), np.max(T)]) + ax.set(xlabel='Time (s)', ylabel=f'Absolut error') + label_ = f'$|{label}_{i + 1} - ' + r'\widetilde{' + label + '}_' + str(i + 1) + '|$' if len(Y) > 1 \ + else f'$|{label} - ' + r'\widetilde{' + label + '}|$' + if zoom is not None: + E = np.abs(Y[i, :zoom] - Y_dmd[i, :zoom]) + else: + E = np.abs(Y[i] - Y_dmd[i]) + + ax.semilogy(T, E, label=label_) + ax.legend() + + if config.save_results: + ax.set_title('') + fig.tight_layout() + filename = f'{config.exp_id}_error_{label}_{i + 1}' + if zoom is not None: + filename += '_zoom' + filename += '.pgf' + fig.savefig(os.path.join(config.figures_path, filename)) + else: + plt.show() + + +def bode_plot(w, lti, lti_dmd, i=None, j=None): + """ + Creates bode plot of two lti systems. + + Parameters + ---------- + w : numpy.ndarray + Frequencies. + lti : pymor.models.iosys.LTIModel + Original lti system. + lti_dmd : pymor.models.iosys.LTIModel + Identified lti system. + i : int, optional + If specified only the bode plot of the `i`,`j` input output combination is plotted. Default `None`. + j : int, optional + If specified only the bode plot of the `i`,`j` input output combination is plotted. Default `None`. + """ + fig = new_fig(fraction=1, ratio=1) + ax = fig.subplots(2 * config.m, config.m, sharex=True, squeeze=False) + artists = lti.transfer_function.bode_plot(w, ax=ax, label='G(s)') + artists_dmd = lti_dmd.transfer_function.bode_plot(w, ax=ax, label='\widetilde{G}(s)', ls='--') + if not config.save_results: + plt.show() # can only show figures when Matplotlib is not using pgf + + if i is not None and j is not None: + # Bode of i,j input pair + + # Magnitude + fig = new_fig() + ax = fig.add_subplot(111) + ax.set(xlim=[np.min(w), np.max(w)]) + ax.loglog(*artists[2 * i, j][0].get_data(), label='$G$') + ax.loglog(*artists_dmd[2 * i, j][0].get_data(), + ls='--', label=r'$\widetilde{G}$') + ax.legend() + ax.set_title(r'Magnitude plot of $G_{' + str(i + 1) + str(j + 1) + '}$') + ax.set_xlabel('Frequency (rad/s)') + ax.set_ylabel('Magnitude') + + if config.save_results: + ax.set_title('') + fig.tight_layout() + fig.savefig(os.path.join(config.figures_path, f'{config.exp_id}_mag_{i + 1}_{j + 1}.pgf')) + else: + plt.show() + + # Phase + fig = new_fig() + ax = fig.add_subplot(111) + ax.set(xlim=[np.min(w), np.max(w)]) + ax.semilogx(*artists[2 * i + 1, j][0].get_data(), label='$G$') + ax.semilogx(*artists_dmd[2 * i + 1, j][0].get_data(), + ls='--', label=r'$\widetilde{G}$') + ax.legend() + ax.set_title(r'Phase plot of $G_{' + str(i + 1) + str(j + 1) + '}$') + ax.set_xlabel('Frequency (rad/s)') + ax.set_ylabel('Phase (deg)') + + if config.save_results: + ax.set_title('') + fig.tight_layout() + fig.savefig(os.path.join(config.figures_path, f'{config.exp_id}_phase_{i + 1}_{j + 1}.pgf')) + else: + plt.show() + + +def magnitude_plot(w, lti): + """ + Creates magnitude bode plot an lti system. + + Parameters + ---------- + w : numpy.ndarray + Frequencies. + lti : pymor.models.iosys.LTIModel + (Error) lti system. + """ + fraction = 1 / 2 if config.save_results else 1 + + fig = new_fig(fraction=fraction) + ax = fig.add_subplot(1, 1, 1) + ax.set(xlim=[np.min(w), np.max(w)]) + lti.transfer_function.mag_plot(w, ax=ax, label='Error') + _ = ax.legend() + if config.save_results: + ax.set_title('') + fig.tight_layout() + fig.savefig(os.path.join(config.figures_path, f'{config.exp_id}_error_mag_plot.pgf')) + else: + plt.show() + + +def poles_plot(lti, lti_dmd=None): + """ + Creates poles bode plot of one or two lti systems. + + Parameters + ---------- + lti : pymor.models.iosys.LTIModel + Original lti system. + lti_dmd : pymor.models.iosys.LTIModel + Identified lti system. + """ + # Eigenvalues + poles = lti.poles() + fig, ax = plt.subplots() + ax.plot(poles.real, poles.imag, '.', label='FOM') + if lti_dmd is not None: + poles_dmd = lti_dmd.poles() + ax.plot(poles_dmd.real, poles_dmd.imag, 'x', label='dmd') + + _ = ax.set_title('Poles') + _ = ax.legend() + ax.set(xlabel='Re($\lambda$)', ylabel='Im($\lambda$)') + + if config.save_results: + ax.set_title('') + fig.tight_layout() + fig.savefig(os.path.join(config.figures_path, f'{config.exp_id}_poles.pgf')) + else: + plt.show()