formulating one overall representation of the system via the so-called embedding principle, i.e., the global form.

In general, global models capture the dynamics of the system w.r.t. the scheduling signal and, for local models, the scheduling signal is used to interpolate (schedule) between the local linear time-invariant models.
The LPV analysis, modelling, identification, and control literature has become mature in recent years. To accelerate the applicability of LPV methodologies to applications, various computational tools have become available, e.g., LPVTools [1]}, Control Systems and System Identification Toolbox in MATLAB [2]}, Predictor-Based Subspace IDentification (PBSID) Toolbox [3]}, and LPV Input/Output Systems Identification Toolbox (LPVIOID) [4]}. The analysis and control methods in the LPVTools and Control Systems toolboxes are mainly focused on a local grid based formulation (Jacobian-linearizations) or a grid based linear fractional representation. The Control Systems Toolbox does include a limited set of tools based on a global form. The PBSID and LPVIOID toolboxes focus on specific algorithms for identification of global LPV models in either SS or IO form, respectively, but with limited options to further utilize these models and with only support for affine dependency on the scheduling signal.
The concept of the initial release of the LPVcore toolbox is to facilitate modelling, simulation, and identification of LPV representations in a global form in MATLAB. Modelling and simulating systems can be performed based on a continuous-time (CT) or discrete-time (DT) LPV representation. To include the parameter-varying nature, the real and complex number based matrix representations in MATLAB are extended to facilitate parameter-variations. These parameter-varying matrix functions form the core of the toolbox to facilitate algebraic operations. In addition, the LPVcore toolbox includes identification methods to identify discrete-time LPV IO, SS and LFR models. The initial release of the toolbox does not cover the comprehensive LPV literature. However, it provides a solid basis to incorporate LPV analysis and control methods, implement realization tools in a computationally efficient manner, and extend the set of identification methods in the future.
The intention of this paper is to provide an overview of the methods in the initial release of the LPVcore toolbox. In Section , the handling of parameter-varying matrix functions and the various LPV representation forms are introduced. The LPV-IO, LPV-SS and LPV-LFR identification methods are described in Section . Section  presents an example in which LPVcore is used to estimate a model of a DC motor with unbalanced disc. In Section , the conclusions are given.

LPV representations

In this section, the LPV representations are described to model the system at hand. The key element in formulating the representations in a global form is the parameter-varying matrix object of the LPVcore toolbox, as discussed in Section REF . With the parameter-varying matrix object, LPV input-output, state-space and linear-fractional representations are constructed, as detailed in Sections REF ,  REF  and REF .

Parameter-varying matrix function

The LPV representations are defined using parameter-varying matrix functions. MATLAB does not have a native implementation of these functions. Therefore, LPVcore introduces the |pmatrix| object, defining the following parameter-varying matrix function:
\((A \diamond p) = A_0 + \sum _{i=1}^{n_{\alpha }} A_i \left(\alpha _{i}\diamond p \right), \) 
where \(p:\mathbb {T}\rightarrow \mathbb {P}\subseteq \mathbb {R}^{{n_\mathrm {p}}}\)  is the scheduling variable, \(A_i\in \mathbb {R}^{k\times l}\)  are the parameters, and \(\alpha _{i}\in {R}\)  are scheduling-dependent coefficient functions. For the CT case, \(\mathbb {T}=\mathbb {R}\)  and, for the DT case, \(\mathbb {T}=\mathbb {Z}\) . \(\mathbb {P}\)  is a compact subset of \(\mathbb {R}^{{n_\mathrm {p}}}\) . The set \({R}\)  is defined as the set of real-analytic functions of the form \(f:\mathbb {R}^{{n_\mathrm {p}}}\times \ldots \times \mathbb {R}^{{n_\mathrm {p}}}\rightarrow \mathbb {R}\) . Furthermore, in the CT case, for a \(p\in \mathbb {P}\) , we define the following notation: if \(f\in {R}\) , then \((f\diamond p)=f(p,\tfrac{\mathrm {d}}{\mathrm {d}t}p,\tfrac{\mathrm {d}^{2}}{\mathrm {d}t^{2}} p,\ldots )\)  where \(\tfrac{\mathrm {d}}{\mathrm {d}t}p\)  is the time derivative of the scheduling signal. Similarly, in the DT case, for a \(p\in \mathbb {P}\) , we define: if \(f\in {R}\) , then \((f\diamond p)=f(p,q^1 p,q^{-1} p,q^{2}p,\ldots )\)  where \(q\)  is the time-shift operator, i.e., \(q p_t=p_{t+1}\) . Hence, the set \({R}\)  allows to include various types of non-linearities and time-varying effects into a representation. For example, the affine basis functions \(\alpha _{i}\in {R}\)  can represent polynomials, cosines, exponentials, logarithms, and rational functions in \(\tfrac{\mathrm {d}^{i}}{\mathrm {d}t^{i}} p\)  or \(p_{t\pm i}\)  with \(i\ge 0\) . For the full mathematical treatment of \({R}\)  see [5]}. Moreover, for notational simplicity, we denote by \({R}^{k\times l}\)  the set of all functions of the form (REF ).
In LPVcore, the |timemap| object is introduced to simplify the handling of dynamic dependency on \(p\) , by introducing the extended scheduling signal \(\rho \)  that includes all the shifted versions (DT) or time derivatives (CT) of \(p\)  that arise from the use of the \(\diamond \)  operator in (REF ). For example, to represent dynamic dependency on both \(p_t\)  and \(p_{t-1}\) , the user can generate the extended scheduling signal using |timemap([0, -1], 'dt')|. The first argument represents the time shifts or derivatives and the second argument denotes the time domain (|'dt'| and |'ct'| for DT and CT, respectively).
Then, to model a parameter-varying matrix function in LPVcore, two more ingredients are required:
the coefficients \(\lbrace A_i\rbrace _{i=0}^{n_\alpha }\) , and

the basis functions \(\lbrace \alpha _i\rbrace _{i=1}^{n_\alpha }\) .

To support simple use of parameter-varying matrix functions for casual users, but also support the more advanced needs of researchers, a flexible |pmatrix| object is defined, which represents a parameter-varying matrix function based on the three ingredients described above. The first argument of |pmatrix| specifies the coefficients, followed by name-value pairs that describe the basis functions and the extended scheduling signal. To demonstrate the degree of flexibility in specifying parameter-varying matrix functions, the use of |pmatrix| is shown in a few examples:

Consider the following affine parameter-varying matrix function:
\(\left( A \diamond p \right)_t = A_0 + A_1 p_t + A_2 p_{t-1}.\) 
The detailed way to represent this function in LPVcore is to first construct a |timemap| object to extend the scheduling signal with the time-shifted version, and then use this object to create a |pmatrix|, i.e.:

rho = timemap([0, -1], 'dt')
A = pmatrix(cat(3, A0, A1, A2), ...
'BasisType', 'affine', ...
'BasisParametrization', 0, 1, 2, ...
'SchedulingTimeMap', rho)

Note that the parameter |'BasisParametrization'| is a cell array of indices of the extended scheduling signal, with 0 representing the constant factor 1. Indices 1 and 2 refer to \(p_t\)  and \(p_{t-1}\) , respectively. As an alternative to |pmatrix|, one can simply use |preal| and assemble (REF ) in one line by algebraic operations (see Example REF ).


Consider the following affine parameter-varying matrix function with two scheduling signals:
\(\left( A \diamond (p, q) \right)_t = p_{t-1} + 2 q_{t-1} + 3 p_t + 4 q_t\) 
Additional scheduling signals can be specified in the |'Name'| option of |timemap|. Furthermore, affine basis type is the default and need not be specified when creating the |pmatrix|:

rho = timemap([-1, 0], 'dt', ...
'Name', 'p', 'q')
A = pmatrix(cat(3, 1, 2, 3, 4), ...
'BasisParametrization', ...
1, 2, 3, 4, ...
'SchedulingTimeMap', rho)



Consider the following polynomial parameter-varying matrix function:
\(\left( A \diamond p \right)_t = p_t^{2} + 2 p_{t} p_{t-2} + p_{t-2}^2.\) 
Setting |'BasisType'| to |'poly'| specifies polynomial dependence, for which the associated parametrization is a vector with the degree of each term of the extended scheduling signal, i.e.:

rho = timemap([-2, 0], 'dt')
A = pmatrix(cat(3, 1, 2, 1), ...
'BasisType', 'poly', ...
'BasisParametrization', ...
[2, 0], [1, 1], [0, 2], ...
'SchedulingTimeMap', rho)

Custom basis functions are supported, too.
In LPVcore, many key operations available in MATLAB for constant matrices, such as taking the sum or product, are extended for |pmatrix| objects, providing an intuitive manner to manipulate these objects. The |pmatrix| object automatically merges different scheduling signals when these operations are used (see Table REF  for an overview).
<TABLE>For convenience, the function |preal| is also available to generate a scalar coefficient. Combined with |pdiff| and |pshift| to differentiate or shift the scheduling in CT or DT respectively, it can be used as an alternative syntax for constructing parameter-varying matrix functions, e.g.:

The parameter-varying matrix function (REF ) can be constructed by a combination of |preal| and algebraic operations in the following way:

p = preal('p', 'dt')
A = A0 + A1 * p + A2 * pshift(p, -1)

As highlighted by the above examples, an important strength of LPVcore is the implementation of the parameter-varying matrix functions and, as shown in the next sections, the ability to model, analyze, identify, and control LPV representations in a global form. Local representational forms can be obtained via conversion functions, which will be introduced for compatibility with LPVTools in the future.

LPV-IO representation

With the parameter-varying matrix function defined, various forms of LPV representations of a system at hand can be formulated. A commonly used representation form is the input-output (IO) representation, given by the following equation:

\((\mathcal {A}(\xi )\diamond p)_t\ y_t=(\mathcal {B}(\xi )\diamond p)_t\ u_t,  \vspace{1.42262pt}\) 
where \(y:\mathbb {T}\rightarrow \mathbb {Y}=\mathbb {R}^{{n_\mathrm {y}}}\)  is the measured output signal and \(u:\mathbb {T}\rightarrow \mathbb {U}=\mathbb {R}^{{n_\mathrm {u}}}\)  denotes the input signal. For the CT case, \(\xi =\tfrac{\mathrm {d}}{\mathrm {d}t}\) , i.e., \(\xi y_t=\tfrac{\mathrm {d}}{\mathrm {d}t}y_t\)  denotes the time derivative of the output; and for the DT case, \(\xi = q^{-1}\)  is the backward time-shift operator. The scheduling-dependent polynomials \({A}(\xi )\) , \({B}(\xi )\)  are polynomials in the indeterminate \(\xi \)  with scheduling-dependent coefficients:

\(\mathcal {A}(\xi )\diamond p&=I_{{n_\mathrm {y}}}+\sum _{i=1}^{n_\mathrm {a}}(A_i\diamond p)\, \xi ^{i}, \\[-2mm]\mathcal {B}(\xi )\diamond p&=\sum _{j=0}^{n_\mathrm {b}}(B_j\diamond p)\, \xi ^{j},\) 


with finite polynomial orders \({n_\mathrm {a}}\ge 0\)  and \({n_\mathrm {b}}\ge 0\) , respectively, and \(A_i\in {R}^{{n_\mathrm {y}}\times {n_\mathrm {y}}}\)  and \(B_j\in {R}^{{n_\mathrm {y}}\times {n_\mathrm {u}}}\)  are scheduling-dependent matrix coefficient functions that are affine in the basis functions \(\lbrace \alpha _{i,l} \rbrace _{i=1,l=1}^{{n_\mathrm {a}},n_{\alpha }}\) , \(\lbrace \beta _{j,l} \rbrace _{j=0,l=1}^{{n_\mathrm {b}},n_{\beta }}\) , respectively, similarly to (REF ). The LPV-IO representation (REF ) can be created by using the |lpvio| command.


LPV-SS representation

The system can also be represented by a first-order differential or difference equation, i.e., state-space representation:
\(\xi x_t &= ({A} \diamond p)_t\,&x_t&+({B}\diamond p)_t&&u_t,  \\y_t &= ({C}\diamond p)_t&x_t&+({D}\diamond p)_t&&u_t. \) 


where \(x:\mathbb {T}\rightarrow \mathbb {X}=\mathbb {R}^{{n_\mathrm {x}}}\)  is the so-called state (latent) variable and \({A}\in {R}^{{n_\mathrm {x}}\times {n_\mathrm {x}}}\) , \({B}\in {R}^{{n_\mathrm {x}}\times {n_\mathrm {u}}}\) , \({C}\in {R}^{{n_\mathrm {y}}\times {n_\mathrm {x}}}\) , and \({D}\in {R}^{{n_\mathrm {y}}\times {n_\mathrm {u}}}\)  are matrix coefficient functions. For the CT case, \(\xi =\tfrac{\mathrm {d}}{\mathrm {d}t}\)  and, for the DT case, \(\xi = q\)  is the forward time-shift operator. The LPV-SS representation (REF ) can be created by using the |lpvss| command.


LPV-LFR representation

An LPV-LFR representation can be used if the dynamics of the system can be decomposed as the linear fractional transform (LFT) of an LTI system \(G\)  with a parameter-varying matrix function \(\Delta \) , as shown in Figure REF . The input-output behavior of the system can then be described by the following equations:
\(\xi x_t &= {A}\,x_t +{B}_\mathrm {w} w_t && + {B}_\mathrm {u} u_t,  \\z_t &= {C}_\mathrm {z} x_t +{D}_\mathrm {zw} w_t && + {D}_\mathrm {zu} u_t,  \\y_t &= {C}_\mathrm {y} x_t +{D}_\mathrm {yw} w_t && + {D}_\mathrm {yu} u_t,  \\w_t &= (\Delta \diamond p)_t z_t.\) 


The LPV-LFR representation can be created by using the |lpvlfr| command. The LPV-LFR representation is a generalization of the LPV-SS representation, as can be seen by setting \({D}_\mathrm {zw} = 0\) . Therefore, LPVcore internally treats LPV-SS models as LPV-LFR models, and every operation defined on LPV-LFR models, such as multiplication, concatenation and feedback interconnection, can directly be used on LPV-SS models. For a complete list of supported operations, see |www.lpvcore.net|.
<FIGURE>

Identification of LPV models

LPVcore includes several system identification methods to estimate DT LPV-IO and LPV-SS models in a global setting, and one method to estimate LPV-LFR models in a local setting. For the global approach, the estimation of the model parameters is performed by minimizing the prediction-error. Hence, in Section REF , the prediction-error minimization (PEM) setting is briefly introduced. Next, the LPV-IO model set and the corresponding identification methods are discussed in Section REF . In Section REF , the model set and PEM for LPV-SS models are discussed. For the local approach, to identify LPV-LFR models, the estimation is based on minimizing the worst-case local \(\mathcal {H}_{\infty }\) -error between the model and a given set of LTI models at the measured operating points. This approach is discussed in Section REF .

The PEM identification setting

This section provides a short overview of PEM, for a more detailed description, see e.g. [6]}. In PEM, model estimates are obtained by minimizing their associated prediction-error with respect to a loss function where the (unknown) model parameters are the optimization variables. Hence, the concept is to obtain a parameter estimate \(\hat{\theta }_N\)  of the parameter \(\theta _\mathrm {o}\)  describing the data-generating system by minimizing the least-squares criterion, i.e., by
\(\hat{\theta }_N=\operatornamewithlimits{\arg \!\min ~}_{\theta \in \mathbb {R}^{n_{\theta }}}V({D}_N,\theta ), \) 
with
\(V({D}_N,\theta )=\frac{1}{N}\sum _{t=1}^{N} \Vert \varepsilon _{t | \theta }\Vert  _2^2= \frac{1}{N} \Vert  \varepsilon _\theta \Vert  _{\ell _2}^2,\) 
where \({D}_N= \lbrace u_t, p_t, y_t\rbrace _{t=1}^N\)  is the observed data sequence of the system and the one-step-ahead prediction-error is \(\varepsilon _{t | \theta }=y_t-\hat{y}_{t | \theta , t-1}\) . In Sections REF  and REF , the one-step-ahead prediction-error \(\varepsilon _{t | \theta }\)  associated with the IO and SS model sets will be introduced.
In system identification, the data \({D}_N\)  is assumed to be corrupted with noise. Hence, it is essential to extend the system representations with a noise process. These extensions are discussed in the next sections.


LPV-IO model identification

Model set:
The IO representation of the data-generating system is captured by the parametrized model \(\mathcal {M}_\theta \) :
\((\mathcal {F}(q^{-1},\theta )\diamond p)_t\ \breve{y}_t &= (\mathcal {B}(q^{-1},\theta )\diamond p)_t\,q^{-\tau _\mathrm {d}} \ u_t,  \\(\mathcal {D}(q^{-1},\theta )\diamond p)_t\ v_t&= (\mathcal {C}(q^{-1},\theta )\diamond p)_t\ \varepsilon _{t|\theta }, \\(\mathcal {A}(q^{-1},\theta )\diamond p)_t\ y_t &=\breve{y}_t+v_t,\) 


where \(\mathcal {A}(\centerdot ),\ldots ,\mathcal {D}(\centerdot ),\mathcal {F}(\centerdot )\)  are polynomials in \(q^{-1}\)  with scheduling-dependent coefficients, which are parametrized in \(\theta \)  using a-priori given basis functions.
The model (REF ) with \(\mathcal {A}\triangleq I_{{n_\mathrm {y}}}\)  is also known as the LPV Box-Jenkins (BJ) model. By considering \(\mathcal {C}\triangleq \mathcal {D} \triangleq \mathcal {F}\triangleq I_{{n_\mathrm {y}}}\)  in (REF ), the LPV version of the so-called auto regressive with exogenous input (ARX) model structure is obtained. Similarly, the auto regressive moving average with exogenous input (ARMAX) model is found by considering \(\mathcal {D}\triangleq \mathcal {F}\triangleq I_{{n_\mathrm {y}}}\) , and the output-error (OE) model by \(\mathcal {A}\triangleq \mathcal {C}\triangleq \mathcal {D}\triangleq I_{{n_\mathrm {y}}}\) .
Each of the above model structures is available in LPVcore using the |lpvidpoly| object, which serves as an LPV extension to the |idpoly| object included in the System Identification Toolbox. A powerful feature of |lpvidpoly| is the ability to exploit known sparsity of \(\theta \) , as demonstrated in the following example.
[LPV-ARX model]
Consider an LPV-ARX model of the form \(A(p,{\theta }) q^{-1} y + y = B{(\theta )} u\) :
\(A(p,{\theta }) = \begin{pmatrix} \theta _1 & \theta _2 \\ 0 & \theta _3 \end{pmatrix} p, \quad B{(\theta )} = \begin{pmatrix} \theta _4 \\ 0 \end{pmatrix}.\) 
This model structure can be created as follows:

p = preal('p', 'dt');
A = [1, 1; 0, 1] * p;
B = [1; 0];
sys = lpvidpoly(eye(2), A, B);


The object |sys| is then used as a template model structure in subsequent LPV-IO identification methods, which use it to extract the scheduling dependence and, in some methods, the initial parameter values. Zero elements are, by default, excluded from the set of identifiable parameters, providing an intuitive method for semi-grey box identification of LPV-IO representations. This behaviour can be disabled: see documentation at |www.lpvcore.net|.

Identification schemes:
The toolbox contains four time-domain methods for identifying IO model structures (REF ): 1) linear regression, 2) pseudo-linear regression, 3) gradient-based search, and 4) instrumental variable method. These methods are based on [7]}, [6]}, [9]}. Linear regression is applicable for the ARX model set and the model estimate is obtained by executing |lpvarx|. For the |lpvarx| method, regularization is included by solving a Tikhonov regression problem, where various weighting matrices can automatically be tuned using generalised cross-validation [10]} or marginal likelihood optimization [11]}. These options can be selected using the option set |lpvarxOptions|. The pseudo-linear regression algorithm is executed for the ARMAX, OE, and BJ model sets by calling |lpvarmax|, |lpvoe|, and |lpvbj|, respectively. The gradient-based search is applicable for all afore-mentioned model types and is executed by calling |lpvpolyest|. The instrumental variable method is applied by the command |lpviv|. At the moment, this function is available for model structures with a single input and a single output (SISO), and will be extended to the general case in a later release.


LPV-SS model identification

Model set:
The model set for LPV-SS representations supports innovation noise, i.e.:
\(x_{t+1} &\!=\! ({A}(\theta )\diamond p)_t&x_t&\!+\!({B}(\theta )\diamond p)_t&&u_t\!+\!({K}(\theta )\diamond p)_t e^\mathrm {o}_t,  \\y_t &\!=\! ({C}(\theta )\diamond p)_t&x_t&\!+\!({D}(\theta )\diamond p)_t&&u_t\!+\!e^\mathrm {o}_t, \) 


where the matrix coefficient functions \({A},\ldots ,{D},{K}\)  are affine in the set of a-priori specified basis functions \(\left\lbrace \psi _i\diamond p\right\rbrace _{i=1}^{n_{\psi }}\)  similar to (REF ). The sequence \(e^\mathrm {o}_t\)  is the sample path realization of the zero-mean stationary process:
\( \mathbf {e}^\mathrm {o}_t\sim \mathcal {N}(0,\Xi ),\) 
with \(\Xi \)  a positive-definite, symmetric real matrix.

Identification schemes:
To identify the state-space model with parametrization (REF ), the toolbox includes an iterative, gradient-based (GB) optimization technique which obtains an estimate of the innovation model. In a future release, a second iterative optimization technique which is based on the expectation-maximization (EM) ([12]}) method will be included. The |lpvssest| function estimates an LPV-SS model with the GB algorithm. Since this method is based on a non-linear optimization, an initial estimate must be provided by the user. This initial estimate can be obtained by identifying an IO model and transforming it into the state-space form (available in a future release of LPVcore) or by subspace identification which is the default option [13]}, [3]}.
The GB method includes several steps to efficiently estimate the parametric LPV-SS model with innovation noise (REF ). Among others, the implementation ensures that the non-linear optimization does not wander among parametrizations of SS models with equivalent input-scheduling-output behavior. The method is explained in detail in [9]}.


LPV-LFR model identification

The local identification scheme for LPV-LFR representations included in LPVcore is based on [16]}. Due to a lack of space, a detailed explanation can be found on |www.lpvcore.net| soon.


Using identified models for control synthesis

LPVcore includes several control synthesis algorithms. LPV-SS models identified using the previously introduced schemes can directly be used for control by synthesizing a stabilizing controller based on given performance specifications in a variety of settings. For a complete overview of the available algorithms, consult the documentation which is included in LPVcore.


Simulation study

Non-linear data-generating system

To illustrate the effectiveness of LPVcore, several of the implemented identification methods are applied to estimate a model of a DC motor with an unbalanced disc. The system consists of a vertically mounted disc with a lumped mass \(m\)  at distance \(l\)  from the center, actuated by an input voltage \(u\) . Due to the gravitational force exerted on the lumped mass, the system exhibits non-linear dynamics, which can be represented in the following differential equation:
\(\ddot{\theta } = -\frac{1}{\tau } \dot{\theta } + \frac{K_\mathrm {m}}{\tau } u - \frac{m g l}{J} \sin (\theta )\) 
where \(\theta \)  is the disc angular position, \(\tau \)  is the lumped back EMF constant, \(K_{m}\)  is the motor constant, \(J\)  is the complete disc inertia and \(g\)  is the acceleration due to gravity. The goal of this section is to identify the dynamics of (REF ) based on simulation data. The simulated output is the angular position \(\theta \)  corrupted with additive white Gaussian noise \(\varepsilon \)  of a specified signal-to-noise ratio SNR.


LPV identification

An LPV embedding of the non-linear dynamics in (REF ) can be obtained by introducing the scheduling signal \(p\) :
\(p := \operatorname{\mathrm {sinc}}(\theta ) = \frac{\sin (\theta )}{\theta }.\) 
Substitution of (REF ) into (REF ) leads to the following LPV-IO representation of the dynamics of the unbalanced disc:
\(\ddot{\theta } = -\frac{1}{\tau } \dot{\theta } - \frac{m g l}{J} p \theta + \frac{K_\mathrm {m}}{\tau } u\) 
Global identification approaches available in LPVcore operate in DT. Therefore, the LPV system (REF ) is discretized using forward Euler, i.e.:
\(\dot{\theta }_k := \dot{\theta }(k T_\mathrm {s}) \approx \frac{\theta _{k+1} - \theta _{k}}{T_\mathrm {s}},\) 
with \(T_\mathrm {s}=0.75\)  ms the sampling time. Applying (REF ) to (REF ) yields the following DT approximation:
\(\theta _{k} + A_1 \theta _{k-1} +\left( A_2 \diamond p \right)_k \theta _{k-2}= \frac{K_\mathrm {m} T_\mathrm {s}^2}{\tau } u_{k-2},\) 
with
\(A_1 = \frac{T_\mathrm {s}}{\tau } - 2, \quad (A_2 \diamond p)_k = 1 - \frac{T_\mathrm {s}}{\tau } + \frac{m g l T_\mathrm {s}^2}{J} p_{k-2}. \\\) 
Figure REF  shows the quality of the approximation in terms of the Best Fit Rate (BFR) for a simulation. The BFR, expressed as a percentage, is defined as follows:
\(\mathrm {BFR} := 100\% \cdot \max \left( 0, 1 - \frac{\left\Vert y - \hat{y}\right\Vert _2}{ \left\Vert  y - \mathrm {mean}(y)\right\Vert } \right).\) 
The values of all numerical parameters used in this example are listed in Table REF .
<TABLE><FIGURE>For the identification procedure, five datasets of 400 samples each are generated: four for estimation with varying levels of additive White Gaussian output noise and one for validation (without noise). These datasets are generated by simulating (REF ) for a random-phase multi-sine input \(u\)  containing 10 frequencies equally spaced in a passband from 0 to 0.75 times the Nyquist frequency. The signal is scaled to have a maximum amplitude of 0.25 V. Zero-order-hold is used to interpolate between the samples. The SNRs for the estimation datasets are 0, 10, 20 and 40 dB. A noise-free example of the estimation dataset and a comparison with the DT LPV-IO approximation (REF ) can be seen in Figure REF . Four LPV-IO model structures are identified: ARX, ARMAX, OE and BJ. It is assumed that the true model structure (REF ) is not known exactly and the following more generic template structure is used instead (see (REF )):

\((\mathcal {A}(q^{-1}) \diamond p)_k = 1 + \left( 1 + p_{k-1} \right) q^{-1} + \left( 1 + p_{k-2} \right) q^{-2}\)

\((\mathcal {B}(q^{-1}) \diamond p)_k = (\mathcal {F}(q^{-1}) \diamond p)_k = (\mathcal {A}(q^{-1}) \diamond p)_k\)

\(\mathcal {C}q^{-1}) = \mathcal {D}(q^{-1}) = 1 + q^{-1}\)

The models are identified using the estimation dataset in the following order:
The LPV-ARX model is identified by calling |lpvarx|.

The LPV-ARMAX and LPV-OE models are identified by calling |lpvpolyest|. The LPV-ARX model is used as initialization.

The LPV-BJ model is identified by calling |lpvpolyest|. The LPV-OE model is used as initialization.

For each method, the default options were used, except the number of gradient descent iterations performed by |lpvpolyest|: it was set to 400. The quality of the resulting models is evaluated using the validation dataset by calculating the BFR between the output of the non-linear CT system (REF ) and the output of the DT model estimate.

Results

The achieved BFR of the model estimates can be seen in Figure REF . For each SNR, the LPV-OE and LPV-BJ models have a substantially higher BFR than the LPV-ARX and LPV-ARMAX models, which is in line with the expectation since output noise was used. The low BFR of the LPV-ARX and LPV-ARMAX models, even for high SNRs, is likely the result of a high parameter sensitivity due to the proximity of the poles of the frozen dynamics of (REF ) to the unit circle. For example, at \(p = 0\)  (\(\theta = \pm \pi \)  rad), the magnitude of the two poles of the resulting LTI dynamics are 1 and 0.92. When \(p = 1\)  (\(\theta = 0\) ), both poles have a magnitude of 0.96. The output of the LPV-ARX and LPV-OE models estimated with SNR = 20 dB are compared to the original non-linear CT dynamics (REF ) in Figure REF .
<FIGURE><FIGURE>

Conclusion

This paper introduces a software package for MATLAB called the LPVcore toolbox to model, simulate, and estimate linear parameter-varying input-output, state-space and linear fractional representations in a global setting. To support the global setting, the package defines parameter-varying matrix functions, which is the key difference to other software suites that use a grid or linear fractional based representation. This paper provides an overview of the core modelling tools combined with the simulation and identification methods. The software and examples are freely available on |www.lpvcore.net|.
The initial toolbox forms a solid basis for additional algorithms. In the near future, algorithms for LPV realization and a wide range of subspace identification methods will be implemented.

