pyvib.pnlss module

class pyvib.pnlss.PNLSS(*system, **kwargs)

Bases: pyvib.statespace.NonlinearStateSpace, pyvib.statespace.StateSpaceIdent

jacobian(x0, weight=False)

nlterms(eq, degree, structure)

Set active nonlinear terms/monomials to be optimized

output(u, t=None, x0=None)

pyvib.pnlss.combinations(n, degrees)

Lists all nonlinear terms in a multivariate polynomial.

Lists the exponents of all possible monomials in a multivariate polynomial with n inputs. Only the nonlinear degrees in degrees are considered.

Parameters:

• n (int) – number of inputs
• degrees (ndarray) – array with the degrees of nonlinearity

Returns:

monomials – matrix of exponents

Return type:

ndarray(ncomb,n)

Examples

A polynomial with all possible quadratic and cubic terms in the variables x and y contains the monomials x*x, x*y, y*y, x*x*x, x*x*y, x*y*y, and y*y*y. >>> combinations(2,[2, 3]) array([[2, 0], # -> x^2 * y^0 = x*x

[1, 1], # -> x^1 * y^1 = x*y [0, 2], # -> x^0 * y^2 = y*y [3, 0], # -> x^3 * y^0 = x*x*x [2, 1], # -> x^2 * y^1 = x*x*y [1, 2], # -> x^1 * y^2 = x*y*y [0, 3]]) # -> x^0 * y^3 = y*y*y

Element (i,j) of out indicates the power to which variable j is raised in monomial i. For example, out[4] = [2,1], which means that the fifth monomial is equal to x^2*y^1 = x*x*y.

pyvib.pnlss.dnlsim(system, u, t=None, x0=None)

Simulate output of a discrete-time nonlinear system.

Calculate the output and the states of a nonlinear state-space model.

x(t+1) = A x(t) + B u(t) + E zeta(x(t),u(t)) y(t) = C x(t) + D u(t) + F eta(x(t),u(t))

where zeta and eta are polynomials whose exponents are given in xpowers and ypowers, respectively. The maximum degree in one variable (a state or an input) in zeta or eta is given in max_nx and max_ny, respectively. The initial state is given in x0.

pyvib.pnlss.element_jacobian(samples, A_Edwdx, C_Fdwdx, active)

Compute Jacobian of the output y wrt. A, B, and E

The Jacobian is calculated by filtering an alternative state-space model

∂x∂Aᵢⱼ(t+1) = Iᵢⱼx(t) + (A + E*∂ζ∂x)*∂x∂Aᵢⱼ(t) ∂y∂Aᵢⱼ(t) = (C + F*∂η∂x)*∂x∂Aᵢⱼ(t)

where JA = ∂y∂Aᵢⱼ

Parameters:

• samples (ndarray) – x, u or zeta corresponding to JA, JB, or JE
• A_Edwdx (ndarray (n,n,NT)) – The result of A + E*∂ζ∂x
• C_Fdwdx (ndarray (p,n,NT)) – The result of C + F*∂η∂x
• active (ndarray) – Array with index of active elements. For JA: np.arange(n**2), JB: n*m or JE: xactive

Returns:

• JA, JB or JE depending on the samples given as input
• See fJNL

pyvib.pnlss.hom_combinations(n, degree)

Lists the exponents of all possible terms in a homogeneous polynomial monomial representation, e.g. [1 2] represents x1*x2**2

Examples

>>> hom_combinations(2,2)
array([[2, 0],
       [1, 1],
       [0, 2]])

pyvib.pnlss.jacobian(x0, system, weight=False)

Compute the Jacobians of a steady state nonlinear state-space model

Jacobians of a nonlinear state-space model

x(t+1) = A x(t) + B u(t) + E zeta(x(t),u(t)) y(t) = C x(t) + D u(t) + F eta(x(t),u(t))

i.e. the partial derivatives of the modeled output w.r.t. the active elements in the A, B, E, F, D, and C matrices, fx: JA = ∂y/∂Aᵢⱼ

x0 : ndarray

flattened array of state space matrices

pyvib.pnlss.select_active(structure, n, m, q, nx)

Select active elements in E or F matrix.

Select the active elements (i.e. those on which optimization will be done) in the E or F matrix. In particular, the linear indices (see also sub2ind and ind2sub) of the active elements in the transpose of the E or F matrix are calculated.

structure: str

string indicating which elements in the E or F matrix are active. ‘diagonal’: active elements in row j of the E matrix are those

corresponding to pure nonlinear terms in state j (only for state equation)

‘inputsonly’ : only terms in inputs ‘statesonly’ : only terms in states ‘nocrossprod’ : no cross-terms ‘affine’ : only terms that are linear in one state ‘affinefull’ : only terms that are linear in one state or constant in

the states

‘full’ : all terms ‘empty’ : no terms ‘nolastinput’ : no terms in last input row_E : only row given by row_E in E matrix is active (only for state

equation)

n : int

number of states

m : int

number of inputs

q : int

number of rows in corresponding E/F matrix

q = n if E matrix is considered, q = p if F matrix is considered

nx : int | list

degrees of nonlinearity in E/F matrix

Returns:active – F matrix Return type:linear indices of the active elements in the transpose of the E or

Examples

>>> n = 2  # Number of states
>>> m = 1  # Number of inputs
>>> p = 1  # Number of outputs
>>> nx = 2 # Degree(s) of nonlinearity
Powers of all possible terms in n+m inputs of degree(s) nx
>>> terms = combinations(n+m,nx)
array([[2, 0, 0],
       [1, 1, 0],
       [1, 0, 1],
       [0, 2, 0],
       [0, 1, 1],
       [0, 0, 2]])

There are six quadratic terms in the two states x1 and x2, and the input u, namely x1^2, x1*x2, x1*u, x2^2, x2*u, and u^2. The matrix E is a 2 x 6 matrix that contains the polynomial coefficients in each of these 6 terms for both state updates. The active elements will be calculated as linear indices in the transpose of E, hence E can be represented as E = [e1 e2 e3 e4 e5 e6;

e7 e8 e9 e10 e11 e12]

The matrix F is a 1 x 6 matrix that contains the polynomial coefficients in each of the 6 terms for the output equation. The matrix F can be represented as F = [f1 f2 f3 f4 f5 f6]

Diagonal structure >>> activeE = select_active(‘diagonal’,n,m,n,nx) array([0, 9]) Only e1 and e10 are active. This corresponds to a term x₁² in the first state equation and a term x₂² in the second state equation.

Inputs only structure >>> activeE = select_active(‘inputsonly’,n,m,n,nx) array([ 5, 11]) Only e6 and e12 are active. This corresponds to a term u² in both state equations. In all other terms, at least one of the states (possibly raised to a certain power) is a factor.

>>> activeF = select_active('inputsonly',n,m,p,nx)
array([5])
Only f6 is active. This corresponds to a term u² in the output equation.

States only structure >>> activeE = select_active(‘statesonly’,n,m,n,nx) array([0, 1, 3, 6, 7, 9])

Only e1, e2, e4, e7, e8, and e10 are active. This corresponds to terms x₁², x₁*x₂, and x₂² in both state equations. In all other terms, the input (possibly raised to a certain power) is a factor.

No cross products structure >>> activeE = select_active(‘nocrossprod’,n,m,n,nx) array([ 0, 3, 5, 6, 9, 11]) Only e1, e4, e6, e7, e10, and e12 are active. This corresponds to terms x₁², x₂², and u² in both state equations. All other terms are crossterms where more than one variable is present as a factor.

State affine structure >>> activeE = select_active(‘affine’,n,m,n,nx) array([ 2, 4, 8, 10]) Only e3, e5, e9, and e11 are active. This corresponds to terms x₁*u and x₂*u in both state equations, since in these terms only one state appears, and it appears linearly.

Full state affine structure >>> activeE = select_active(‘affinefull’,n,m,n,nx) array([ 2, 4, 5, 8, 10, 11]) Only e3, e5, e6, e9, e11, and e12 are active. This corresponds to terms x₁*u, x₂*u and u² in both state equations, since in these terms at most one state appears, and if it appears, it appears linearly.

Full structure >>> activeE = select_active(‘full’,n,m,n,nx) array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]) All elements in the E matrix are active.

Empty structure >>> activeE = select_active(‘empty’,n,m,n,nx) array([], dtype=int64) None of the elements in the E matrix are active.

** One row in E matrix structure** row_E = 1 # Select which row in E is active >>> activeE = select_active(row_E,n,m,n,nx) array([ 6, 7, 8, 9, 10, 11]) Only the elements in the second row of E are active

No terms in last input structure This is useful in a PNLSS model when considering the initial state as a parameter. The state at time one can be estimated by adding an extra input u_art(t) that is equal to one at time zero and zero elsewhere. Like this, an extended PNLSS model is estimated, where the last column in its B matrix corresponds to the state at time one in the original PNLSS model. To ensure that the optimization is only carried out on the parameters of the original PNLSS model, only the corresponding coefficients in the E/F matrix should be selected as active.

Powers of all possible terms with one extra input >>> terms_extended = combinations(n+m+1,nx) array([[2, 0, 0, 0],

[1, 1, 0, 0], [1, 0, 1, 0], [1, 0, 0, 1], [0, 2, 0, 0], [0, 1, 1, 0], [0, 1, 0, 1], [0, 0, 2, 0], [0, 0, 1, 1], [0, 0, 0, 2]])

The nonlinear terms in the extra input should not be considered for optimization. >>> activeE_extended = select_active(‘nolastinput’,n,m+1,n,nx) array([ 0, 1, 2, 4, 5, 7, 10, 11, 12, 14, 15, 17])

Only the terms where the last input is raised to a power zero are active. This corresponds to the elif structure == where all terms in the original PNLSS model are active.

The example below illustrates how to combine a certain structure in the original model (e.g. ‘nocrossprod’) with the estimation of the initial state.

>>> activeE_ext = select_active('nolastinput',n,m+1,n,nx)
>>> activeE_ext = activeE_ext[select_active('nocrossprod',n,m,n,nx)]
array([ 0,  4,  7, 10, 14, 17])

This corresponds to the terms x₁², x₂², and u₁² in both rows of the E_extended matrix, and thus to all terms in the original model, except for the cross terms.

Note that an alternative approach is to include the initial state in the parameter vector (TODO see also fLMnlssWeighted_x0u0).