A Python/NumPy implementation of the Loewner framework for data-driven rational interpolation of complex-valued data and dynamical systems.
Jump to a section
Background (Explanation & Some Theory)
Getting Started With This Module
It doesn't take much to start building interpolants.
import numpy as np
import loewner_framework as lf
import loewner_framework.linear_daes as ld
Lambda = np.diag([0.0, complex(imag=1.0), complex(imag=-1.0)])
R = np.array([[1, 1, 1]])
W = np.array([[complex(real=5.0), complex(real=0.5, imag=-0.1), complex(real=0.5, imag=0.1)]])
rtd = lf.RightTangentialData(Lambda, R, W)
M = np.diag([1.0, complex(imag=10.0), complex(imag=-10.0)])
L = np.array([ [1],
[1],
[1]])
V = np.array([ [1],
[complex(real=1.0, imag=-0.5)],
[complex(real=1.0, imag=0.5)]])
ltd = lf.LeftTangentialData(M, L, V)
loewnertd = lf.LoewnerTangentialData(rtd, ltd)
interpolant_builder = lf.InterpolantFactory(loewnertd)
interpolant = interpolant_builder.minimal_order(label="Interpolant")
And you can easily view the results on a Bode plot.
bode_plot = ld.BodePlot()
bode_plot.add_system(interpolant, color='b')
bode_plot.add_data_tick(at_frequency=10**-2, pin_to_system_num=0, pin_to_output_num=1, pin_to_input_num=1, markersize=12, color="red")
bode_plot.add_data_tick(at_frequency=10**0, pin_to_system_num=0, pin_to_output_num=1, pin_to_input_num=1, markersize=12, color="red")
bode_plot.add_data_tick(at_frequency=10**1, pin_to_system_num=0, pin_to_output_num=1, pin_to_input_num=1, markersize=12, color="red")
bode_plot.show(w_start=-2, w_end=2, w_num_points=10000)
The "interpolant" is returned as a generalized state-space / DAE model in the form of a LinearDAE object, with matrices
interpolant.E, interpolant.A, interpolant.B, interpolant.C, interpolant.D
The actual rational function that interpolates the right and left tangential data sets is the transfer function of the LinearDAE object,
This transfer function can be evaluated directly, e.g.
interpolant.tf(complex(real=0.0, imag=1.0))
The result is a numpy.ndarray of dimension 2. This is maintained even for SISO systems with a scalar transfer function for the purposes of compatibility with MIMO systems having matrix transfer functions.
See a mistake? Have a feature request? Want to contribute?
Submit an issue, or reach out to joeldsimard@gmail.com.
The background exposition in this section has been paraphrased from the PhD dissertation J. Simard, Interpolation and model reduction of nonlinear systems in the Loewner framework. Diss. Imperial College London, 2023.
For further details and specific references, see the Some Sources section. You can also check out the MOR Wiki which covers some basic concepts, however this page is somewhat out of date.
The objective of the rational interpolation problem for linear time-invariant (LTI) systems is to construct an internal or external description of a system the transfer matrix of which is consistent with finite sets of tangential data in the complex plane, and the generalized realization problem for LTI systems entails the construction of a minimal state-space model of a system given the transfer matrix of the system.
The Loewner framework was developed as a toolset used for the solution of the rational interpolation problem and the generalized realization problem for LTI DAEs.
The core object of the Loewner framework is the Loewner matrix, a divided-difference matrix constructed from two sets of points in the complex plane: the right tangential data and the left tangential data. The sets of tangential data are obtained by sampling, along particular directions, a rational matrix function evaluated at particular points in the complex plane. The Loewner matrix has a fine structure that allows its factorization into the tangential generalized controllability and observability matrices. Once obtained, the tangential generalized controllability and observability matrices can be used to construct LTI descriptor system realizations interpolating the sets of tangential data.
While the Loewner framework can be interpreted entirely from the complex rational interpolation point of view, it is both convenient and useful to consider the dynamical systems perspective. In particular, we deal with LTI DAEs of the form
with generalized state
The purpose of the Loewner framework is to determine a DAE of this form such that its transfer function interpolates sets of right and left tangential interpolation data. Therefore, if one just wants to get a rational function that interpolates sets of complex data we can simply build the DAE in the Loewner framework and then retrieve its transfer function.
Tangential data consist of two disjoint sets of data: the right tangential interpolation data, and the left tangential interpolation data. The right tangential interpolation data are represented by the set
and the left tangential interpolation data are represented by the set
These data sets act as the interpolation constraints for the DAE that will be produced by the Loewner framework. Basically, the transfer function of the DAE will simultaneously satisfy
and
For the systems perspective we represent the data via right tangential data matrices
and the left tangential data matrices
To ensure the existence of a unique Loewner matrix associated to the data matrices, it is required that the matrices
You can also construct real-valued equivalent data matrices for the purpose of building rational interpolants with real-valued coefficients, however this is not yet discussed here.
The primary tools used to accomplish the interpolation objective in the Loewner framework are the Loewner matrix,
for rows
The matrices are also obtained as the solutions to Sylvester equations,
which is the preferred method of construction when using low-dimensional equivalent real-valued tangential data matrices.
These matrices satisfy
There are some additional matrices and relationships associated to the systems perspective, such as the tangential generalized controllability and observability matrices, however these are not yet discussed here. The module builds and utilizes these objects, however it is not required for the user to deal with them.
If the number of right tangential interpolation points is equal to the number of left tangential interpolation points, i.e.
which has the transfer function
This is not very flexible as there are no degrees of freedom in the model to further influence its properties beyond achieving interpolation.
A feedforward term can be leveraged to provide a family of differential-algebraic interpolants matching the tangential data. So long as it is regular, then for any choice of matrix
with transfer function
interpolates the tangential interpolation data. Since the matrix
The most general scenario for interpolant construction also considers Loewner matrices that don't have to be square, and interpolants that are allowed to have dimension greater than or equal to
with transfer function
is an interpolant of the tangential interpolation data.
For any desired system dimension, the matrices
Furthermore, for any dimension greater than or equal to
This section describes the details about some of the options provided for building interpolants in this module. The subsections are presented in order based upon the actual construction of an interpolant, with later subsections depending on those that come before.
To set up the tangential interpolation data sets, we instantiate RightTangentialData and LeftTangentialData objects. The constructors of these classes must be provided the data points in some form. The matrices
To build the right and left data sets directly from frequency data we can prepare lists of complex points for the right tangential data set
lambda_f = [0.1*j, 10.0*j, 1000.0*j]
r_f = [[1, 1, 1]]
w_f = [[2, -1, 0.25]]
and for the left tangential data set
mu_f = [0.0, 1.0*j, 100.0*j]
ell_f = [ [1],
[1],
[1]]
v_f = [ [-0.1],
[-2],
[0.5]]
These lists are used to prepare numpy.ndarrays with dimension two and then provided to the RightTangentialData and LeftTangentialData class constructors. This is done via
Lambda = np.diag(lambda_f)
R = np.array(r_f)
W = np.array(w_f)
rtd = lf.RightTangentialData(Lambda, R, W)
and
M = np.diag(mu_f)
L = np.array(ell_f)
V = np.array(v_f)
ltd = lf.LeftTangentialData(M, L, V)
To use equivalent real-valued data matrices we can directly define the numpy.ndarrays via
Lambda = np.array([ [0, 0, 0],
[0, 0, 1],
[0, -1, 0]])
R = np.array([[1, 1, 0]])
W = np.array([[1, 2, 1]])
rtd = lf.RightTangentialData(Lambda, R, W)
and
M = np.array([ [1, 0, 0],
[0, 0, -100],
[0, 100, 0]])
L = np.array([ [1],
[0],
[1]])
V = np.array([ [-2],
[-1],
[-1]])
ltd = lf.LeftTangentialData(M, L, V)
Instead of providing the tangential data matrices
n, m, p = 10, 1, 1
E = np.eye(n)
A = np.random.random_sample((n,n))
B = np.random.random_sample((n,m))
C = np.random.random_sample((p,n))
D = np.random.random_sample((p,m))
system = ld.LinearDAE(A, B, C, D, label="System")
We set the frequencies and tangential directions to evaluate this system's transfer function at (in equivalent real-valued matrix form) via
Lambda = np.array([ [0, 0, 0],
[0, 0, 1],
[0, -1, 0]])
R = np.array([[1, 1, 0]])
and
M = np.array([ [1, 0, 0],
[0, 0, -100],
[0, 100, 0]])
L = np.array([ [1],
[0],
[1]])
Finally, we provide the LinearDAE object to the constructors in place of the
rtd = lf.RightTangentialData(Lambda, R, system)
ltd = lf.LeftTangentialData(M, L, system)
When a system is used to generate the data, you can retrieve the associated tangential generalized controllability and observability matrices
rtd.X
ltd.Y
If
To check that the RightTangentialData and LeftTangentialData objects contain data this is useable for interpolation, check the isComplete property
rtd.isComplete
>>> True (if data is valid), False (if data is not valid)
ltd.isComplete
>>> True (if data is valid), False (if data is not valid)
You can retrieve all the data matrices via system properties
rtd.Lambda, rtd.R, rtd.W
ltd.M, ltd.L, ltd.V
To guarantee that the Loewner matrices can be generated, we need that both data sets are valid (isComplete == True). In addition we need that
To see if the RightTangentialData and LeftTangentialData objects were generated using an underlying system, check the isFromSystem property
rtd.isFromSystem
>>> True (if from LinearDAE), False (if from matrix W)
ltd.isFromSystem
>>> True (if from LinearDAE), False (if from matrix W)
This ensures that the
To get the Loewner and shifted Loewner matrices we instantiate a LoewnerTangentialData object. The constructor must be provided with valid RightTangentialData and LeftTangentialData objects (isComplete being True for both). The associated
loewnertd = lf.LoewnerTangentialData(rtd, ltd)
The isFromSystem property is True only if both the RightTangentialData and LeftTangentialData isFromSystem properties are True
loewnertd.isFromSystem
>>> True (if from LinearDAE), False (if from matrix W)
To ensure that unique Loewner and shifted Loewner matrices were successfully determined, check the isComplete property.
loewnertd.isComplete
>>> True (if data is valid), False (if data is not valid)
If this property is True, then we can use this data to start building interpolants.
You can retrieve the Loewner and shifted Loewner matrices via properties
loewnertd.Loewner, loewnertd.shiftedLoewner
To build interpolants, initialize an InterpolantFactory object. The constructor must be provided a LoewnerTangentialData object with the isComplete property being True
interpolant_builder = lf.InterpolantFactory(loewnertd)
If the isComplete property of the LoewnerTangentialData object is True, then all should be good. But this can be further checked via the InterpolantFactory object
interpolant_builder.isComplete
If the Loewner matrix is square we can build interpolants of minimal order without feedforward term
interp1 = interpolant_builder.minimal_order(label="Interpolant Without Feedforward")
We can also provide a feedforward term to get different models of minimal dimension
D = np.random.random_sample((interpolant_builder.p, interpolant_builder.m))
interp2 = interpolant_builder.minimal_order(D, label="Interpolant With Feedforward")
The returned objects are LinearDAE objects, so their state-space matrices can be retrieved from properties, and the interpolating transfer function can be constructed or evaluated
interp1.E, interp1.A, interp1.B, interp1.C, interp1.D
interp1.tf(complex(real=0.0, imag=1.0))
To generate interpolants more generally (with non-square Loewner matrices and/or non-minimal dimension) we must provided free parameters
We can get a dictionary containing the parameter shapes required to produce an interpolant of a desired dimension using the InterpolantFactory's parameter_dimensions method
desired_dimension = max(interpolant_builder.nu, interpolant_builder.rho) + 1
dimension_possible, shapes_dict: tuple[bool,dict] = interpolant_builder.parameter_dimensions(desired_dimension)
Due to returned dictionary format, if the desired interpolant dimension is possible to achieve then we can generate random free parameter matrices of exactly the correct sizes (or set them to None if not usable) using the following idiom:
if dimension_possible:
params_dict = { "D": np.random.random_sample(shapes_dict["D"]) if shapes_dict["D"] else None,
"P": np.random.random_sample(shapes_dict["P"]) if shapes_dict["P"] else None,
"Q": np.random.random_sample(shapes_dict["Q"]) if shapes_dict["Q"] else None,
"G": np.random.random_sample(shapes_dict["G"]) if shapes_dict["G"] else None,
"T": np.random.random_sample(shapes_dict["T"]) if shapes_dict["T"] else None,
"H": np.random.random_sample(shapes_dict["H"]) if shapes_dict["H"] else None,
"F": np.random.random_sample(shapes_dict["F"]) if shapes_dict["F"] else None}
We can also check that some set of free parameters can be used to generate an interpolant using the InterpolantFactory's check_consistent_shapes method
parameters_consistent, associated_dimension: tuple[bool,int] = interpolant_builder.check_consistent_shapes(**params_dict)
This also determines the dimension of the interpolant associated to the provided free parameters, if the interpolant is possible.
Note that if everything is as expected, then we should also have
associated_dimension == desired_dimension
>>> True
The reason this may not be be true is because when generating an interpolant of minimal-order with the parameterization method for a square Loewner matrix it is possible to pass a valid parameter dictionary where each value in params_dict is equal to None (and the D matrix is interpreted to be zeros), however the parameter_dimensions method (used beforehand) will return a shapes dictionary with all values equal to None when the desired dimension is not possible. Whenever a desired dimension is possible, the value in the shapes dictionary associated to the D key will not be None. Hence, to avoid this issue we can check one of the following conditionals:
associated_dimension == desired_dimension
or
not shapes_dict["D"] is None
or
dimension_possible and parameters_consistent
Finally, once we are sure that the free parameters
if dimension_possible and parameters_consistent:
interp: ld.LinearDAE = interpolant_builder.parameterization(**params_dict, label="Interpolant")
Just like using the minimal_order method, we can retrieve the state-space matrices and reconstruct or evaluate the interpolation rational transfer function
interp.E, interp.A, interp.B, interp.C, interp.D
interp.tf(complex(real=0.0, imag=1.0))
For a square Loewner matrix, interpolants with dimension
desired_poles = np.arange(-6, 0, 1) * 1
poleplaced_interp = interpolant_builder.double_order_pole_placed(desired_poles, D=None, label="Pole-Placed Interpolant")
The free parameters
This module is built on top of the modules generalized-sylvester and py-linear-DAEs. There is no need to import them, they are built-in; however the documentation for these modules may be useful for using the LinearDAE class, the BodePlot class, and the generalized_sylvester.solve function.
Currently the Loewner objects are determined via solving generalized Sylvester equations. This allows for the easy construction of real-valued interpolants as the equivalent real-valued matrices can be used directly in the generalized Sylvester equation approach. However, using real-valued data matrices directly means that we can't build the Loewner matrices entry-wise; this would require transforming to a "diagonalized" classical complex data representation, which can be done but is not yet implemented here. The generalized Sylvester equation approach doesn't scale as well as the entry-wise approach as the quantity of data increases, so for very large amounts of data or very large systems the current approach may be slow or use a large amount of memory. The entry-wise approach will be added eventually, along with functionality for converting between classical complex data and real-valued matrix data.
In scenarios where there is "redundant data" the minimal order interpolant may actually have a dimension less than max(rho, nu), in which case an SVD procedure can be used to extract a lower order interpolant. This procedure is not implemented here, and for now the minimal order interpolant is considered to have dimension max(rho, nu) with the assumption that there is no redundant data, or that we simply don't care that we have redundant data.
loewner_framework.RightTangentialData.__init__(self, Lambda: numpy.ndarray, R: numpy.ndarray, W: numpy.ndarray | linear_daes.LinearDAE | None = None)
RightTangentialData.isComplete
RightTangentialData.isFromSystem
RightTangentialData.rho
RightTangentialData.m
RightTangentialData.p
RightTangentialData.Lambda
RightTangentialData.R
RightTangentialData.W
RightTangentialData.X
RightTangentialData.__repr__(self) -> str
RightTangentialData.__bool__(self) -> bool
loewner_framework.LeftTangentialData.__init__(self, M: numpy.ndarray, L: numpy.ndarray, V: numpy.ndarray | linear_daes.LinearDAE | None = None)
LeftTangentialData.isComplete
LeftTangentialData.isFromSystem
LeftTangentialData.nu
LeftTangentialData.p
LeftTangentialData.m
LeftTangentialData.M
LeftTangentialData.L
LeftTangentialData.V
LeftTangentialData.Y
LeftTangentialData.__repr__(self) -> str
LeftTangentialData.__bool__(self) -> bool
loewner_framework.LoewnerTangentialData.__init__(self, rtd: loewner_framework.RightTangentialData, ltd: loewner_framework.LeftTangentialData)
LoewnerTangentialData.isComplete
LoewnerTangentialData.isFromSystem
LoewnerTangentialData.rtd
LoewnerTangentialData.ltd
LoewnerTangentialData.Loewner
LoewnerTangentialData.shiftedLoewner
LoewnerTangentialData.rho
LoewnerTangentialData.nu
LoewnerTangentialData.m
LoewnerTangentialData.p
LoewnerTangentialData.Lambda
LoewnerTangentialData.R
LoewnerTangentialData.W
LoewnerTangentialData.M
LoewnerTangentialData.L
LoewnerTangentialData.V
LoewnerTangentialData.__repr__(self) -> str
LoewnerTangentialData.__bool__(self) -> bool
loewner_framework.InterpolantFactory.__init__(self, loewnertd : loewner_framework.LoewnerTangentialData)
InterpolantFactory.isComplete
InterpolantFactory.loewnertd
InterpolantFactory.rho
InterpolantFactory.nu
InterpolantFactory.m
InterpolantFactory.p
InterpolantFactory.Lambda
InterpolantFactory.R
InterpolantFactory.W
InterpolantFactory.M
InterpolantFactory.L
InterpolantFactory.V
InterpolantFactory.Loewner
InterpolantFactory.shiftedLoewner
InterpolantFactory.minimal_order(self, D: numpy.ndarray | None = None, label: str = "") -> linear_daes.LinearDAE | None
InterpolantFactory.parameter_dimensions(self, total_dimension: int) -> tuple[bool,dict]
InterpolantFactory.check_consistent_shapes(self, D: numpy.ndarray | None = None, P: numpy.ndarray | None = None,
Q: numpy.ndarray | None = None, G: numpy.ndarray | None = None,
T: numpy.ndarray | None = None, H: numpy.ndarray | None = None,
F: numpy.ndarray | None = None) -> tuple[bool,int]
InterpolantFactory.parameterization(self, D: numpy.ndarray | None = None, P: numpy.ndarray | None = None,
Q: numpy.ndarray | None = None, G: numpy.ndarray | None = None,
T: numpy.ndarray | None = None, H: numpy.ndarray | None = None,
F: numpy.ndarray | None = None, label: str = "") -> linear_daes.LinearDAE
InterpolantFactory.double_order_pole_placed(self, desired_poles: _np.ndarray, D: _np.ndarray | None = None, label: str = "") -> _ld.LinearDAE
InterpolantFactory.__repr__(self) -> str
InterpolantFactory.__bool__(self) -> bool