pymor.models.iosys

Module Contents

class pymor.models.iosys.BilinearModel(A, N, B, C, D, E=None, sampling_time=0, error_estimator=None, visualizer=None, name=None)

Bases: pymor.models.interface.Model

Class for bilinear systems.

This class describes input-output systems given by

\[\begin{split}E x'(t) & = A x(t) + \sum_{i = 1}^m{N_i x(t) u_i(t)} + B u(t), \\ y(t) & = C x(t) + D u(t),\end{split}\]

if continuous-time, or

\[\begin{split}E x(k + 1) & = A x(k) + \sum_{i = 1}^m{N_i x(k) u_i(k)} + B u(k), \\ y(k) & = C x(k) + D u(t),\end{split}\]

if discrete-time, where \(E\), \(A\), \(N_i\), \(B\), \(C\), and \(D\) are linear operators and \(m\) is the number of inputs.

Parameters:

• A – The Operator A.

• N – The tuple of Operators N_i.

• B – The Operator B.

• C – The Operator C.

• D – The Operator D or None (then D is assumed to be zero).

• E – The Operator E or None (then E is assumed to be identity).

• sampling_time – 0 if the system is continuous-time, otherwise a positive number that denotes the sampling time (in seconds).

• error_estimator – An error estimator for the problem. This can be any object with an estimate_error(U, mu, model) method. If error_estimator is not None, an estimate_error(U, mu) method is added to the model which will call error_estimator.estimate_error(U, mu, self).

• visualizer – A visualizer for the problem. This can be any object with a visualize(U, model, ...) method. If visualizer is not None, a visualize(U, *args, **kwargs) method is added to the model which forwards its arguments to the visualizer’s visualize method.

• name – Name of the system.

order

The order of the system (equal to A.source.dim).

dim_input

The number of inputs.

dim_output

The number of outputs.

A

The Operator A.

N

The tuple of Operators N_i.

B

The Operator B.

C

The Operator C.

D

The Operator D.

E

The Operator E.

cache_region = 'memory'

class pymor.models.iosys.ImpulseFunction(dim_input, component, sampling_time)

Bases: pymor.analyticalproblems.functions.Function

Interface for Parameter dependent analytical functions.

Every Function is a map of the form

f(μ): Ω ⊆ R^d -> R^(shape_range)

The returned values are NumPy arrays of arbitrary (but fixed) shape. Note that NumPy distinguishes between one-dimensional arrays of length 1 (with shape (1,)) and zero-dimensional scalar arrays (with shape ()). In pyMOR, we usually expect scalar-valued functions to have shape_range == ().

While the function might raise an error if it is evaluated for an argument not in the domain Ω, the exact behavior is left undefined.

Functions are vectorized in the sense, that if x.ndim == k, then

f(x, μ)[i0, i1, ..., i(k-2)] == f(x[i0, i1, ..., i(k-2)], μ).

In particular, f(x, μ).shape == x.shape[:-1] + shape_range.

dim_domain

The dimension d > 0.

shape_range

The shape of the function values.

Methods

evaluate

Evaluate the function for given argument x and parameter values mu.

evaluate(x, mu=None)

Evaluate the function for given argument x and parameter values mu.

class pymor.models.iosys.LTIModel(A, B, C, D=None, E=None, sampling_time=0, T=None, initial_data=None, time_stepper=None, num_values=None, presets=None, solver_options=None, ast_pole_data=None, error_estimator=None, visualizer=None, name=None)

Bases: pymor.models.interface.Model

Class for linear time-invariant systems.

This class describes input-state-output systems given by

\[\begin{split}E(\mu) \dot{x}(t, \mu) & = A(\mu) x(t, \mu) + B(\mu) u(t), \\ y(t, \mu) & = C(\mu) x(t, \mu) + D(\mu) u(t),\end{split}\]

if continuous-time, or

\[\begin{split}E(\mu) x(k + 1, \mu) & = A(\mu) x(k, \mu) + B(\mu) u(k), \\ y(k, \mu) & = C(\mu) x(k, \mu) + D(\mu) u(k),\end{split}\]

if discrete-time, where \(A\), \(B\), \(C\), \(D\), and \(E\) are linear operators.

All methods related to the transfer function (e.g., frequency response calculation and Bode plots) are attached to the transfer_function attribute.

Parameters:

• A – The Operator A.

• B – The Operator B.

• C – The Operator C.

• D – The Operator D or None (then D is assumed to be zero).

• E – The Operator E or None (then E is assumed to be identity).

• sampling_time – 0 if the system is continuous-time, otherwise a positive number that denotes the sampling time (in seconds).

• T – The final time T.

• initial_data – The initial data x_0. Either a VectorArray of length 1 or (for the Parameter-dependent case) a vector-like Operator (i.e. a linear Operator with source.dim == 1) which applied to NumpyVectorArray(np.array([1])) will yield the initial data for given parameter values. If None, it is assumed to be zero.

• time_stepper – The time-stepper to be used by solve.

• num_values – The number of returned vectors of the solution trajectory. If None, each intermediate vector that is calculated is returned.

• presets – A dict of preset attributes or None. The dict must only contain keys that correspond to attributes of LTIModel such as poles, c_lrcf, o_lrcf, c_dense, o_dense, hsv, h2_norm, hinf_norm, l2_norm and linf_norm. Additionally, the frequency at which the \(\mathcal{H}_\infty/\mathcal{L}_\infty\) norm is attained can be preset with fpeak.

• solver_options – The solver options to use to solve matrix equations.

• ast_pole_data –

  Used in get_ast_spectrum. Can be:

  • dictionary of parameters for eigs,

  • list of anti-stable eigenvalues (scalars),

  • tuple (lev, ew, rev) where ew contains the sorted anti-stable eigenvalues

    and lev and rev are VectorArrays representing the eigenvectors.

  • None if anti-stable eigenvalues should be computed via dense methods.

• error_estimator – An error estimator for the problem. This can be any object with an estimate_error(U, mu, model) method. If error_estimator is not None, an estimate_error(U, mu) method is added to the model which will call error_estimator.estimate_error(U, mu, self).

• visualizer – A visualizer for the problem. This can be any object with a visualize(U, model, ...) method. If visualizer is not None, a visualize(U, *args, **kwargs) method is added to the model which forwards its arguments to the visualizer’s visualize method.

• name – Name of the system.

order

The order of the system.

dim_input

The number of inputs.

dim_output

The number of outputs.

A

The Operator A.

B

The Operator B.

C

The Operator C.

D

The Operator D.

E

The Operator E.

transfer_function

The transfer function.

Methods

from_abcde_files	Create LTIModel from matrices stored in .[ABCDE] files.
from_files	Create LTIModel from matrices stored in separate files.
from_mat_file	Create LTIModel from matrices stored in a .mat file.
from_matrices	Create LTIModel from matrices.
get_ast_spectrum	Compute anti-stable subset of the poles of the LTIModel.
gramian	Compute a Gramian.
h2_norm	Compute the \(\mathcal{H}_2\)-norm of the LTIModel.
hankel_norm	Compute the Hankel-norm of the LTIModel.
hinf_norm	Compute the \(\mathcal{H}_\infty\)-norm of the LTIModel.
hsv	Hankel singular values.
impulse_resp	Compute impulse response from all inputs.
l2_norm	Compute the \(\mathcal{L}_2\)-norm of the LTIModel.
linf_norm	Compute the \(\mathcal{L}_\infty\)-norm of the LTIModel.
moebius_substitution	Create a transformed LTIModel by applying an arbitrary Moebius transformation.
poles	Compute system poles.
step_resp	Compute step response from all inputs.
to_abcde_files	Save operators as matrices to .[ABCDE] files in Matrix Market format.
to_abcde_matrices	Return A, B, C, D, and E operators as matrices.
to_continuous	Converts a discrete-time LTIModel to a continuous-time LTIModel.
to_discrete	Converts a continuous-time LTIModel to a discrete-time LTIModel.
to_files	Write operators to files as matrices.
to_mat_file	Save operators as matrices to .mat file.
to_matrices	Return operators as matrices.

cache_region = 'memory'

classmethod from_abcde_files(files_basename, sampling_time=0, T=None, time_stepper=None, num_values=None, presets=None, solver_options=None, error_estimator=None, visualizer=None, name=None)

Create LTIModel from matrices stored in .[ABCDE] files.

Parameters:

• files_basename – The basename of files containing A, B, C, and optionally D and E.

• sampling_time – 0 if the system is continuous-time, otherwise a positive number that denotes the sampling time (in seconds).

• T – The final time T.

• time_stepper – The time-stepper to be used by solve.

• num_values – The number of returned vectors of the solution trajectory. If None, each intermediate vector that is calculated is returned.

• presets – A dict of preset attributes or None. See LTIModel.

• solver_options – The solver options to use to solve the Lyapunov equations.

• error_estimator – An error estimator for the problem. This can be any object with an estimate_error(U, mu, model) method. If error_estimator is not None, an estimate_error(U, mu) method is added to the model which will call error_estimator.estimate_error(U, mu, self).

• visualizer – A visualizer for the problem. This can be any object with a visualize(U, model, ...) method. If visualizer is not None, a visualize(U, *args, **kwargs) method is added to the model which forwards its arguments to the visualizer’s visualize method.

• name – Name of the system.

Returns:

lti – The LTIModel with operators A, B, C, D, and E.

classmethod from_files(A_file, B_file, C_file, D_file=None, E_file=None, sampling_time=0, T=None, initial_data_file=None, time_stepper=None, num_values=None, presets=None, solver_options=None, error_estimator=None, visualizer=None, name=None)

Create LTIModel from matrices stored in separate files.

Parameters:

• A_file – The name of the file (with extension) containing A.

• B_file – The name of the file (with extension) containing B.

• C_file – The name of the file (with extension) containing C.

• D_file – None or the name of the file (with extension) containing D.

• E_file – None or the name of the file (with extension) containing E.

• sampling_time – 0 if the system is continuous-time, otherwise a positive number that denotes the sampling time (in seconds).

• T – The final time T.

• initial_data_file – None or the name of the file (with extension) containing the initial data.

• time_stepper – The time-stepper to be used by solve.

• num_values – The number of returned vectors of the solution trajectory. If None, each intermediate vector that is calculated is returned.

• presets – A dict of preset attributes or None. See LTIModel.

• solver_options – The solver options to use to solve the Lyapunov equations.

• error_estimator – An error estimator for the problem. This can be any object with an estimate_error(U, mu, model) method. If error_estimator is not None, an estimate_error(U, mu) method is added to the model which will call error_estimator.estimate_error(U, mu, self).

• visualizer – A visualizer for the problem. This can be any object with a visualize(U, model, ...) method. If visualizer is not None, a visualize(U, *args, **kwargs) method is added to the model which forwards its arguments to the visualizer’s visualize method.

• name – Name of the system.

Returns:

lti – The LTIModel with operators A, B, C, D, and E.

classmethod from_mat_file(file_name, sampling_time=0, T=None, time_stepper=None, num_values=None, presets=None, solver_options=None, error_estimator=None, visualizer=None, name=None)

Create LTIModel from matrices stored in a .mat file.

Supports the format used in the SLICOT benchmark collection.

Parameters:

• file_name – The name of the .mat file (extension .mat does not need to be included) containing A, B, and optionally C, D, and E.

• sampling_time – 0 if the system is continuous-time, otherwise a positive number that denotes the sampling time (in seconds).

• T – The final time T.

• time_stepper – The time-stepper to be used by solve.

• num_values – The number of returned vectors of the solution trajectory. If None, each intermediate vector that is calculated is returned.

• presets – A dict of preset attributes or None. See LTIModel.

• solver_options – The solver options to use to solve the Lyapunov equations.

• error_estimator – An error estimator for the problem. This can be any object with an estimate_error(U, mu, model) method. If error_estimator is not None, an estimate_error(U, mu) method is added to the model which will call error_estimator.estimate_error(U, mu, self).

• visualizer – A visualizer for the problem. This can be any object with a visualize(U, model, ...) method. If visualizer is not None, a visualize(U, *args, **kwargs) method is added to the model which forwards its arguments to the visualizer’s visualize method.

• name – Name of the system.

Returns:

lti – The LTIModel with operators A, B, C, D, and E.

classmethod from_matrices(A, B, C, D=None, E=None, sampling_time=0, T=None, initial_data=None, time_stepper=None, num_values=None, presets=None, solver_options=None, error_estimator=None, visualizer=None, name=None)

Create LTIModel from matrices.

Parameters:

• A – The NumPy array or SciPy spmatrix A.

• B – The NumPy array or SciPy spmatrix B.

• C – The NumPy array or SciPy spmatrix C.

• D – The NumPy array or SciPy spmatrix D or None (then D is assumed to be zero).

• E – The NumPy array or SciPy spmatrix E or None (then E is assumed to be identity).

• sampling_time – 0 if the system is continuous-time, otherwise a positive number that denotes the sampling time (in seconds).

• T – The final time T.

• initial_data – The initial data x_0 as a NumPy array. If None, it is assumed to be zero.

• time_stepper – The time-stepper to be used by solve.

• num_values – The number of returned vectors of the solution trajectory. If None, each intermediate vector that is calculated is returned.

• presets – A dict of preset attributes or None. See LTIModel.

• solver_options – The solver options to use to solve the Lyapunov equations.

• error_estimator – An error estimator for the problem. This can be any object with an estimate_error(U, mu, model) method. If error_estimator is not None, an estimate_error(U, mu) method is added to the model which will call error_estimator.estimate_error(U, mu, self).

• visualizer – A visualizer for the problem. This can be any object with a visualize(U, model, ...) method. If visualizer is not None, a visualize(U, *args, **kwargs) method is added to the model which forwards its arguments to the visualizer’s visualize method.

• name – Name of the system.

Returns:

lti – The LTIModel with operators A, B, C, D, and E.

get_ast_spectrum(mu=None)

Compute anti-stable subset of the poles of the LTIModel.

Parameters:

mu – Parameter.

Returns:

• lev – VectorArray of left eigenvectors.

• ew – One-dimensional NumPy array of anti-stable eigenvalues sorted from smallest to largest.

• rev – VectorArray of right eigenvectors.

gramian(typ, mu=None)

Compute a Gramian.

Parameters:

• typ –

  The type of the Gramian:

  • 'c_lrcf': low-rank Cholesky factor of the controllability Gramian,

  • 'o_lrcf': low-rank Cholesky factor of the observability Gramian,

  • 'c_dense': dense controllability Gramian,

  • 'o_dense': dense observability Gramian,

  • 'bs_c_lrcf': low-rank Cholesky factor of the Bernoulli stabilized controllability Gramian,

  • 'bs_o_lrcf': low-rank Cholesky factor of the Bernoulli stabilized observability Gramian,

  • 'lqg_c_lrcf': low-rank Cholesky factor of the “controllability” LQG Gramian,

  • 'lqg_o_lrcf': low-rank Cholesky factor of the “observability” LQG Gramian,

  • ('br_c_lrcf', gamma): low-rank Cholesky factor of the “controllability” bounded real Gramian,

  • ('br_o_lrcf', gamma): low-rank Cholesky factor of the “observability” bounded real Gramian.

  • 'pr_c_lrcf': low-rank Cholesky factor of the “controllability” positive real Gramian,

  • 'pr_o_lrcf': low-rank Cholesky factor of the “observability” positive real Gramian.

  Note

  For '*_lrcf' types, the method assumes the system is asymptotically stable. For '*_dense' types, the method assumes that the underlying Lyapunov equation has a unique solution, i.e. no pair of system poles adds to zero in the continuous-time case and no pair of system poles multiplies to one in the discrete-time case. Additionally, for 'pr_c_lrcf' and 'pr_o_lrcf', it is assumed that D + D^T is invertible.

• mu – Parameter values.

Returns:

• If typ ends with '_lrcf', then the Gramian factor as a VectorArray from self.A.source.

• If typ ends with '_dense', then the Gramian as a NumPy array.

h2_norm(mu=None)

Compute the \(\mathcal{H}_2\)-norm of the LTIModel.

Note

Assumes the system is asymptotically stable.

Parameters:

mu – Parameter values.

Returns:

norm – \(\mathcal{H}_2\)-norm.

hankel_norm(mu=None)

Compute the Hankel-norm of the LTIModel.

Note

Assumes the system is asymptotically stable.

Parameters:

mu – Parameter values.

Returns:

norm – Hankel-norm.

hinf_norm(mu=None, return_fpeak=False, ab13dd_equilibrate=False, tol=1e-10)

Compute the \(\mathcal{H}_\infty\)-norm of the LTIModel.

Note

Assumes the system is asymptotically stable. Under this is assumption the \(\mathcal{H}_\infty\)-norm is equal to the \(\mathcal{L}_\infty\)-norm. Accordingly, this method calls linf_norm.

Parameters:

• mu – Parameter values.

• return_fpeak – Whether to return the frequency at which the maximum is achieved.

• ab13dd_equilibrate – Whether slycot.ab13dd should use equilibration.

• tol – Tolerance in norm computation.

Returns:

• norm – \(\mathcal{H}_\infty\)-norm.

• fpeak – Frequency at which the maximum is achieved (if return_fpeak is True).

hsv(mu=None)

Hankel singular values.

Note

Assumes the system is asymptotically stable.

Parameters:

mu – Parameter values.

Returns:

sv – One-dimensional NumPy array of singular values.

impulse_resp(mu=None, return_solution=False)

Compute impulse response from all inputs.

Parameters:

• mu – Parameter values for which to solve.

• return_solution – If True, the model solution for the given parameter values mu is returned.

Returns:

• output – Impulse response as a 3D NumPy array where output.shape[0] is the number of outputs, and output.shape[1] is the number of time steps, output.shape[2] is the number of inputs.

• solution – The tuple of solution VectorArrays for every input. Returned only when return_solution is True.

l2_norm(mu=None)

Compute the \(\mathcal{L}_2\)-norm of the LTIModel.

The \(\mathcal{L}_2\)-norm of an LTIModel is defined via the integral

\[\lVert H \rVert_{\mathcal{L}_2} = \left( \frac{1}{2 \pi} \int_{-\infty}^{\infty} \lVert H(\boldsymbol{\imath} \omega) \rVert_{\operatorname{F}}^2 \operatorname{d}\!\omega \right)^{\frac{1}{2}}.\]

Parameters:

mu – Parameter.

Returns:

norm – \(\mathcal{L}_2\)-norm.

linf_norm(mu=None, return_fpeak=False, ab13dd_equilibrate=False, tol=1e-10)

Compute the \(\mathcal{L}_\infty\)-norm of the LTIModel.

The \(\mathcal{L}_\infty\)-norm of an LTIModel is defined via

\[\lVert H \rVert_{\mathcal{L}_\infty} = \sup_{\omega \in \mathbb{R}} \lVert H(\boldsymbol{\imath} \omega) \rVert_2.\]

Parameters:

• mu – Parameter.

• return_fpeak – Whether to return the frequency at which the maximum is achieved.

• ab13dd_equilibrate – Whether slycot.ab13dd should use equilibration.

• tol – Tolerance in norm computation.

Returns:

• norm – \(\mathcal{L}_\infty\)-norm.

• fpeak – Frequency at which the maximum is achieved (if return_fpeak is True).

moebius_substitution(M, sampling_time=0)

Create a transformed LTIModel by applying an arbitrary Moebius transformation.

This method returns a transformed LTIModel such that the transfer function of the original and transformed LTIModel are related by a Moebius substitution of the frequency argument:

\[H(s)=\tilde{H}(M(s)),\]

where

\[M(s) = \frac{as+b}{cs+d}\]

is a Moebius transformation. See [CCA96] for details.

Parameters:

• M – The MoebiusTransformation that defines the frequency mapping.

• sampling_time – The sampling time of the transformed system (in seconds). 0 if the system is continuous-time, otherwise a positive number. Defaults to zero.

Returns:

sys – The transformed LTIModel.

poles(mu=None)

Compute system poles.

Note

Assumes the systems is small enough to use a dense eigenvalue solver.

Parameters:

mu – Parameter values for which to compute the systems poles.

Returns:

One-dimensional |NumPy array| of system poles.

step_resp(mu=None, return_solution=False)

Compute step response from all inputs.

Parameters:

• mu – Parameter values for which to solve.

• return_solution – If True, the model solution for the given parameter values mu is returned.

Returns:

• output – Step response as a 3D NumPy array where output.shape[0] is the number of outputs, and output.shape[1] is the number of time steps, output.shape[2] is the number of inputs.

• solution – The tuple of solution VectorArrays for every input. Returned only when return_solution is True.

to_abcde_files(files_basename, mu=None)

Save operators as matrices to .[ABCDE] files in Matrix Market format.

Parameters:

• files_basename – The basename of files containing the operators.

• mu – The parameter values for which to write the operators to files.

to_abcde_matrices(format=None, mu=None)

Return A, B, C, D, and E operators as matrices.

Parameters:

• format – Format of the resulting matrices: NumPy array if ‘dense’, otherwise the appropriate SciPy spmatrix. If None, a choice between dense and sparse format is automatically made.

• mu – The parameter values for which to convert the operators.

Returns:

• A – The NumPy array or SciPy spmatrix A.

• B – The NumPy array or SciPy spmatrix B.

• C – The NumPy array or SciPy spmatrix C.

• D – The NumPy array or SciPy spmatrix D or None (if D is a ZeroOperator).

• E – The NumPy array or SciPy spmatrix E or None (if E is an IdentityOperator).

to_continuous(method='Tustin', w0=0)

Converts a discrete-time LTIModel to a continuous-time LTIModel.

Parameters:

• method – A string that defines the transformation method. At the moment only Tustin’s method is supported.

• w0 – If method=='Tustin', this parameter can be used to specify the prewarping-frequency. Defaults to zero.

Returns:

sys – Continuous-time LTIModel.

to_discrete(sampling_time, method='Tustin', w0=0)

Converts a continuous-time LTIModel to a discrete-time LTIModel.

Parameters:

• sampling_time – A positive number that denotes the sampling time of the resulting system (in seconds).

• method – A string that defines the transformation method. At the moment only Tustin’s method is supported.

• w0 – If method=='Tustin', this parameter can be used to specify the prewarping-frequency. Defaults to zero.

Returns:

sys – Discrete-time LTIModel.

to_files(A_file, B_file, C_file, D_file=None, E_file=None, mu=None)

Write operators to files as matrices.

Parameters:

• A_file – The name of the file (with extension) containing A.

• B_file – The name of the file (with extension) containing B.

• C_file – The name of the file (with extension) containing C.

• D_file – The name of the file (with extension) containing D or None if D is a ZeroOperator.

• E_file – The name of the file (with extension) containing E or None if E is an IdentityOperator.

• mu – The parameter values for which to write the operators to files.

to_mat_file(file_name, mu=None)

Save operators as matrices to .mat file.

Parameters:

• file_name – The name of the .mat file (extension .mat does not need to be included).

• mu – The parameter values for which to write the operators to files.

to_matrices(format=None, mu=None)

Return operators as matrices.

Parameters:

• format – Format of the resulting matrices: NumPy array if ‘dense’, otherwise the appropriate SciPy spmatrix. If None, a choice between dense and sparse format is automatically made.

• mu – The parameter values for which to convert the operators.

Returns:

• A – The NumPy array or SciPy spmatrix A.

• B – The NumPy array or SciPy spmatrix B.

• C – The NumPy array or SciPy spmatrix C.

• D – The NumPy array or SciPy spmatrix D or None (if D is a ZeroOperator).

• E – The NumPy array or SciPy spmatrix E or None (if E is an IdentityOperator).

class pymor.models.iosys.LinearDelayModel(A, Ad, tau, B, C, D=None, E=None, sampling_time=0, error_estimator=None, visualizer=None, name=None)

Bases: pymor.models.interface.Model

Class for linear delay systems.

This class describes input-state-output systems given by

\[\begin{split}E x'(t) & = A x(t) + \sum_{i = 1}^q{A_i x(t - \tau_i)} + B u(t), \\ y(t) & = C x(t) + D u(t),\end{split}\]

if continuous-time, or

\[\begin{split}E x(k + 1) & = A x(k) + \sum_{i = 1}^q{A_i x(k - \tau_i)} + B u(k), \\ y(k) & = C x(k) + D u(k),\end{split}\]

if discrete-time, where \(E\), \(A\), \(A_i\), \(B\), \(C\), and \(D\) are linear operators.

All methods related to the transfer function (e.g., frequency response calculation and Bode plots) are attached to the transfer_function attribute.

Parameters:

• A – The Operator A.

• Ad – The tuple of Operators A_i.

• tau – The tuple of delay times (positive floats or ints).

• B – The Operator B.

• C – The Operator C.

• D – The Operator D or None (then D is assumed to be zero).

• E – The Operator E or None (then E is assumed to be identity).

• sampling_time – 0 if the system is continuous-time, otherwise a positive number that denotes the sampling time (in seconds).

• error_estimator – An error estimator for the problem. This can be any object with an estimate_error(U, mu, model) method. If error_estimator is not None, an estimate_error(U, mu) method is added to the model which will call error_estimator.estimate_error(U, mu, self).

• visualizer – A visualizer for the problem. This can be any object with a visualize(U, model, ...) method. If visualizer is not None, a visualize(U, *args, **kwargs) method is added to the model which forwards its arguments to the visualizer’s visualize method.

• name – Name of the system.

order

The order of the system (equal to A.source.dim).

dim_input

The number of inputs.

dim_output

The number of outputs.

q

The number of delay terms.

tau

The tuple of delay times.

A

The Operator A.

Ad

The tuple of Operators A_i.

B

The Operator B.

C

The Operator C.

D

The Operator D.

E

The Operator E.

transfer_function

The transfer function.

cache_region = 'memory'

class pymor.models.iosys.LinearStochasticModel(A, As, B, C, D=None, E=None, sampling_time=0, error_estimator=None, visualizer=None, name=None)

Bases: pymor.models.interface.Model

Class for linear stochastic systems.

This class describes input-state-output systems given by

\[\begin{split}E \mathrm{d}x(t) & = A x(t) \mathrm{d}t + \sum_{i = 1}^q{A_i x(t) \mathrm{d}\omega_i(t)} + B u(t) \mathrm{d}t, \\ y(t) & = C x(t) + D u(t),\end{split}\]

if continuous-time, or

\[\begin{split}E x(k + 1) & = A x(k) + \sum_{i = 1}^q{A_i x(k) \omega_i(k)} + B u(k), \\ y(k) & = C x(k) + D u(t),\end{split}\]

if discrete-time, where \(E\), \(A\), \(A_i\), \(B\), \(C\), and \(D\) are linear operators and \(\omega_i\) are stochastic processes.

Parameters:

• A – The Operator A.

• As – The tuple of Operators A_i.

• B – The Operator B.

• C – The Operator C.

• D – The Operator D or None (then D is assumed to be zero).

• E – The Operator E or None (then E is assumed to be identity).

• sampling_time – 0 if the system is continuous-time, otherwise a positive number that denotes the sampling time (in seconds).

• error_estimator – An error estimator for the problem. This can be any object with an estimate_error(U, mu, model) method. If error_estimator is not None, an estimate_error(U, mu) method is added to the model which will call error_estimator.estimate_error(U, mu, self).

• visualizer – A visualizer for the problem. This can be any object with a visualize(U, model, ...) method. If visualizer is not None, a visualize(U, *args, **kwargs) method is added to the model which forwards its arguments to the visualizer’s visualize method.

• name – Name of the system.

order

The order of the system (equal to A.source.dim).

dim_input

The number of inputs.

dim_output

The number of outputs.

q

The number of stochastic processes.

A

The Operator A.

As

The tuple of Operators A_i.

B

The Operator B.

C

The Operator C.

D

The Operator D.

E

The Operator E.

cache_region = 'memory'

class pymor.models.iosys.PHLTIModel(J, R, G, P=None, S=None, N=None, E=None, Q=None, solver_options=None, error_estimator=None, visualizer=None, name=None)

Bases: LTIModel

Class for (continuous) port-Hamiltonian linear time-invariant systems.

This class describes input-state-output systems given by

\[\begin{split}E(\mu) \dot{x}(t, \mu) & = (J(\mu) - R(\mu)) Q(\mu) x(t, \mu) + (G(\mu) - P(\mu)) u(t), \\ y(t, \mu) & = (G(\mu) + P(\mu))^T Q(\mu) x(t, \mu) + (S(\mu) - N(\mu)) u(t),\end{split}\]

where \(H(\mu) = Q(\mu)^T E(\mu)\),

\[\begin{split}\Gamma(\mu) = \begin{bmatrix} J(\mu) & G(\mu) \\ -G(\mu)^T & N(\mu) \end{bmatrix}, \text{ and } \mathcal{W}(\mu) = \begin{bmatrix} R(\mu) & P(\mu) \\ P(\mu)^T & S(\mu) \end{bmatrix}\end{split}\]

satisfy \(H(\mu) = H(\mu)^T \succ 0\), \(\Gamma(\mu)^T = -\Gamma(\mu)\), and \(\mathcal{W}(\mu) = \mathcal{W}(\mu)^T \succcurlyeq 0\).

A dynamical system of this form, together with a given quadratic (energy) function \(\mathcal{H}(x, \mu) = \tfrac{1}{2} x^T H(\mu) x\), typically called Hamiltonian, is called a port-Hamiltonian system.

All methods related to the transfer function (e.g., frequency response calculation and Bode plots) are attached to the transfer_function attribute.

Parameters:

• J – The Operator J.

• R – The Operator R.

• G – The Operator G.

• P – The Operator P or None (then P is assumed to be zero).

• S – The Operator S or None (then S is assumed to be zero).

• N – The Operator N or None (then N is assumed to be zero).

• E – The Operator E or None (then E is assumed to be identity).

• Q – The Operator Q or None (then Q is assumed to be identity).

• solver_options – The solver options to use to solve the Lyapunov equations.

• name – Name of the system.

order

The order of the system.

dim_input

The number of inputs.

dim_output

The number of outputs.

J

The Operator J.

R

The Operator R.

G

The Operator G.

P

The Operator P.

S

The Operator S.

N

The Operator N.

E

The Operator E.

Q

The Operator Q.

transfer_function

The transfer function.

Methods

from_matrices	Create PHLTIModel from matrices.
from_passive_LTIModel	Convert a passive LTIModel to a PHLTIModel.
to_berlin_form	Convert the PHLTIModel into its Berlin form.
to_matrices	Return operators as matrices.

cache_region = 'memory'

classmethod from_matrices(J, R, G, P=None, S=None, N=None, E=None, Q=None, solver_options=None, error_estimator=None, visualizer=None, name=None)

Create PHLTIModel from matrices.

Parameters:

• J – The NumPy array or SciPy spmatrix J.

• R – The NumPy array or SciPy spmatrix R.

• G – The NumPy array or SciPy spmatrix G.

• P – The NumPy array or SciPy spmatrix P or None (then P is assumed to be zero).

• S – The NumPy array or SciPy spmatrix S or None (then S is assumed to be zero).

• N – The NumPy array or SciPy spmatrix N or None (then N is assumed to be zero).

• E – The NumPy array or SciPy spmatrix E or None (then E is assumed to be identity).

• Q – The NumPy array or SciPy spmatrix Q or None (then Q is assumed to be identity).

• solver_options – The solver options to use to solve the Lyapunov equations.

• error_estimator – An error estimator for the problem. This can be any object with an estimate_error(U, mu, model) method. If error_estimator is not None, an estimate_error(U, mu) method is added to the model which will call error_estimator.estimate_error(U, mu, self).

• visualizer – A visualizer for the problem. This can be any object with a visualize(U, model, ...) method. If visualizer is not None, a visualize(U, *args, **kwargs) method is added to the model which forwards its arguments to the visualizer’s visualize method.

• name – Name of the system.

Returns:

phlti – The PHLTIModel with operators J, R, G, P, S, N, and E.

classmethod from_passive_LTIModel(model)

Convert a passive LTIModel to a PHLTIModel.

Note

The method uses dense computations and converts model to dense matrices.

Parameters:

• model – The passive LTIModel to convert.

• generalized – If True, the resulting PHLTIModel will have \(Q=I\).

to_berlin_form()

Convert the PHLTIModel into its Berlin form.

Returns a PHLTIModel with \(Q=I\), by left multiplication with \(Q^T\).

Returns:

model – PHLTIModel with \(Q=I\).

to_matrices(format=None, mu=None)

Return operators as matrices.

Parameters:

• format – Format of the resulting matrices: NumPy array if ‘dense’, otherwise the appropriate SciPy spmatrix. If None, a choice between dense and sparse format is automatically made.

• mu – The parameter values for which to convert the operators.

Returns:

• J – The NumPy array or SciPy spmatrix J.

• R – The NumPy array or SciPy spmatrix R.

• G – The NumPy array or SciPy spmatrix G.

• P – The NumPy array or SciPy spmatrix P or None (if P is a ZeroOperator).

• S – The NumPy array or SciPy spmatrix S or None (if S is a ZeroOperator).

• N – The NumPy array or SciPy spmatrix N or None (if N is a ZeroOperator).

• E – The NumPy array or SciPy spmatrix E or None (if E is an IdentityOperator).

• Q – The NumPy array or SciPy spmatrix Q or None (if Q is an IdentityOperator).

class pymor.models.iosys.SecondOrderModel(M, E, K, B, Cp, Cv=None, D=None, sampling_time=0, solver_options=None, error_estimator=None, visualizer=None, name=None)

Bases: pymor.models.interface.Model

Class for linear second order systems.

This class describes input-output systems given by

\[\begin{split}M(\mu) \ddot{x}(t, \mu) + E(\mu) \dot{x}(t, \mu) + K(\mu) x(t, \mu) & = B(\mu) u(t), \\ y(t, \mu) & = C_p(\mu) x(t, \mu) + C_v(\mu) \dot{x}(t, \mu) + D(\mu) u(t),\end{split}\]

if continuous-time, or

\[\begin{split}M(\mu) x(k + 2, \mu) + E(\mu) x(k + 1, \mu) + K(\mu) x(k, \mu) & = B(\mu) u(k), \\ y(k, \mu) & = C_p(\mu) x(k, \mu) + C_v(\mu) x(k + 1, \mu) + D(\mu) u(k),\end{split}\]

if discrete-time, where \(M\), \(E\), \(K\), \(B\), \(C_p\), \(C_v\), and \(D\) are linear operators.

All methods related to the transfer function (e.g., frequency response calculation and Bode plots) are attached to the transfer_function attribute.

Parameters:

• M – The Operator M.

• E – The Operator E.

• K – The Operator K.

• B – The Operator B.

• Cp – The Operator Cp.

• Cv – The Operator Cv or None (then Cv is assumed to be zero).

• D – The Operator D or None (then D is assumed to be zero).

• sampling_time – 0 if the system is continuous-time, otherwise a positive number that denotes the sampling time (in seconds).

• solver_options – The solver options to use to solve the Lyapunov equations.

• error_estimator – An error estimator for the problem. This can be any object with an estimate_error(U, mu, model) method. If error_estimator is not None, an estimate_error(U, mu) method is added to the model which will call error_estimator.estimate_error(U, mu, self).

• visualizer – A visualizer for the problem. This can be any object with a visualize(U, model, ...) method. If visualizer is not None, a visualize(U, *args, **kwargs) method is added to the model which forwards its arguments to the visualizer’s visualize method.

• name – Name of the system.

order

The order of the system (equal to M.source.dim).

dim_input

The number of inputs.

dim_output

The number of outputs.

M

The Operator M.

E

The Operator E.

K

The Operator K.

B

The Operator B.

Cp

The Operator Cp.

Cv

The Operator Cv.

D

The Operator D.

transfer_function

The transfer function.

Methods

from_files	Create LTIModel from matrices stored in separate files.
from_matrices	Create a second order system from matrices.
gramian	Compute a second-order Gramian.
h2_norm	Compute the \(\mathcal{H}_2\)-norm.
hankel_norm	Compute the Hankel-norm.
hinf_norm	Compute the \(\mathcal{H}_\infty\)-norm.
poles	Compute system poles.
psv	Position singular values.
pvsv	Position-velocity singular values.
to_files	Write operators to files as matrices.
to_lti	Return a first order representation.
to_matrices	Return operators as matrices.
vpsv	Velocity-position singular values.
vsv	Velocity singular values.

cache_region = 'memory'

classmethod from_files(M_file, E_file, K_file, B_file, Cp_file, Cv_file=None, D_file=None, sampling_time=0, solver_options=None, error_estimator=None, visualizer=None, name=None)

Create LTIModel from matrices stored in separate files.

Parameters:

• M_file – The name of the file (with extension) containing A.

• E_file – The name of the file (with extension) containing E.

• K_file – The name of the file (with extension) containing K.

• B_file – The name of the file (with extension) containing B.

• Cp_file – The name of the file (with extension) containing Cp.

• Cv_file – None or the name of the file (with extension) containing Cv.

• D_file – None or the name of the file (with extension) containing D.

• sampling_time – 0 if the system is continuous-time, otherwise a positive number that denotes the sampling time (in seconds).

• solver_options – The solver options to use to solve the Lyapunov equations.

• error_estimator – An error estimator for the problem. This can be any object with an estimate_error(U, mu, model) method. If error_estimator is not None, an estimate_error(U, mu) method is added to the model which will call error_estimator.estimate_error(U, mu, self).

• visualizer – A visualizer for the problem. This can be any object with a visualize(U, model, ...) method. If visualizer is not None, a visualize(U, *args, **kwargs) method is added to the model which forwards its arguments to the visualizer’s visualize method.

• name – Name of the system.

Returns:

some – The SecondOrderModel with operators M, E, K, B, Cp, Cv, and D.

classmethod from_matrices(M, E, K, B, Cp, Cv=None, D=None, sampling_time=0, solver_options=None, error_estimator=None, visualizer=None, name=None)

Create a second order system from matrices.

Parameters:

• M – The NumPy array or SciPy spmatrix M.

• E – The NumPy array or SciPy spmatrix E.

• K – The NumPy array or SciPy spmatrix K.

• B – The NumPy array or SciPy spmatrix B.

• Cp – The NumPy array or SciPy spmatrix Cp.

• Cv – The NumPy array or SciPy spmatrix Cv or None (then Cv is assumed to be zero).

• D – The NumPy array or SciPy spmatrix D or None (then D is assumed to be zero).

• sampling_time – 0 if the system is continuous-time, otherwise a positive number that denotes the sampling time (in seconds).

• solver_options – The solver options to use to solve the Lyapunov equations.

• error_estimator – An error estimator for the problem. This can be any object with an estimate_error(U, mu, model) method. If error_estimator is not None, an estimate_error(U, mu) method is added to the model which will call error_estimator.estimate_error(U, mu, self).

• visualizer – A visualizer for the problem. This can be any object with a visualize(U, model, ...) method. If visualizer is not None, a visualize(U, *args, **kwargs) method is added to the model which forwards its arguments to the visualizer’s visualize method.

• name – Name of the system.

Returns:

lti – The SecondOrderModel with operators M, E, K, B, Cp, Cv, and D.

gramian(typ, mu=None)

Compute a second-order Gramian.

Parameters:

• typ –

  The type of the Gramian:

  • 'pc_lrcf': low-rank Cholesky factor of the position controllability Gramian,

  • 'vc_lrcf': low-rank Cholesky factor of the velocity controllability Gramian,

  • 'po_lrcf': low-rank Cholesky factor of the position observability Gramian,

  • 'vo_lrcf': low-rank Cholesky factor of the velocity observability Gramian,

  • 'pc_dense': dense position controllability Gramian,

  • 'vc_dense': dense velocity controllability Gramian,

  • 'po_dense': dense position observability Gramian,

  • 'vo_dense': dense velocity observability Gramian.

  Note

  For '*_lrcf' types, the method assumes the system is asymptotically stable. For '*_dense' types, the method assumes that the underlying Lyapunov equation has a unique solution, i.e. no pair of system poles adds to zero in the continuous-time case and no pair of system poles multiplies to one in the discrete-time case.

• mu – Parameter values.

Returns:

• If typ is 'pc_lrcf', 'vc_lrcf', 'po_lrcf' or 'vo_lrcf', then the Gramian factor as a

• VectorArray from self.M.source.

• If typ is 'pc_dense', 'vc_dense', 'po_dense' or 'vo_dense', then the Gramian as a

• |NumPy array|.

h2_norm(mu=None)

Compute the \(\mathcal{H}_2\)-norm.

Note

Assumes the system is asymptotically stable.

Parameters:

mu – Parameter values.

Returns:

norm – \(\mathcal{H}_2\)-norm.

hankel_norm(mu=None)

Compute the Hankel-norm.

Note

Assumes the system is asymptotically stable.

Parameters:

mu – Parameter values.

Returns:

norm – Hankel-norm.

hinf_norm(mu=None, return_fpeak=False, ab13dd_equilibrate=False, tol=1e-10)

Compute the \(\mathcal{H}_\infty\)-norm.

Note

Assumes the system is asymptotically stable.

Parameters:

• mu – Parameter values.

• return_fpeak – Should the frequency at which the maximum is achieved should be returned.

• ab13dd_equilibrate – Should slycot.ab13dd use equilibration.

• tol – Tolerance in norm computation.

Returns:

• norm – \(\mathcal{H}_\infty\).

• fpeak – Frequency at which the maximum is achieved (if return_fpeak is True).

poles(mu=None)

Compute system poles.

Note

Assumes the systems is small enough to use a dense eigenvalue solver.

Parameters:

mu – Parameter values.

Returns:

One-dimensional |NumPy array| of system poles.

psv(mu=None)

Position singular values.

Note

Assumes the system is asymptotically stable.

Parameters:

mu – Parameter values.

Returns:

One-dimensional |NumPy array| of singular values.

pvsv(mu=None)

Position-velocity singular values.

Note

Assumes the system is asymptotically stable.

Parameters:

mu – Parameter values.

Returns:

One-dimensional |NumPy array| of singular values.

to_files(M_file, E_file, K_file, B_file, Cp_file, Cv_file=None, D_file=None, mu=None)

Write operators to files as matrices.

Parameters:

• M_file – The name of the file (with extension) containing M.

• E_file – The name of the file (with extension) containing E.

• K_file – The name of the file (with extension) containing K.

• B_file – The name of the file (with extension) containing B.

• Cp_file – The name of the file (with extension) containing Cp.

• Cv_file – The name of the file (with extension) containing Cv or None if D is a ZeroOperator.

• D_file – The name of the file (with extension) containing D or None if D is a ZeroOperator.

• mu – The parameter values for which to write the operators to files.

to_lti()

Return a first order representation.

The first order representation

\[\begin{split}\begin{bmatrix} I & 0 \\ 0 & M \end{bmatrix} \frac{\mathrm{d}}{\mathrm{d}t}\! \begin{bmatrix} x(t) \\ \dot{x}(t) \end{bmatrix} & = \begin{bmatrix} 0 & I \\ -K & -E \end{bmatrix} \begin{bmatrix} x(t) \\ \dot{x}(t) \end{bmatrix} + \begin{bmatrix} 0 \\ B \end{bmatrix} u(t), \\ y(t) & = \begin{bmatrix} C_p & C_v \end{bmatrix} \begin{bmatrix} x(t) \\ \dot{x}(t) \end{bmatrix} + D u(t)\end{split}\]

is returned.

Returns:

lti – LTIModel equivalent to the second-order model.

to_matrices(mu=None)

Return operators as matrices.

Parameters:

mu – The parameter values for which to convert the operators.

Returns:

• M – The NumPy array or SciPy spmatrix M.

• E – The NumPy array or SciPy spmatrix E.

• K – The NumPy array or SciPy spmatrix K.

• B – The NumPy array or SciPy spmatrix B.

• Cp – The NumPy array or SciPy spmatrix Cp.

• Cv – The NumPy array or SciPy spmatrix Cv or None (if Cv is a ZeroOperator).

• D – The NumPy array or SciPy spmatrix D or None (if D is a ZeroOperator).

vpsv(mu=None)

Velocity-position singular values.

Note

Assumes the system is asymptotically stable.

Parameters:

mu – Parameter values.

Returns:

One-dimensional |NumPy array| of singular values.

vsv(mu=None)

Velocity singular values.

Note

Assumes the system is asymptotically stable.

Parameters:

mu – Parameter values.

Returns:

One-dimensional |NumPy array| of singular values.

class pymor.models.iosys.StepFunction(dim_input, component, sampling_time)

Bases: pymor.analyticalproblems.functions.Function

Interface for Parameter dependent analytical functions.

Every Function is a map of the form

f(μ): Ω ⊆ R^d -> R^(shape_range)

The returned values are NumPy arrays of arbitrary (but fixed) shape. Note that NumPy distinguishes between one-dimensional arrays of length 1 (with shape (1,)) and zero-dimensional scalar arrays (with shape ()). In pyMOR, we usually expect scalar-valued functions to have shape_range == ().

While the function might raise an error if it is evaluated for an argument not in the domain Ω, the exact behavior is left undefined.

Functions are vectorized in the sense, that if x.ndim == k, then

f(x, μ)[i0, i1, ..., i(k-2)] == f(x[i0, i1, ..., i(k-2)], μ).

In particular, f(x, μ).shape == x.shape[:-1] + shape_range.

dim_domain

The dimension d > 0.

shape_range

The shape of the function values.

Methods

evaluate

Evaluate the function for given argument x and parameter values mu.

evaluate(x, mu=None)

Evaluate the function for given argument x and parameter values mu.

pymor.models.iosys.sparse_min_size(value=1000)

Return minimal sparse problem size for which to warn about converting to dense.