# Improved Parameter Estimation Techniques for Induction Motors using Hybrid
Algorithms

Julius Susanto, and Syed Islam Corresponding author: J. Susanto, e-mail:
<EMAIL_ADDRESS>
###### Abstract

The performance of Newton-Raphson, Levenberg-Marquardt, Damped Newton-Raphson
and genetic algorithms are investigated for the estimation of induction motor
equivalent circuit parameters from commonly available manufacturer data. A new
hybrid algorithm is then proposed that combines the advantages of both descent
and natural optimisation algorithms. Through computer simulation, the hybrid
algorithm is shown to significantly outperform the conventional algorithms in
terms of convergence and squared error rates. All of the algorithms are tested
on a large data set of 6,380 IEC (50Hz) and NEMA (60Hz) motors.

##  I Introduction

The three-phase induction motor is arguably the workhorse of modern industry,
found in almost all industrial settings from manufacturing to mining.
Equivalent circuit parameters of induction machines are essential for time-
domain simulations where the dynamic interactions between the machine(s) and
the power system need to be analysed, for example:

  * •

Motor starting and re-acceleration

  * •

Bus transfer studies

  * •

Changes in motor loading

  * •

Motor behaviour during faults

  * •

Dynamic voltage stability

However, motor manufacturers do not tend to provide the equivalent circuit
parameters for their machines. As the parameters are motor specific, typical
values found in the literature are often not sufficiently accurate. Moreover,
power system studies involving motors are normally performed during the design
phases of projects, where the motors themselves have not yet been ordered and
on-site testing is not possible.

It is therefore desirable to estimate motor equivalent circuit parameters from
the data that manufacturers make available in their catalogues, data sheets
and technical brochures (e.g. performance parameters such as breakdown torque,
locked rotor torque, full-load power factor, full-load efficiency, etc).

A number of parameter estimation techniques have been proposed in the
literature based on manufacturer data (for example [2], [3], [4], [5] and
[6]). The de facto approach that has emerged, and which has been adopted by
the majority of commercial software packages, has been to use an algorithm
based on the Newton-Raphson method. However, it has been observed that the
Newton-Raphson based algorithms can have poor convergence and error
performance. Therefore, parameter estimation algorithms with improved
performance are desired.

In this article, the performance of a number of parameter estimation
algorithms based on readily available manufacturer data is investigated for
the double-cage induction motor model. Section II describes the double-cage
motor model and the formulation of the parameter estimation problem. In
Section III, descent algorithms such as the Newton-Raphson method are
discussed, followed by an overview of natural optimisation algorithms in
Section IV. The properties of descent and natural optimisation algorithms are
combined to overcome some of the limitations of descent algorithms and a new
hybrid algorithm is proposed in Section V. The algorithms are then tested via
computer simulation in Section VI and finally, the conclusions of this study
are presented in Section VII.

##  II Parameter Estimation Problem for Induction Motors

###  II-A Double Cage Induction Motor Model

It has been previously shown that the double-cage model is appropriate for
squirrel-cage induction motors in order to capture both the starting and
breakdown performance of the motor [1] [6].

The double cage equivalent circuit model with eight slip-invariant parameters
valid over the full range of slip values (i.e. from 0 to 1 pu) is shown in
Figure 1.

Figure 1: Double cage equivalent circuit (eight parameter model)

In the equivalent circuit, the inner cage leakage reactance $X_{r1}$ is always
higher than the outer cage leakage reactance $X_{r2}$, but the outer cage
impedance is typically higher than the inner cage impedance on starting. These
conditions can be resolved by including the following two inequality
constraints in the model [6]:

  * •

$X_{r1}>X_{r2}$

  * •

$R_{r2}>R_{r1}$

In order to estimate motor efficiency, the core (and mechanical) losses also
need to be included in the model. For simplicity, the core (and mechanical)
losses are represented as a single shunt resistance $R_{c}$ at the input of
the equivalent circuit [8].

###  II-B Parameter Estimation Problem Formulation

The characteristics of an induction motor are normally provided by
manufacturers in the form of standard performance parameters, with the
following parameters being the most common:

  * •

Nominal voltage, $U_{n}$ (V)

  * •

Nominal frequency, $f$ (Hz)

  * •

Rated asynchronous speed, $n_{fl}$ (rpm)

  * •

Rated (stator) current, $I_{s,fl}$ (A)

  * •

Rated mechanical power, $P_{m,fl}$ (kW)

  * •

Rated torque, $T_{n}$ (Nm)

  * •

Full load power factor, $\cos{\phi_{fl}}$ (pu)

  * •

Full load efficiency, $\eta_{fl}$ (pu)

  * •

Breakdown torque, $T_{b}/T_{n}$ (normalised)

  * •

Locked rotor torque, $T_{lr}/T_{n}$ (normalised)

  * •

Locked rotor current, $I_{lr}/I_{s,fl}$) (pu)

Given a set of performance parameters that contain features on the torque-
speed and current-speed curves (e.g. breakdown torque, locked-rotor current,
etc), the goal of a parameter estimation algorithm is to determine equivalent
circuit parameters that yield these features.

While all of the performance parameters in the list above can be used in the
estimation procedure, there are actually only six indpendent magnitudes that
can be formed from them: $P_{m,fl}$, $Q_{fl}$, $T_{b}$, $T_{lr}$, $I_{lr}$ and
$\eta_{fl}$ [6].

The six independent magnitudes can be used to formulate the parameter
estimation problem in terms of a non-linear least squares problem, with a set
of non-linear equations of the form $\mathbf{F}(\boldsymbol{x})=\mathbf{0}$:

$\displaystyle f_{1}(\boldsymbol{x})=P_{m,fl}-P(s_{f})=0$ (1) $\displaystyle
f_{2}(\boldsymbol{x})=Q_{fl}-Q(s_{f})=0$ (2) $\displaystyle
f_{3}(\boldsymbol{x})=T_{b}-T(s_{max})=0$ (3) $\displaystyle
f_{4}(\boldsymbol{x})=T_{lr}-T(s=1)=0$ (4) $\displaystyle
f_{5}(\boldsymbol{x})=I_{lr}-I(s=1)=0$ (5) $\displaystyle
f_{6}(\boldsymbol{x})=\eta{fl}-\eta(s_{f})=0$ (6)

where $\mathbf{F}=(f_{1},f_{2},f_{3},f_{4},f_{5},f_{6})$ and
$\boldsymbol{x}=(R_{s},X_{s},X_{m},R_{r1},X_{r1},R_{r2},X_{r2},R_{c})$

In the formulation above, there are six independent equations in $\mathbf{F}$,
but eight unknown parameters in $\boldsymbol{x}$. The non-linear system is
therefore underdetermined, which leads to one of the key limitations in
descent algorithms (refer to Section III for further discussion).

###  II-C Classes of Parameter Estimation Algorithms

The parameter estimation problem formulated previously can be solved by a
variety of non-linear least squares solver algorithms. As is usually the case
with non-linear least squares problems, a closed form solution is not
available and iterative algorithms are used to converge on a solution by
minimising error residuals.

Induction motor parameter estimation algorithms generally fall under two broad
classes:

  1. 1.

Descent Methods: are the class of algorithms based on variations of Newton’s
method for convergence to a solution, e.g. Newton-Raphson, Levenberg-
Marquardt, etc

  2. 2.

Natural Optimisation Methods: are the class of algorithms based on processes
found in nature where successive randomised trials are filtered for ”fitness”
at each iteration, e.g. genetic algorithm, particle swarm optimisation, ant
colony optimisation, simulated annealing, etc

These classes of parameter estimation algorithms are discussed in more detail
in the following sections.

##  III Descent Algorithms

###  III-A Requirement for Linear Restrictions

It was mentioned in Section II that the non-linear system of equations in the
problem forumulation is underdetermined. Therefore, in order to make the
systems exactly determined and solvable with descent algorithms, one must
either:

  1. 1.

Fix two parameters a priori (i.e. parameters are ”known”), or;

  2. 2.

Impose two constraints on the problem formulations, e.g. linear restrictions

The use of linear restrictions was found to be superior to fixed parameters
and therefore, all of the descent algorithms in this paper will include two
linear restrictions by default.

It was shown in [6] that the stator resistance $R_{s}$ was the least sensitive
parameter in the equivalent circuit, i.e. variations in the value of $R_{s}$
had the least significant effect on the resulting torque-speed and current-
speed curves. Therefore, $R_{s}$ can be subject to a linear restriction by
linking it to the inner cage rotor resistance $R_{r1}$. Moreover, it is
assumed that the outer cage rotor reactance $X_{r2}$ is linearly related to
the stator reactance $X_{s}$. This leads to the following linear restrictions:

  * •

$R_{s}=k_{r}R_{r1}$

  * •

$X_{r2}=k_{x}X_{s}$

Where $k_{r}$ and $k_{x}$ are constant linear restrictions

###  III-B Newton-Raphson Algorithm

Of the class of descent methods used to solve non-linear least squares
problems, the Newton-Raphson (NR) algorithm is probably the most
straightforward. The NR algorithm is an iterative method where each iteration
is calculated as follows:

$\boldsymbol{x}^{k+1}=\boldsymbol{x}^{k}-h_{n}\mathbf{J}^{-1}\mathbf{F}(\boldsymbol{x}^{k})$
(7)

where $\boldsymbol{x}^{k+1}$ is the solution at the $(k+1)$th iteration,
$\boldsymbol{x}^{k}$ is the solution at the $k$th iteration, $h_{n}$ is the
step-size coefficient and $\mathbf{J}$ is the Jacobian matrix evaluated with
the parameters at the $k$th iteration, $\boldsymbol{x}^{k}$.

The Jacobian matrix $\mathbf{J}$ has the general form:

$\mathbf{J}=\left(\begin{array}[]{ccc}\frac{\partial f_{1}}{\partial
x_{1}}&\ldots&\frac{\partial f_{1}}{\partial x_{6}}\\\ \vdots&\ddots&\vdots\\\
\frac{\partial f_{6}}{\partial x_{1}}&\ldots&\frac{\partial f_{6}}{\partial
x_{6}}\\\ \end{array}\right)$ (8)

For systems where it is impractical to compute the exact partial derivatives
analytically, a numerical approximation may be used with finite difference
equations:

$\frac{\partial f_{i}}{\partial
x_{j}}\approx\frac{f_{i}(\mathbf{x}+\boldsymbol{\delta_{j}}h)-f_{i}(\mathbf{x})}{h}$
(9)

where $\boldsymbol{\delta_{j}}$ is vector of zeros with a single non-zero
value of 1 at the j-th element and $h$ is a constant with a very small
absolute value (e.g. $1\times 10^{-6}$).

In this paper, a modified form of the NR algorithm proposed by Pedra in [8]
for the double cage model is used as the default NR solver. This algorithm was
selected because of its completeness, numerical accuracy and robustness
compared to previously proposed methods (for example, in [2], [3] and [4]).
Furthermore, the algorithm can be applied using commonly available
manufacturer data, whereas other algorithms require more detailed data that
may not be readily available (for example, multiple points on the torque-speed
curve in [5]).

###  III-C Levenberg-Marquardt Algorithm

The Levenberg-Marquardt (LM) algorithm, sometimes called the damped least-
squares algorithm, is another popular technique for solving least-squares
problems [9] [10]. In the LM algorithm, each iteration is calculated as
follows:

$\boldsymbol{x}^{k+1}=\boldsymbol{x}^{k}-\left[\mathbf{J}^{T}\mathbf{J}+\lambda\,\mathbf{diag}(\mathbf{J}^{T}\mathbf{J})\right]^{-1}\mathbf{J}^{T}\mathbf{F}(\boldsymbol{x}^{k})$
(10)

where $\boldsymbol{x}^{k+1}$ is the solution at the $(k+1)$th iteration,
$\boldsymbol{x}^{k}$ is the solution at the $k$th iteration, $\lambda$ is the
damping parameter (more on this later) and $\mathbf{J}$ is the Jacobian matrix
evaluated with the parameters at the $k$th iteration, $\boldsymbol{x}^{k}$ (as
described previously in Equation 8).

####  III-C1 Choice of Damping Parameter

The selection of the damping parameter $\lambda$ affects both the direction
and magnitude of an iteration step. If the damping parameter is large, then
the algorithm will move at short steps in the steepest descent direction. This
is good when the iteration is far away from the solution. On the other hand,
if the damping parameter is small, then the algorithm approaches a Gauss-
Newton type method, which exhibits good convergence in the neighbourhood of
the solution.

The damping parameter should thus be updated at each iteration depending on
whether the algorithm is far or close to the solution. Two methods for
adjusting the damping parameter are described below.

##### Gain Ratio Adjustment

Marquardt suggested updating the damping parameter based on a ”gain ratio”
[10]:

$\rho=\frac{\mathbf{F}(\boldsymbol{x}^{k})-\mathbf{F}(\boldsymbol{x}^{k+1})}{\frac{1}{2}\mathbf{\Delta
x}^{T}(\lambda\mathbf{\Delta
x}-\mathbf{J}^{T}\mathbf{F}(\boldsymbol{x}^{k}))}$ (11)

Where $\mathbf{\Delta
x}=-\left[\mathbf{J}^{T}\mathbf{J}+\lambda\,\mathbf{diag}(\mathbf{J}^{T}\mathbf{J})\right]^{-1}\mathbf{J}^{T}\mathbf{F}(\boldsymbol{x}^{k})$
is the correction step at iteration k.

The damping parameter is adjusted depending on the value of the gain ratio as
follows:

$\lambda=\begin{cases}\lambda\times\beta,&\text{if }\rho<\rho_{1}\\\
\frac{\lambda}{\gamma},&\text{if }\rho>\rho_{2}\end{cases}$ (12)

Where $\rho_{1}$, $\rho_{2}$, $\beta$ and $\gamma$ are algorithm control
parameters.

##### Error Term Adjustment

An alternative to using the gain ratio is to adjust the damping parameter
based only on the error term (i.e. the numerator of the gain ratio). The
damping parameter is therefore updated as follows:

$\lambda=\begin{cases}\lambda\times\beta,&\text{if
}\mathbf{F}(\boldsymbol{x}^{k})-\mathbf{F}(\boldsymbol{x}^{k+1})<0\\\
\frac{\lambda}{\gamma},&\text{if
}\mathbf{F}(\boldsymbol{x}^{k})-\mathbf{F}(\boldsymbol{x}^{k+1})>0\end{cases}$
(13)

Where $\rho_{1}$, $\rho_{2}$, $\beta$ and $\gamma$ are algorithm control
parameters.

Computer simulations showed that the error term adjustment yielded
significantly better performance than the gain ratio adjustment. Therefore in
these studies, the error term adjustment is used exclusively along with the
following algorithm control parameters: $\rho_{1}=0.25$, $\rho_{2}=0.75$,
$\beta=3$ and $\gamma=3$.

###  III-D Damped Newton-Raphson Algorithm

The damped Newton-Raphson algorithm is a variation of the conventional NR
algorithm where a damping factor helps to get around problems with near-
singular and/or ill-conditioned Jacobian matrices. In the damped NR algorithm,
each iteration is calculated as follows:

$\boldsymbol{x}^{k+1}=\boldsymbol{x}^{k}-h_{n}(\mathbf{J}^{-1}+\lambda
I)\mathbf{F}(\boldsymbol{x}^{k})$ (14)

Where the damping parameter $\lambda$ is adjusted at each iteration based on
the error term as follows:

$\lambda=\begin{cases}\lambda\times\beta,&\text{if
}\mathbf{F}(\boldsymbol{x}^{k})-\mathbf{F}(\boldsymbol{x}^{k+1})<0\\\
\frac{\lambda}{\gamma},&\text{if
}\mathbf{F}(\boldsymbol{x}^{k})-\mathbf{F}(\boldsymbol{x}^{k+1})>0\end{cases}$
(15)

##  IV Natural Optimisation Algorithms

###  IV-A Genetic Algorithm

In the context of motor parameter estimation, the genetic algorithm (GA) is
used to minimise the squared error of the problem formulation vector
$\mathbf{F}$. In GA terminology, the squared error is called the ”fitness
function” and is calculated as follows:

$fitness=\mathbf{F}\mathbf{F}^{\prime}$ (16)

where $\mathbf{F}=(f_{1},f_{2},f_{3},f_{4},f_{5},f_{6})$

Genetic algorithms can be binary coded where the solution paramaters are
quantized into binary strings (for example, in [11], [12] and [13]). However,
the equivalent circuit parameters in a motor are continuous parameters and not
naturally quantized. Thus, binary coding necessarily imposes limits on the
precision of the parameters (i.e. due to the chosen length of the binary
string). For this reason, a continuous parameter algorithm is used in this
study. A description of the genetic algorithm as it applies to motor parameter
estimation follows.

An initial population of $n_{pop}$ parameter estimates are randomly sampled
from a uniform distribution with upper and lower limits as shown in Table I.

TABLE I: Range of initial parameter estimates Parameter | Range of Initial Estimate (pu)
---|---
Lower Bound | Upper Bound
$R_{s}$ | 0 | 0.15
$X_{s}$ | 0 | 0.15
$X_{m}$ | 0 | 5
$R_{r1}$ | 0 | 0.15
$X_{r1}$ | 0 | 0.30
$R_{r2}$ | 0 | 0.15
$X_{r2}$ | 0 | 0.15
$R_{c}$ | 0 | 100

The fitness of each member in the population is then calculated and ranked.
The lowest fitness members are discarded and the rest are retained to form the
mating pool for the next generation (there are $n_{pool}$ members in the
mating pool).

The fittest $n_{e}$ members in the mating pool are retained for the next
generation as elite children.

Of the remaining $n_{pop}-n_{e}$ children to be created for the next
generation, $c_{f}$% will be produced by crossover and the rest ($1-c_{f}$%)
by mutation. The proportion $c_{f}$ is called the crossover fraction.

  1. 1.

Crossover: in the crossover process, two members of the mating pool are
randomly selected and combined by taking a random blend of each member’s
parameters.

  2. 2.

Mutation: in the mutation process, a member of the mating pool is randomly
selected and its parameters are mutated by adding Gaussian noise with
parameter-dependent standard deviations (see Table II).

TABLE II: Standard deviations for mutation noise Parameter | Standard Deviation ($\sigma$)
---|---
$R_{s}$ | 0.01
$X_{s}$ | 0.01
$X_{m}$ | 0.33
$R_{r1}$ | 0.01
$X_{r1}$ | 0.01
$R_{r2}$ | 0.01
$X_{r2}$ | 0.01
$R_{c}$ | 6.67

The fitness of the next generation is then calculated and the process repeats
itself for $n_{gen}$ generations.

The default settings for the genetic algorithm implemented in this project for
motor parameter estimation are shown in Table III.

TABLE III: Default settings for genetic algorithm Setting | Setting Description | Default Value
---|---|---
$n_{pop}$ | Population of each generation | 20
$n_{pool}$ | Number of members in the mating pool | 15
$n_{e}$ | Number of elite children | 2
$c_{f}$ | Crossover fraction | 80%

###  IV-B Other Natural Optimisation Algorithms

Since the popularisation of the genetic algorithm in the 1980s, a number of
other natural optimisation algorithms have been introduced. All of these
algorithms are inspired by natural processes and all share the basic
methodology of injecting randomness and selecting for fitness in order to
iteratively reach an optimal point in the search-space.

The literature has numerous examples of natural optimisation algorithms
adapted for the estimation of induction motor parameters, for example:

  * •

Simulated annealing [14]

  * •

Particle swarm optimisation [15] [16]

  * •

Artificial immune system [17]

  * •

Bacterial foraging technique [18]

  * •

Ant colony optimisation [19]

  * •

Harmony search [20]

The results of these investigations indicate that while some of the
alternative algorithms lead to faster convergence than the genetic algorithm,
they tend to converge to the same solution and error rates. The aim of this
study is to improve error rates and for this reason, alternative natural
optimisation algorithms are not explored further.

##  V Proposed Hybrid Algorithm

In the descent algorithms, linear restrictions are imposed on $R_{s}$ and
$X_{r2}$ in order to make the underdetermined system of equations solvable. It
was shown in [7] that the double cage model with core losses has 8 minimum
independent variables (MIVs). Therefore, by constraining $R_{s}$ and $X_{r2}$
with linear restrictions, the solution space is also constrained by two
degrees of freedom. Therefore, a solution could potentially exist to an
otherwise non-converging problem were the linear restrictions removed.

The proposed hybrid algorithm attempts to overcome this limitation of descent
algorithms by applying a genetic algorithm to select $R_{s}$ and $X_{r2}$. In
other words, a baseline descent algorithm (e.g. NR, Damped NR, LM, etc) is run
with fixed values for $R_{s}$ and $X_{r2}$, which are in turn iteratively
selected using a genetic algorithm in an outer loop. A flowchart of the
proposed hybrid algorithm is shown in Figure 2. A more detailed description of
the proposed algorithm follows.

Figure 2: Flowchart for hybrid algorithm (with natural selection of $R_{s}$
and $X_{r2}$)

An initial population of $n_{pop}$ estimates for $R_{s}$ and $X_{r2}$ are
randomly sampled from a uniform distribution with upper and lower limits as
shown in Table IV. Each pair of $R_{s}$ and $X_{r2}$ estimates is referred to
as a member of the population.

TABLE IV: Range of initial parameter estimates Parameter | Range of Initial Estimate (pu)
---|---
Lower Bound | Upper Bound
$R_{s}$ | 0 | 0.15
$X_{r2}$ | 0 | 0.15

The descent algorithm is then run on each member of the population. The
fitness of each member (in terms of the squared error
$\mathbf{F}^{\prime}\mathbf{F}$) is calculated and ranked. As in the
conventional genetic algorithm, the lowest fitness members are discarded and
the rest are retained to form the mating pool for the next generation.

The fittest $n_{e}$ members in the mating pool are retained for the next
generation as elite children. Of the remaining $n_{pop}-n_{e}$ children to be
created for the next generation, $c_{f}$% will be produced by crossover and
the rest ($1-c_{f}$%) by mutation. The proportion $c_{f}$ is called the
crossover fraction.

  1. 1.

Crossover: in the crossover process, two members of the mating pool are
randomly selected and combined by taking a random blend of each member’s
parameters, e.g. the crossover of parameter $R_{s}$:

$R_{s,child}=\alpha R_{s,parent1}+(1-\alpha)R_{s,parent2}$ (17)

where $\alpha$ is a random variable selected from a uniform distribution over
the interval $[0,1]$

  2. 2.

Mutation: in the mutation process, a member of the mating pool is randomly
selected and its parameters are mutated by adding Gaussian noise with standard
deviations of 0.01.

The descent algorithm is then run for the next generation of estimates for
$R_{s}$ and $X_{r2}$. The fitness is calculated and the process repeats itself
for $n_{gen}$ generations. If at any point during the process the descent
algorithm converges, then the hybrid algorithm stops and selects the parameter
estimates from the converged descent algorithm as the solution. Otherwise, the
parameter estimates yielding the best fitness after $n_{gen}$ generations are
selected.

TABLE V: Default settings for proposed hybrid algorithms Setting | Setting Description | Default Value
---|---|---
$n_{pop}$ | Population of each generation | 15
$n_{pool}$ | Number of members in the mating pool | 10
$n_{e}$ | Number of elite children | 2
$c_{f}$ | Crossover fraction | 80%
$n_{gen}$ | Maximum number of generations | 10

The following hybrid algorithms are investigated in this paper:

  * •

NR-GA - Conventional NR and GA selection of $R_{s}$ and $X_{r2}$

  * •

LM-GA - Levenberg-Marquardt (with an error term lambda adjustment) and GA
selection of $R_{s}$ and $X_{r2}$

  * •

DNR-GA - Damped NR (with a maximum number of 30 iterations for each damped NR
step) and GA selection of $R_{s}$ and $X_{r2}$

The default settings applied for the proposed hybrid algorithms are shown in
Table V.

##  VI Computer Simulations

###  VI-A Motor Data Set

The algorithms presented in this paper were tested by computer simulation on a
large data set from the EuroDEEM and MotorMaster databases (version 1.0.17 - 4
April 2007) [21]. The original data set included over 15,000 IEC and NEMA type
motors. From the original set, the motor data was conditioned by eliminating
duplicate records, removing motors without power factor, efficiency or torque
data and removing motors with strange or inconsistent data (e.g. full load
torque greater than breakdown torque, asynchronous speed greater than
synchronous speed, etc). After data cleansing, the final data set consisted of
motors with nominal ratings from 0.37kW to 1000kW, and the following total
quantities:

  * •

4,002 IEC 50Hz motors

  * •

2,378 NEMA 60Hz motors

###  VI-B Simulation Results

A summary of the computer simulation results is shown in Table VI. The table
depicts the convergence rate, average squared error and maximum squared error
for each of the algorithms described in this paper. Note that in this study,
the convergence criterion is a squared error value of $1\times 10^{-5}$ and
the maximum number of iterations for the descent algorithms is 30.

TABLE VI: Simulation results for all algorithms Case | IEC Motors | NEMA Motors
---|---|---
Convergence | Squared Error |  Convergence | Squared Error
Average | Maximum | Average | Maximum
Newton-Raphson ($k_{x}=1$, $k_{r}=0.5$) |  685 (17.1%) | 0.5411 | 5.9996 |  751 (31.6%) | 0.2514 | 5.9902
Levenberg-Marquardt ($k_{x}=1$, $k_{r}=0.5$) |  740 (18.5%) | 0.9114 | 6 |  770 (32.4%) | 0.2867 | 6
Damped NR ($k_{x}=1$, $k_{r}=0.5$) |  628 (15.7%) | 0.2058 | 9.2314 |  568 (23.9%) | 0.0899 | 4.6188
Newton-Raphson ($k_{x}=0.5$, $k_{r}=1$) |  974 (24.3%) | 0.9261 | 120.08 |  934 (39.3%) | 0.1425 | 5.999
Levenberg-Marquardt ($k_{x}=0.5$, $k_{r}=1$) |  1035 (25.9%) | 2.691 | 6 |  945 (39.7%) | 1.7268 | 6
Damped NR ($k_{x}=0.5$, $k_{r}=1$) |  1006 (25.1%) | 0.04 | 4.4584 |  935 (39.3%) | 0.05054 | 4.6574
Genetic Algorithm (Max gens = 30) |  0 (0.0%) | 0.0471 | 3.6459 |  0 (0.0%) | 0.04029 | 0.50916
Genetic Algorithm (Max gens = 50) |  0 (0.0%) | 0.0281 | 2.031 |  0 (0.0%) | 0.03279 | 0.45363
Genetic Algorithm (Max gens = 100) |  0 (0.0%) | 0.01861 | 0.8062 |  0 (0.0%) | 0.02676 | 0.37716
Hybrid NR-GA |  1363 (34.1%) | 0.0625 | 5.9998 |  1159 (48.7%) | 0.0282 | 0.9682
Hybrid LM-GA |  1388 (34.7%) | 1.0941 | 6 |  1181 (49.7%) | 0.53387 | 6
Hybrid DNR-GA |  1373 (34.3%) | 0.0174 | 5.9995 |  1168 (49.1%) | 0.01948 | 0.42282

In terms of convergence and error rates, it can be seen from Table VI that the
hybrid algorithms are superior to the other algorithms. However, there is no
hybrid algorithm that clearly stands out as the best option. The LM-GA has the
highest number of converging solutions, but when it fails, it yields poor
results (as evidenced by the high average error rate). The DNR-GA is more
consistent in terms of low error rates, but the convergence rate is also
lower.

The genetic algorithms with maximum number of generations $\geq$30 also yield
low average and maximum error rates, but never low enough to qualify for
convergence (as defined by a squared error of $<1\times 10^{-5}$). However,
the performance of the genetic algorithm is unaffected by the choice of linear
restrictions / initial conditions.

On the other hand, the results show that the descent algorithms are
significantly affected by the choice of linear restrictions $k_{r}$ and
$k_{x}$, and yield inferior convergence and error rates when compared to the
proposed hybrid algorithms.

The simulations suggest that the LM algorithm can lead to a higher convergence
rate compared to the conventional NR algorithm, but at the cost of a higher
average squared error. The LM algorithm works well in the neighbourhood of the
solution, but does not perform very well at the early stages, particularly
when the initial estimates are far from the solution. The LM algorithm can
also produce spectacularly bad results when the Jacobian matrix is ill-
conditioned or near-singular.

The damped NR algorithm is intended to help address the issue of ill-
conditioned and near-singular Jacobian matrices. Adding a damping factor
$\lambda I$ to the Jacobian matrix makes it more likely to be invertible.
However, while the damped NR algorithm leads to lower average squared error,
it also takes more iterations to converge.

###  VI-C Algorithm Computation Time

Indicative computation times for the different algorithms in this project are
shown in Table VII. The computation times were obtained from simulations
performed on a 2.1GHz Intel dual core processor with 2GB RAM and are presented
here primarily for comparison.

TABLE VII: Average algorithm solution time  Algorithm | Solution Time (s)
---|---
Average | Maximum
Newton-Raphson | 0.257 | 0.742
Levenberg-Marquardt | 0.162 | 0.332
Damped NR | 0.241 | 0.427
Genetic Algorithm (10 Gens) | 0.257 | 0.328
Genetic Algorithm (30 Gens) | 0.778 | 0.947
Genetic Algorithm (50 Gens) | 1.356 | 1.758
Genetic Algorithm (100 Gens) | 2.753 | 3.036
Hybrid NR-GA | 24.395 | 53.290
Hybrid LM-GA | 14.289 | 42.198
Hybrid DNR-GA | 29.050 | 64.568

From Table VII, it can be seen that the hybrid algorithms have average
solution times between 50 and 100 times slower than the conventional descent
algorithms (i.e. NR, LM and DNR algorithms). This is because the evolutionary
part of the hybrid algorithm must run the descent algorithm multiple times for
each generation. For example, based on the default settings as shown in Table
V, the hybrid algorithm may have to perform up to $n_{pop}\times
n_{gen}=10\times 15=150$ descent algorithms. This would occur in the worst
case condition when the hybrid algorithm fails to converge.

The genetic algorithm has average solution times that are dependent on the
maximum number of generations to be simulated. For a low number of generations
(e.g. 10), the GA solution times are comparable to that of the descent
algorithms. The solution times increase more or less linearly as the maximum
number of generations is increased.

###  VI-D Algorithm Performance and Motor Rated Power

The performance of the algorithms are analysed with the data sets broken down
by motor rated power. Table VIII presents the breakdown of the IEC and NEMA
motor data sets, showing the quantity of motors for various nominal power
ranges.

TABLE VIII: Breakdown of motor data sets by motor rated power Motor Rating | No. IEC Motors | No. NEMA Motors
---|---|---
0.37kW - 3.6kW | 1208 | 630
4kW - 15kW | 963 | 598
18.5kW - 75kW | 973 | 741
90kW - 185kW | 477 | 284
200kW - 630kW | 355 | 123
$>$630kW | 26 | 2
TOTAL | 4002 | 2378
TABLE IX: Algorithm performance broken down by rated power (IEC motors) Case | 0.37 - 3.6kW | 4 - 15kW | 18.5 - 75kW | 90 - 185kW | 200 - 630kW |  $>$630kW
---|---|---|---|---|---|---
Conv-ergence |  Average Error2 |  Conv-ergence |  Average Error2 |  Conv-ergence |  Average Error2 |  Conv-ergence |  Average Error2 |  Conv-ergence |  Average Error2 |  Conv-ergence |  Average Error2
Newton-Raphson ($k_{x}=1$, $k_{r}=0.5$) |  4 (0.33%) |  0.8323 |  52 (0.33%) |  0.4891 |  305 (31.4%) |  0.5880 |  166 (34.8%) |  0.1788 |  137 (38.6%) |  0.0887 |  21 (80.8%) |  0.0092
Levenberg-Marquardt ($k_{x}=1$, $k_{r}=0.5$) |  14 (1.16%) |  1.9423 |  62 (6.44%) |  0.6743 |  336 (34.5%) |  0.1055 |  190 (39.8%) |  0.3229 |  122 (34.4%) |  0.9954 |  16 (61.5%) |  1.6178
Damped NR ($k_{x}=1$, $k_{r}=0.5$) |  4 (0.33%) |  0.5279 |  36 (3.7%) |  0.1132 |  262 (26.9%) |  0.0447 |  160 (33.5%) |  0.0402 |  143 (40.3%) |  0.0401 |  23 (88.5%) |  0.0047
Genetic Algorithm (Max gens = 30) |  0 (0.0%) |  0.0830 |  0 (0.0%) |  0.0196 |  0 (0.0%) |  0.0266 |  0 (0.0%) |  0.0469 |  0 (0.0%) |  0.0550 |  0 (0.0%) |  0.0569
Genetic Algorithm (Max gens = 50) |  0 (0.0%) |  0.0446 |  0 (0.0%) |  0.0157 |  0 (0.0%) |  0.0180 |  0 (0.0%) |  0.0251 |  0 (0.0%) |  0.0362 |  0 (0.0%) |  0.0313
Hybrid NR-GA |  48 (3.97%) |  0.1842 |  245 (25.4%) |  0.0219 |  559 (57.5%) |  0.0035 |  281 (58.9%) |  0.0032 |  207 (58.3%) |  0.0040 |  23 (88.5%) |  0.0014
Hybrid LM-GA |  149 (12.3%) |  2.147 |  274 (28.5%) |  1.196 |  405 (41.6%) |  0.5132 |  303 (63.5%) |  0.1923 |  240 (67.6%) |  0.1185 |  17 (65.4%) |  0.0005
Hybrid DNR-GA |  140 (11.6%) |  0.037 |  265 (27.5%) |  0.014 |  404 (41.5%) |  0.0104 |  303 (63.5%) |  0.0021 |  244 (68.7%) |  0.0014 |  17 (65.4%) |  0.0014
TABLE X: Algorithm performance broken down by rated power (NEMA motors) Case | 0.37 - 3.6kW | 4 - 15kW | 18.5 - 75kW | 90 - 185kW | 200 - 630kW |  $>$630kW
---|---|---|---|---|---|---
Conv-ergence |  Average Error2 |  Conv-ergence |  Average Error2 |  Conv-ergence |  Average Error2 |  Conv-ergence |  Average Error2 |  Conv-ergence |  Average Error2 |  Conv-ergence |  Average Error2
Newton-Raphson ($k_{x}=1$, $k_{r}=0.5$) |  9 (1.43%) |  0.3225 |  100 (16.7%) |  0.3254 |  415 (56.0%) |  0.1455 |  162 (57.0%) |  0.1052 |  65 (52.9%) |  0.0469 |  0 (0.0%) |  0.2338
Levenberg-Marquardt ($k_{x}=1$, $k_{r}=0.5$) |  83 (13.2%) |  0.2788 |  158 (26.4%) |  0.2833 |  395 (53.3%) |  0.3086 |  102 (35.9%) |  0.3288 |  31 (25.2%) |  0.1653 |  1 (50.0%) |  0.0765
Damped NR ($k_{x}=1$, $k_{r}=0.5$) |  58 (9.21%) |  0.0820 |  110 (18.4%) |  0.0751 |  297 (40.1%) |  0.0722 |  74 (26.1%) |  0.1683 |  28 (22.8%) |  0.0863 |  1 (50.0%) |  0.2711
Genetic Algorithm (Max gens = 30) |  0 (0.0%) |  0.0366 |  0 (0.0%) |  0.0362 |  0 (0.0%) |  0.0449 |  0 (0.0%) |  0.0515 |  0 (0.0%) |  0.0493 |  0 (0.0%) |  0.0926
Genetic Algorithm (Max gens = 50) |  0 (0.0%) |  0.0306 |  0 (0.0%) |  0.0303 |  0 (0.0%) |  0.0374 |  0 (0.0%) |  0.0405 |  0 (0.0%) |  0.0344 |  0 (0.0%) |  0.0872
Hybrid NR-GA |  161 (25.6%) |  0.0317 |  311 (52.0%) |  0.0245 |  512 (69.1%) |  0.0226 |  123 (43.3%) |  0.0294 |  51 (41.5%) |  0.0059 |  1 (50.0%) |  0.0416
Hybrid LM-GA |  237 (37.6%) |  0.519 |  308 (51.5%) |  0.5374 |  413 (55.7%) |  0.4339 |  153 (53.9%) |  0.5693 |  69 (56.1%) |  0.5659 |  1 (50.0%) |  0.0554
Hybrid DNR-GA |  235 (37.3%) |  0.02 |  307 (51.3%) |  0.0181 |  405 (54.7%) |  0.0192 |  151 (53.2%) |  0.0141 |  69 (56.1%) |  0.0125 |  1 (50.0%) |  0.0586

Tables IX and X show the convergence and average squared error rates for the
IEC and NEMA motor data sets respectively, with the data sets subdivided by
rated power. It is observed that the convergence and average error rates are
not uniformly distributed across the full range of motor rated powers.

Of interest is the poor performance of all algorithms on smaller motors,
particularly motors rated below 4kW, where convergence rates are between 0.3%
and 12.3% for IEC motors and between 1.4% and 37.3% for NEMA motors.
Performance begins to improve for all algorithms as the motor size is
increased. For motors $\geq$90kW, the convergence rates of the hybrid
algorithms improve to $>$60% for the IEC motors and $>$50% for the NEMA
motors. It should be noted that dynamic modelling is least likely to be
performed on individual small motors, since they are often aggregated as
lumped loads in power system studies.

##  VII Conclusion

In this paper, a number of algorithms were investigated for the estimation of
induction motor parameters based on manufacturer data. Simulations on a large
data set of IEC and NEMA motors showed that the conventional Newton-Raphson
algorithm performs quite poorly. Hybrid algorithms were introduced as an
alternative to the conventional NR algorithm and computer simulations
suggested that the proposed hybrid algorithms show promise as a parameter
estimation tool, with large improvements in convergence and error rates over
the conventional algorithms.

The key drawback for the hybrid algorithms is their computation time, which
depending on the algorithm settings, can be significantly slower than
conventional descent or genetic algorithms. In any case, it can be argued that
motor parameter estimation for the purpose of power system studies is not
particularly time critical and a slow computation time can be tolerated in
return for better algorithm performance.

## References

  * [1] P. Vas, _Electric Machines and Drives. A Space-Vector Theory Approach_ , Clarendon Press, London, 1993
  * [2] B. K. Johnson and J. R. Willis, _Tailoring induction motor analytical models to fit known performance characteristics and satisfy particular study needs_ ,  IEEE Transactions on Power Systems, vol. 6, no. 3, 1991.
  * [3] G. Rogers and D. Shirmohammadi, _Induction machine modelling for electromagnetic transient program_ ,  IEEE Transactions on Energy Conversion, vol. EC-2, no. 4, 1987.
  * [4] S. S. Waters and R. D. Willoughby, _Modeling induction motors for system studies_ ,  IEEE Transactions on Industry Applications, vol. 1A-19, no. 5, 1983\.
  * [5] D. Lindenmeyer, H. W. Dommel, A. Moshref, and P. Kundur, _An induction motor parameter estimation method_ ,  Electrical Power and Energy Systems, vol. 23, pp. 251–262, 2001\.
  * [6] J. Pedra, _Estimation of induction motor double-cage model parameters from manufacturer data_ ,  IEEE Transactions on Energy Conversion, vol. 19, no. 2, 2004.
  * [7] F. Corcoles, J. Pedra, M. Salichs, and L. Sainz, _Analysis of the induction machine parameter identification_ ,  IEEE Transactions on Energy Conversion, vol. 17, no. 2, 2002.
  * [8] J. Pedra, _On the determination of induction motor parameters from manufacturer data for electromagnetic transient programs_ ,  IEEE Transactions on Power Systems, vol. 23, no. 4, 2008.
  * [9] K. Levenberg, _A method for the solution of certain non-linear problems in least squares_ ,  The Quarterly of Applied Mathematics, vol. 2, 1944.
  * [10] D.W. Marquardt, _An algorithm for least-squares estimation of non-linear parameters_ ,  Journal of the Society for Industrial and Applied Mathematics, vol. 11, no. 2, 1963.
  * [11] P. Nangsue, P. Pillay, and S. E. Conry, _Evolutionary algorithms for induction motor parameter determination_ ,  IEEE Transactions on Energy Conversion, vol. 14, no. 3, 1999.
  * [12] H.H. Weatherford and C.W. Brice, _Estimation of induction motor parameters by a genetic algorithm_ ,  in Conference Record of the 2003 Annual Pulp and Paper Industry Technical Conference, 2003., 2003, pp. 21–28.
  * [13] R Nolan, P Pillay, and T Haque, _Application of genetic algorithms to motor parameter determination_ ,  in Proceedings of 1994 IEEE Industry Applications Society Annual Meeting, 1994, pp. 47–54.
  * [14] R. Bhuvaneswari and S. Subramanian, _Optimization of three-phase induction motor design using simulated annealing algorithm_ ,  Electric Power Components and Systems, vol. 33, no. 9, pp. 947–956, 2005.
  * [15] V.P. Sakthivel, R. Bhuvaneswari, and S. Subramanian, _Multi-objective parameter estimation of induction motor using particle swarm optimization_ ,  Engineering Applications of Artificial Intelligence, vol. 23, no. 3, pp. 302–312, 2010.
  * [16] V.P. Sakthivel, R. Bhuvaneswari, and S Subramanian, _An Improved Particle Swarm Optimization for Induction Motor Parameter Determination_ ,  International Journal of Computer Applications, vol. 1, no. 2, pp. 71–76, 2010.
  * [17] V.P. Sakthivel, R. Bhuvaneswari, and S. Subramanian, _Artificial immune system for parameter estimation of induction motor_ ,  Expert Systems with Applications, vol. 37, no. 8, pp. 6109–6115, 2010.
  * [18] V. P. Sakthivel, R. Bhuvaneswari, and S. Subramanian, _Bacterial Foraging Technique Based Parameter Estimation of Induction Motor from Manufacturer Data_ ,  Electric Power Components and Systems, vol. 38, no. 6, pp. 657–674, 2010.
  * [19] Zhenfeng Chen, Yanru Zhong, and Jie Li, _Parameter identification of induction motors using ant colony optimization_ ,  in Evolutionary Computation, 2008. CEC 2008. (IEEE World Congress on Computational Intelligence). IEEE Congress on, 2008, pp. 1611–1616.
  * [20] J.R. Marques, I.F. Machado, and J.R. Cardoso, _Induction motor parameter determination using the harmony search algorithm to power, torque and speed estimation_ ,  in IECON 2012 - 38th Annual Conference on IEEE Industrial Electronics Society, 2012, pp. 1835–1840.
  * [21] EuroDEEM, _Eurodeem: The european database of efficient electric motors_ , http://sunbird.jrc.it/energyefficiency/eurodeem/, 2007.

Julius Susanto (M’11) received the B.Eng and M.Sc degrees in electrical
engineering from Curtin University, Perth, Australia, in 2001 and 2013
respectively. Since 2002 he has worked in industry for WorleyParsons Ltd,
DIgSILENT GmbH, TransGrid and Synergy Engineering focusing on power system
studies. He is a Chartered Professional Engineer in Australia.
---
Syed Islam (M’83, SM’93) received the B.Sc. from Bangladesh University of
Engineering and Technology, Bangladesh, M.Sc. and PhD degrees from King Fahd
University of Petroleum and Minerals, Saudi Arabia, all in electrical power
engineering in 1979, 1983 and 1988 respectively. He is currently the Chair
Professor in Electrical Power Engineering at Curtin University, Perth,
Australia. He received the IEEE T Burke Haye s Faculty Recognition award in
2000. His research interests are in Condition Monitoring of Transformers, Wind
Energy Conversion and Power Systems. He is regular reviewer for the IEEE
Transactions on Energy Conversion, Power Systems and Power Delivery. Prof.
Islam is an editor of the IEEE Transaction on Sustainable Energy.
---


# Identifying Single-Input Linear System Dynamics from Reachable Sets

Taha Shafa  Roy Dong  Melkior Ornik TS and MO are with the Department of
Aerospace Engineering at the University of Illinois Urbana-Champaign. RD is
with the Department of Electrical and Computer Engineering at the University
of Illinois Urbana-Champaign. {tahaas2, roydong<EMAIL_ADDRESS>
###### Abstract

This paper is concerned with identifying linear system dynamics without the
knowledge of individual system trajectories, but from the knowledge of the
system’s reachable sets observed at different times. Motivated by a scenario
where the reachable sets are known from partially transparent manufacturer
specifications or observations of the collective behavior of adversarial
agents, we aim to utilize such sets to determine the unknown system’s
dynamics. This paper has two contributions. Firstly, we show that the sequence
of the system’s reachable sets can be used to uniquely determine the system’s
dynamics for asymmetric input sets under some generic assumptions, regardless
of the system’s dimensions. We also prove the same property holds up to a sign
change for two-dimensional systems where the input set is symmetric around
zero. Secondly, we present an algorithm to determine these dynamics. We apply
and verify the developed theory and algorithms on an unknown band-pass filter
circuit solely provided the unknown system’s reachable sets over a finite
observation period.

##  I Introduction

This paper aims to determine whether it is possible to use a control system’s
reachable sets obtained at different time instances to calculate the system’s
dynamics. In certain instances, we may be able to determine an approximation
of a system’s reachable sets over a finite observation period. The purpose of
this paper is to show that such information can be utilized to arrive at a
dynamic model for an unknown system. Practical applications may include system
identification of high-density drone and missile swarms [1, 2] where the
reachable set can be found by observing multiple agents collectively, but
without the capability of distinguishing them. Other applications include
predicting macro-level population behaviors, e.g., determining how crowd
behavior changes under social or economic events like the introduction of a
new population or changes in the stock market [3]. We may also be able to
model internal body functions on the cellular level [4, 5], namely
understanding how cells change their identity and behavior in living systems.

We must first show that model identification using reachable sets will
uniquely determine an unknown system’s true dynamics. After uniqueness is
proven, we develop a method to identify a linear model of an unknown system’s
behavior using its reachable sets. Previous research in system identification
presents the most closely related contributions to the method presented in
this paper. However, previous work on system identification classically relies
on frequency response techniques induced by randomized actuator inputs [6, 7].
More sophisticated system identification techniques involve neural networks
[8]. Single-layer and multi-layer neural networks have also been applied with
the use of parameter estimation algorithms using a single hidden layer [9] and
$H_{\infty}$ control-induced excitations for robust identification of system
nonlinearities [10]. More recent work involves using recurrent neural networks
[11, 12] with Long Short-Term Memory Units (LSTM) and fractional order neural
networks (FONN) [13, 14] to identify and control dynamic systems. These
methods, however, cannot be used unless one has access to a system’s actuators
or individual trajectories. The significant difference of our novel method is
that it does not require control of any actuators to model an unknown system
nor observations of individual trajectories.

On a high level, the problem in this paper involves identifying the behaviors
or capabilities of an observed system under limited information. While there
exist other methods for adversarial behavior recognition, those works are
focused on determining adversarial agent goals by matching actions of an agent
against a plan library [15, 16, 17]. More recent work [18, 19] proposes using
evolving fuzzy systems and artificial intelligence to adaptively predict agent
behavior. In contrast, our method is starkly different since it is not
primarily concerned with predicting adversarial behavior, but determining all
possible actions of an adversary within a time horizon. Thus, instead of using
a library of finite predetermined adversarial actions, our method uses
reachable sets to produce a dynamic model of an unknown system.

The outline of this paper is as follows: in Section II, we discuss the problem
statement, namely posing the question of whether linear dynamics can be
uniquely recovered given an unknown system’s sequence of reachable sets and
how to recover said dynamics. In Section III, we address the question of
whether the system dynamics are uniquely determined by the system’s reachable
sets. We show that under generic assumptions, the system dynamics are indeed
unique under asymmetric input sets. For unknown systems with input sets
symmetric around zero, uniqueness modulo a sign has been proved in the two-
dimensional case; we conjecture the same holds for higher dimensions. In
Section IV, we propose a procedure using knowledge of the reachable sets to
calculate the system dynamics. In Section V, we illustrate by example how to
implement this procedure to identify the models of an unknown band-pass filter
circuit and an additional dynamical system with a symmetric input set.

###  I-A Notation

We denote the set of all $n\times m$ real and complex matrices by
$\mathbb{R}^{n\times m}$ and $\mathbb{C}^{n\times m}$ respectively; for
$M\in\mathbb{R}^{n\times m}$, we let $M^{T}\in\mathbb{R}^{m\times n}$ denote
its transpose. Vectors $e_{1},\ldots,e_{n}$ will denote the canonical basis
vectors in $\mathbb{R}^{n}$. We let $\mathbb{N}$ denote the set of all natural
numbers, $\mathbb{Z}_{\geq 0}$ denote the set of non-negative integers, and
$GL(n)$ denote the set of invertible square matrices of dimension
$n\in\mathbb{N}$. Let $\mathcal{S}$ be a set of points in $\mathbb{R}^{n}$.
Then $\mathrm{Conv}(\mathcal{S})$ denotes the convex hull of $\mathcal{S}$.
Notation $B\mathcal{X}$ where $B\in\mathbb{R}^{n\times m}$ and
$\mathcal{X}\subset\mathbb{R}^{m}$ denotes the set
$B\mathcal{X}=\\{Bx~{}|~{}x\in\mathcal{X}\\}$. Given two sets
$\mathcal{A},\,\mathcal{B}\in\mathbb{R}^{n}$, we denote
$\mathcal{A}\oplus\mathcal{B}=\\{a+b~{}|~{}a\in\mathcal{A},\,b\in\mathcal{B}\\}$
as their Minkowski sum. Similarly,
$\mathcal{A}\ominus\mathcal{B}=\\{c\in\mathbb{R}^{n}~{}|~{}c\oplus\mathcal{B}\subseteq\mathcal{A}\\}$
denotes the Minkowski difference. We also define
$\mathcal{A}+b=\\{a+b~{}|~{}a\in\mathbb{R}^{n}\\}$ as the translation of
$\mathcal{A}$ by $b\in\mathbb{R}^{n}$.

##  II Problem Statement

We consider the discrete-time, single-input linear system

$x[i+1]=Ax[i]+bu[i],\quad x[0]=0,$ (1)

where all $i\in\mathbb{Z}_{\geq 0}$, $x\in\mathbb{R}^{n}$,
$A\in\mathbb{R}^{n\times n}$, $b\in\mathbb{R}^{n}$ and
$u~{}\in~{}\mathcal{U}~{}\subset~{}\mathbb{R}$ where
$\mathcal{U}=[\underline{u},\overline{u}]$ such that
$\underline{u}\neq\overline{u}$. We assume $b\neq 0$ since the system’s
reachable sets are trivial otherwise. We also assume $x[0]=0$; by a shift in
coordinates, the case of $x[0]\neq 0$ is equivalent to that of an affine
system $x[i+1]=Ax[i]+bu[i]+c$ with initial state at the origin. Solving the
problem in this setting can likely be approached by reproducing similar
calculations in subsequent sections, but we leave such an effort for future
work.

Our goal is to establish whether the dynamics of (1), i.e., matrices $A$ and
$b$, can be determined using the system’s reachable sets. We now formally
define said reachable sets.

###### Definition 1

For $i\in\mathbb{Z}_{\geq 0}$, the (forward) reachable set of system (1) at
time $i$ is

$\mathcal{R}(i,x[0])=\\{\phi_{u}(i;x[0])~{}|~{}u:\mathbb{Z}_{\geq
0}\to\mathcal{U}\\},$

where $\phi_{u}(\cdot;x[0])$ denotes the controlled trajectory of system (1)
with control signal $u$.

We present the problem of whether the system dynamics are uniquely determined
by the system’s reachable sets.

###### Problem 1

Given a sequence of sets $\\{\mathcal{R}(i,0)\\}_{i\in\mathbb{N}}$ which is
generated by (1) for some $(A,b)$, determine whether $(A,b)$ can be uniquely
recovered from $\\{\mathcal{R}(i,0)\\}_{i\in\mathbb{N}}$.

Notice that we explicitly assume the knowledge of all reachable sets at all
times. Such an assumption might not always be realistic. We will show that we
often need only the first $n+1$ reachable sets to uniquely recover the
dynamics. We leave the more general case – where only reachable sets at
different time steps are available – for future work.

The first step to solving Problem 1 is to derive a simple relationship between
the system matrices and $\mathcal{R}(i,0)$. Given system (1), we naturally
utilize Minkowski sums and the Minkowski difference [20] to produce such a
relationship for all $i\in\mathbb{N}$.

###### Theorem 1

Let $\mathcal{R}(i,0)$ be the reachable set at time $i$ of (1). Then

$A^{i-1}b\mathcal{U}=\mathcal{R}(i,0)\ominus\mathcal{R}(i-1,0).$ (2)

###### Proof:

By (1) it is clear that $\mathcal{R}(1,0)=b\mathcal{U}$. Since

$x[i]=A^{i}x[0]+A^{i-1}bu[0]+\ldots+bu[i-1],$

clearly

$\mathcal{R}(i,0)=A^{i-1}b\mathcal{U}\oplus\ldots\oplus b\mathcal{U}$

and hence

$\mathcal{R}(i,0)=A^{i-1}b\mathcal{U}\oplus\mathcal{R}(i-1,0).$

We recall that the Minkowski sum of two convex sets is also convex [21]. Since
all sets $A^{i-1}b\mathcal{U}$ are convex by the definition of $\mathcal{U}$,
all sets $\mathcal{R}(i,0)$ are convex by induction. Hence, the appropriate
Minkowski difference [22] can be calculated to arrive at (2). ∎

Theorem 1 implies that we can obtain
$\\{A^{i-1}b\mathcal{U}\\}_{i\in\mathbb{N}}$ using the reachable sets
$\mathcal{R}(i,0)$. We will prove that when $\mathcal{U}\neq[-c,c]$ with
$c\in\mathbb{R}$, matrices $A$ and $b$ are indeed generically uniquely defined
from $\\{A^{i-1}b\mathcal{U}\\}_{i\in\mathbb{N}}$, that is, uniquely defined
under the assumptions formally written in Theorem 2 shown to be generic in a
topological sense in Lemma 1. When $\mathcal{U}=[-c,c]$ for some
$c\in\mathbb{R}$, we can show that $(A,b)$ are not uniquely defined, but
conjecture that they are unique up to a change in sign. We prove that this
property holds for $n=2$. We shall refer to solutions for cases with such a
set $\mathcal{U}$ as $\pm$-unique, which is explicitly defined in the next
section.

Following Problem 1, which seeks to determine whether system dynamics are
uniquely defined from reachable sets, we present the second problem, which
aims to explicitly determine such dynamics.

###### Problem 2

Develop a method to recover at least one pair $(A,b)$ which generates
$\\{\mathcal{R}(i,0)\\}_{i\in\mathbb{N}}$.

Based on methods in [20] for calculating Minkowski differences, we can
calculate $\\{A^{i-1}b\mathcal{U}\\}_{i\in\mathbb{N}}$. We show in Section IV
that the results of these Minkowski differences and knowledge of $\mathcal{U}$
are sufficient for calculating $\\{A^{i-1}b\\}_{i\in\mathbb{N}}$, which in
turn can be utilized to calculate the matrix pair $(A,b)$ for controllable
systems. We first tackle Problem 1.

##  III Uniqueness of the Derived System Model

We wish to determine when any pair $(A,b)$ uniquely defines the dynamics of
(1). It can be easily shown that the answer is generally negative. Consider an
unknown system (1) where

$A=\begin{bmatrix}0&0\\\ 0&1\end{bmatrix},\quad b=\begin{bmatrix}0\\\
1\end{bmatrix}$ (3)

and $\mathcal{U}=[0,1]$. By equation (2) of Theorem 1, we see that if
$A^{\prime}=I$, then the reachable sets of (1) with matrix pairs $(A,b)$ and
$(A^{\prime},b)$ are equivalent. Thus, we begin by determining sufficient
conditions which guarantee whether $(A,b)$ can be uniquely recovered as stated
in Problem 1. We will show uniqueness under several technical assumptions;
Lemma 1 shows said assumptions are generic in a topological sense.

###### Lemma 1

Let $\mathcal{N}\subset\mathbb{R}^{n\times n}$ be the set of all matrices such
that if $A\in\mathcal{N}$, then $A^{2}$ has distinct eigenvalues. Let
$b\in\mathbb{R}^{n}\backslash\\{0\\}$ and $\mathcal{O}\in\mathbb{R}^{n\times
n}$ be the set of all matrices such that, if $A\in\mathcal{O}$ and
$\eta\in\mathbb{C}^{n}$ is any left eigenvector of $A$, $b^{T}\eta\neq 0$.
Then, $GL(n)\cap\mathcal{N}\cap\mathcal{O}$ is an open and dense set.

###### Proof:

It is a well known result that the set of all matrices with distinct
eigenvalues and the set $GL(n)$ are both open and dense [23]. Clearly,
openness of the former set implies $\mathcal{N}$ is open. To show
$\mathcal{N}$ is also dense, we would follow similar steps as part of the
proof to show $GL(n)\cap\mathcal{N}\cap\mathcal{O}$ is dense. For
succinctness, we prove $GL(n)\cap\mathcal{N}\cap\mathcal{O}$ is open and dense
and leave the proof that $\mathcal{N}$ is dense to the reader.

Openness of $GL(n)\cap\mathcal{N}\cap\mathcal{O}$ can be trivially concluded
by the continuity of eigenvectors [24], meaning if we consider a matrix $A(t)$
whose elements are a continuous function of $t$, any eigenvectors $v_{i}(t)$
and left eigenvectors $\eta_{i}(t)$ of norm $1$ of $A(t)$ are continuous
function of $t$.

We now prove denseness. In other words, we will show that for any arbitrary
matrix $A$ and any $\epsilon>0$, there exists a matrix $A^{\prime\prime}\in
GL(n)\cap\mathcal{N}\cap\mathcal{O}$ such that
$\|A~{}-~{}A^{\prime\prime}\|~{}<~{}\epsilon$. Let $\eta_{i}$ be the left
eigenvectors of $A$ so that $\eta_{i}^{T}A~{}=~{}\lambda_{i}\eta_{i}^{T}$. By
the denseness of $GL(n)$, for any $\delta>0$ we can find vectors
$\eta_{i}^{\prime T}$ such that $\|\eta_{i}^{T}~{}-~{}\eta_{i}^{\prime
T}\|~{}<~{}\delta$ for all $i$ and $\mathrm{det}([\eta_{1}^{\prime}\textrm{
}\eta_{2}^{\prime}\textrm{ }\cdots\textrm{ }\eta_{n}^{\prime}])~{}\neq~{}0$.
By the continuity of determinants and because $b\neq 0$, we can slightly
perturb one element of $\eta_{i}^{\prime T}$ to obtain $\eta_{i}^{\prime\prime
T}$ such that $\mathrm{det}([\eta_{1}^{\prime\prime}\textrm{
}\eta_{2}^{\prime\prime}\textrm{ }\cdots\textrm{
}\eta_{n}^{\prime\prime}])~{}\neq~{}0,\,\|\eta_{i}^{T}-\eta_{i}^{\prime\prime
T}\|<\delta$ for all $i$, and $b^{T}\eta_{i}\neq 0$. We now let
$\eta_{i}^{\prime\prime}$ form a basis in $\mathbb{C}^{n}$, and define a
matrix $A^{\prime}$ such that $\eta_{i}^{\prime\prime
T}A^{\prime}=\lambda_{i}\eta_{i}^{\prime\prime T}$ and
$A^{\prime}\in\mathcal{O}$. If the perturbations above are performed in a way
that ensure that perturbations of real eigenvectors remain real, and
perturbations of complex conjugate vectors remain complex conjugates, matrix
$A^{\prime}$ is real [25].

Since $\eta_{i}^{\prime\prime}$ form a basis in $\mathbb{C}^{n}$, we can
represent any vector $x\in\mathbb{R}^{n}$ as
$x=\sum_{i=1}^{n}\beta_{i}(x)\eta_{i}^{\prime\prime T}$ where
$\beta_{i}(x)\in\mathbb{R}$. We can compute $\beta_{i}(x)$ as a continuous
function of $x$. Recall that $\|A\|=\max_{\|x\|=1}\|Ax\|$. We consider $x$
such that $\|x\|=1$. Then $\beta_{i}(x)$ is a continuous function on a compact
space and thus has a maximum. Let
$\alpha_{i}=\max\\{|\beta_{i}(x)|~{}|~{}\|x\|=1\\}$. Note that
$x^{T}A^{\prime}=\sum_{i=1}^{n}\lambda_{i}\beta_{i}(x)\eta_{i}^{\prime\prime
T}$. It follows that

$\|x^{T}A-x^{T}A^{\prime}\|\leq\sum_{i=1}^{n}\|(\beta_{i}(x)\eta_{i}^{\prime\prime
T})A-\beta_{i}(x)\lambda_{i}\eta_{i}^{\prime\prime T}\|$
$=\sum_{i=1}^{n}\|\beta_{i}(x)(\eta_{i}^{T}A-(\eta_{i}^{T}-\eta_{i}^{\prime\prime
T})A)-\beta_{i}(x)\lambda_{i}\eta_{i}^{\prime\prime T}\|$
$=\sum_{i=1}^{n}\|\beta_{i}(x)((\eta_{i}^{\prime\prime
T}-\eta_{i}^{T})A+\lambda_{i}(\eta_{i}^{T}-\eta_{i}^{\prime\prime T}))\|$
$<\sum_{i=1}^{n}(\|\alpha_{i}A\|+\|\alpha_{i}\lambda_{i}\|)\delta$

and so if we set
$\delta=\epsilon/(2\sum_{i=1}^{n}\|\alpha_{i}A\|+\|\alpha_{i}\lambda_{i}\|)$,
then $\|x^{T}A-x^{T}A^{\prime}\|<\epsilon/2$.

Given $\lambda_{1},\ldots,\lambda_{n}$, for any $\rho>0$ we can obviously find
a set $\\{\lambda_{1}^{\prime},\ldots,\lambda_{n}^{\prime}\\}$ such that
$|\lambda_{i}-\lambda_{i}^{\prime}|<\rho$ for all $i$,
$\lambda_{i}^{\prime}\neq 0$ for all $i$, and
$\lambda_{i}=\overline{\lambda_{j}}$ implies
$\lambda_{i}^{\prime}=\overline{\lambda_{j}^{\prime}}$. Now, define
$A^{\prime\prime}$ by $\eta_{i}^{\prime\prime
T}A^{\prime\prime}=\lambda_{i}^{\prime}\eta_{i}^{\prime\prime T}$ and
$A^{\prime\prime}\in\mathcal{M}\cap\mathcal{N}\cap\mathcal{O}$. As before, if
the perturbation of eigenvalues is performed in such a way that real
eigenvalues remain real and complex conjugates remain conjugate,
$A^{\prime\prime}$ is real. It follows that

$\|x^{T}A^{\prime}-x^{T}A^{\prime\prime}\|\leq\sum_{i=1}^{n}\|\beta_{i}(x)\lambda_{i}\eta_{i}^{\prime\prime
T}-\beta_{i}(x)\lambda_{i}^{\prime}\eta_{i}^{\prime\prime T}\|$
$=\sum_{i=1}^{n}\|\beta_{i}(x)\eta_{i}^{\prime\prime
T}(\lambda_{i}-\lambda_{i}^{\prime})\|<\sum_{i=1}^{n}\|\alpha_{i}\eta_{i}^{\prime\prime
T}\|\rho.$

If we set $\rho=\epsilon/(2\sum_{i=1}^{n}\|\alpha_{i}\eta_{i}^{\prime\prime
T}\|)$, then $\|x^{T}A^{\prime}-x^{T}A^{\prime\prime}\|<\epsilon/2$. Finally
we have
$\|x^{T}Ax-x^{T}A^{\prime\prime}\|=\|x^{T}A-x^{T}A^{\prime\prime}+x^{T}A^{\prime}-x^{T}A^{\prime}\|\leq\|x^{T}A-x^{T}A^{\prime}\|+\|x^{T}A^{\prime}-x^{T}A^{\prime\prime}\|<\epsilon/2+\epsilon/2=\epsilon$.
Since this inequality holds for all $x$ such that $\|x\|=1$, indeed
$\|A-A^{\prime\prime}\|<\epsilon$, and the claim is proven. ∎

We emphasize that many well-known linear controllable systems, such as the
discrete double integrator, RLC circuit, and linearized pendulum [26], contain
$A$ matrices which satisfy the conditions of Lemma 1. Also, these generic
assumptions are not necessary, but sufficient to guarantee uniqueness. For
example, a row perturbation of $A$ in (3) clearly does not satisfy the generic
assumptions in Lemma 1, but the reachable sets of (1) with this new matrix can
be used to uniquely generate the dynamics, which implies this method can be
applied to a larger set of systems. Finding such non-generic assumptions which
guarantee uniqueness is a highly involved problem and remains for future work.
In the proof below, we will use the assumptions in Lemma 1 to prove that the
dynamics derived from reachable sets are generically unique, at least for an
asymmetric input set.

###### Theorem 2

Let $\mathcal{U}=[c,d]$, where $c\neq\pm d$. Let
$\eta_{i}~{}\in~{}\mathbb{C}^{n}$ for $i\in\\{1,\ldots,n\\}$ be the left
eigenvectors of $A$. Let the sequence
$\\{\mathcal{R}(j,0)\\}_{j\in\mathbb{N}}$ be generated by system (1) for
system matrices $(A,b)$ and $(A^{\prime},b^{\prime})$, where
$A,\,A^{\prime}\in GL(n)$, $A$ and $A^{\prime}$ have $n$ distinct eigenvalues,
and $b^{T}\eta_{i}\neq 0$ for all $i$. Then, $(A,b)=(A^{\prime},b^{\prime})$.

###### Proof:

If $(A,b)$ and $(A^{\prime},b^{\prime})$ for system (1) produce an identical
sequence $\\{\mathcal{R}(j,0)\\}_{j\in\mathbb{N}}$, then
$\mathcal{R}(1,0)=b\mathcal{U}=b^{\prime}\mathcal{U}$, i.e., there are two
options: (i) $bc=b^{\prime}c$ and $bd=b^{\prime}d$ or (ii) $bc=b^{\prime}d$
and $bd=b^{\prime}c$. If the latter option is true, then
$bcd=b^{\prime}d^{2}=b^{\prime}c^{2}$, so $b^{\prime}=0$. In that case, $b=0$,
so $b=b^{\prime}$. If the former option is true, because at least one of $c$
or $d$ is non-zero, again $b=b^{\prime}$.

Let us now perform a coordinate transformation $z=Mx$, where $M$ is chosen so
that $Mb=e_{1}$. Such an $M$ exists since $b\neq 0$. Then,
$\dot{z}=MAx+Mbu=MAM^{-1}z+e_{1}u$. If we define $\hat{A}=MAM^{-1}$, by our
assumptions $\hat{A}$ is invertible and $\hat{A}$ has distinct eigenvalues.
Additionally, it is trivially verified that left eigenvectors of $\hat{A}$ are
$(M^{-1})^{T}\eta_{i}$. Since $M^{-1}e_{1}=b$, the assumption
$b^{T}\eta_{i}\neq 0$ is equivalent to the first element of the left
eigenvectors of $\hat{A}$ being non-zero. To simplify the notation, by a
standard abuse we now let $(A,e_{1})$, $(A^{\prime},e_{1})$ represent the
system matrices after performing the above transformation. By the above
discussion, we are then assuming that $A$ and $A^{\prime}$ are invertible,
have distinct eigenvalues, and that $\eta_{i1}\neq 0$ for all $i$.

Noting that the two systems produce the same reachable sets, by (2) it follows
that $A^{k}e_{1}=A^{\prime k}e_{1}$ for all $k\in\mathbb{N}$. By the same
logic as in the first paragraph of the proof, we see that since $c\neq-d$,
then $A^{k}ce_{1}=A^{\prime k}ce_{1}$ and $A^{k}de_{1}=A^{\prime k}de_{1}$ is
satisfied for all $k\in\mathbb{N}$, giving us the relation

$A^{k}e_{1}=A^{\prime k}e_{1}\quad\forall\,k\in\mathbb{Z}_{\geq 0}.$ (4)

Equation (4) implies $A^{k-1}A^{\prime}e_{1}=A^{\prime k-1}Ae_{1}$ and
$A^{\prime k-2}Ae_{1}=A^{k-1}e_{1}$ for all $k\geq 2$. We have

$A^{k-1}A^{\prime}e_{1}=A^{\prime k-1}Ae_{1}=A^{\prime}A^{\prime
k-2}Ae_{1}=A^{\prime}A^{k-1}e_{1}.$

Hence, $A^{k-1}A^{\prime}e_{1}=A^{\prime}A^{k-1}e_{1}$; since $A^{\prime}$ is
invertible,

$A^{k}e_{1}=A^{\prime(-1)}A^{k}A^{\prime}e_{1}\quad\forall\,k\in\mathbb{Z}_{\geq
0}.$ (5)

Let $v_{i}$ denote the right eigenvectors of $A$ and
$v_{i}^{\prime},\,\eta_{i}^{\prime}$ denote the right and left eigenvectors of
$A^{\prime(-1)}AA^{\prime}$ respectively. Since $A$ and
$A^{\prime(-1)}AA^{\prime}$ are similar matrices, their eigenvalues are equal
[25]. Let $A=VDV^{-1}$ and
$A^{\prime(-1)}AA^{\prime}=V^{\prime}DV^{\prime(-1)}$ where the rows of
$V^{-1}$ and $V^{\prime-1}$ are $\eta_{i}^{T}$ and $\eta_{i}^{\prime T}$
respectively and the columns of $V$ and $V^{\prime}$ are $v_{i}$ and
$v_{i}^{\prime}$ respectively. By our assumptions, $\eta_{i1}\neq 0$, so we
can now scale the $\eta$’s so that $\eta_{i1}=1$. We then redefine $v_{i}$ to
be the newly scaled right eigenvectors such that $\eta_{i1}=1$. Next, we write
(5) in tensor notation [25] and get

$\sum_{i}\lambda_{i}^{k}v_{i}=\sum_{i}\lambda_{i}^{k}v_{i}^{\prime}\eta_{i1}^{\prime
T}\quad\forall\,k\in\mathbb{Z}_{\geq 0}$

which implies

$\sum_{i}\lambda_{i}^{k}(v_{i}-v_{i}^{\prime}\eta_{i1}^{\prime
T})=0\quad\forall\,k\in\mathbb{Z}_{\geq 0}.$ (6)

Taking $k\in\\{0,\ldots,n-1\\}$ we have a series of $n$ equations. For the
$j$-th element of any $v_{i}$ and $v_{i}^{\prime}$, we have

$\Lambda S_{j}=\begin{bmatrix}1&\ldots&1\\\ \lambda_{1}&\ldots&\lambda_{n}\\\
\vdots&\vdots&\vdots\\\
\lambda_{1}^{n-1}&\ldots&\lambda_{n}^{n-1}\end{bmatrix}\begin{bmatrix}v_{1j}-v_{1j}^{\prime}\eta_{11}^{\prime
T}\\\ v_{2j}-v_{2j}^{\prime}\eta_{21}^{\prime T}\\\ \vdots\\\
v_{nj}-v_{nj}^{\prime}\eta_{n1}^{\prime T}\end{bmatrix}=\begin{bmatrix}0\\\
0\\\ \vdots\\\ 0\end{bmatrix}$

for any $j\in\\{1,\ldots,n\\}$. Notice that $\Lambda\in\mathbb{C}^{n\times n}$
is the square Vandermonde matrix [27]. Recall that the Vandermonde matrix is
invertible if elements $\lambda_{i}$ are distinct for all $i$, which holds by
assumption. If $\eta_{i1}^{\prime}=0$ for any $i$, then $v_{i}=0$, which
contradicts the assumption that $A$ is diagonalizable. Consequently,
$\eta_{i1}^{\prime}\neq 0$ for all $i$, so similar to the previous step, we
can scale $v_{i}^{\prime}$ and $\eta_{i1}^{\prime}$ such that
$\eta_{i1}^{\prime}=1$ for all $i$. It follows that $v_{ij}=v_{ij}^{\prime}$
for all $i,j$ since $\Lambda$ is invertible. Therefore,
$A=A^{\prime(-1)}AA^{\prime}$.

Recall that we assumed that all eigenvalues of $A$ are distinct. Thus, since
$A$ and $A^{\prime}$ commute, we can conclude that $A$ and $A^{\prime}$ have
the same eigenvectors [28]. Recall that $A$ and $A^{\prime}$ are both
diagonalizable. If we take the eigenvalue expansion of $A$ and $A^{\prime}$
and multiply both on the left by $V^{-1}$, then equation (4) implies

$D^{k}\begin{bmatrix}\eta_{11}\\\ \eta_{21}\\\ \vdots\\\
\eta_{n1}\end{bmatrix}=D^{\prime k}\begin{bmatrix}\eta_{11}\\\ \eta_{21}\\\
\vdots\\\ \eta_{n1}\end{bmatrix}\quad\forall~{}k\in\mathbb{N},$

where $D^{\prime}$ is the diagonal matrix with eigenvalues of $A^{\prime}$ on
the diagonal. Subtracting the right hand side from both sides reveals that (4)
implies

$(\lambda_{i}^{k}-\lambda_{i}^{\prime
k})\eta_{i1}=0\quad\forall~{}k\in\mathbb{N}.$

By assumption, $\eta_{i1}\neq 0$ for all $i$, so

$\lambda_{i}^{k}=\lambda_{i}^{\prime k}\quad\forall~{}k\in\mathbb{N}.$

Therefore, both $A$ and $A^{\prime}$ have the same eigenvectors and
eigenvalues, hence $A=A^{\prime}$. ∎

Theorem 2 proves that given reachable sets of generic system (1), the pair
$(A,b)$, i.e., the system dynamics, are uniquely defined when the set of
control inputs is not symmetric around $0$. We now want to address the
degenerate case where $\mathcal{U}=[-c,c]$. It can be easily seen that in such
a case, system (1) with $(A,b)$ and $(-A,-b)$ will produce the same reachable
sets. To discuss a relaxed notion of system uniqueness, we provide a formal
definition of $\pm$-uniqueness.

###### Definition 2

The system dynamics $(A,b)$ of (1) are $\pm$-unique if $(A,b)$ and $-(A,b)$
generate the same reachable sets, but there do not exist other pairs
$(A^{\prime},b^{\prime})$ which generate the same reachable sets.

We conjecture that in the case when $\mathcal{U}$ is symmetric around $0$ – a
scenario common in many controls applications [29] – the dynamics are
$\pm$-unique.

###### Conjecture 1

Let $\mathcal{U}=[-c,c]$. Let the sequence
$\\{\mathcal{R}(i,0)\\}_{i\in\mathbb{N}}$ be generated by $(A,b)$, where
$A^{2}$ has distinct eigenvalues and $(A,b)$ are known to satisfy the
assumptions of Theorem 2. Then, $(A,b)$ is $\pm$-unique.

Proving the conjecture above requires extensive theoretical developments and
remains for future work. As an illustration, we formally prove the conjecture
to be true in the two-dimensional case.

###### Theorem 3

Let $n=2$. Then, Conjecture 1 is correct.

###### Proof:

Similarly to the proof of Theorem 2, we have two options: $bc=b^{\prime}c$ or
$bc=-b^{\prime}c$. In the former case, we reach the same result as before,
namely $b=b^{\prime}$. In the latter case, we obtain $b=-b^{\prime}$.
Altogether, we get $b=(-1)^{p(0)}b^{\prime}$ where $p(0)\in\\{0,1\\}$.

As in Theorem 2, through a coordinate transformation, we assume without loss
of generality that $b^{\prime}=e_{1}$. Then
$b\mathcal{U}=(-1)^{p(0)}b^{\prime}\mathcal{U}=(-1)^{p(0)}[-c,c]e_{1}$.
Following the same steps as in the beginning of the proof in Theorem 2, with a
standard abuse of notation, we let $A$, $A^{\prime}$ represent the system
dynamics in this new basis where $A$ and $A^{\prime}$ satisfy our assumptions.
Also, we find that if $\mathcal{U}=[-c,c]$, then we arrive at the relation

$A^{k}e_{1}=(-1)^{p(k)}A^{\prime k}e_{1}\quad\forall\,k\in\mathbb{Z}_{\geq
0}.$ (7)

When $k=2$, we see that regardless of $p(1)$,
$AA^{\prime}e_{1}=(-1)^{p(2)}A^{\prime}Ae_{1}$. Using this fact along with
equation (7) implies $A^{k-1}A^{\prime}e_{1}=(-1)^{p(k)}A^{\prime k-1}Ae_{1}$
for all $k\geq 1$ and $A^{k-2}A^{\prime}e_{1}=(-1)^{p(k-1)}A^{\prime
k-2}Ae_{1}$ for all $k\geq 2$. We have

$(-1)^{p(k)}A^{\prime k-1}Ae_{1}=(-1)^{p(k)}A^{\prime}A^{\prime k-2}Ae_{1}$
$=(-1)^{p(k)}(-1)^{p(k-1)}(-1)^{p(1)}A^{\prime}A^{k-1}e_{1}.$

Hence,
$A^{k-1}A^{\prime}e_{1}=(-1)^{p(k)}(-1)^{p(k-1)}(-1)^{p(1)}A^{\prime}A^{k-1}e_{1}$;
since $A^{\prime}$ is invertible,

$A^{k-1}e_{1}=\frac{A^{\prime(-1)}A^{k-1}A^{\prime}e_{1}}{(-1)^{p(k)}(-1)^{p(k-1)}(-1)^{p(1)}}.$

We define $q(k)\in\\{0,1\\}$ by
$(-1)^{q(k)}~{}=~{}((-1)^{p(k)}(-1)^{p(k-1)}(-1)^{p(1)})^{-1}=(-1)^{p(k)}(-1)^{p(k-1)}(-1)^{p(1)}$.
We then have

$A^{k-1}e_{1}=(-1)^{q(k)}A^{\prime(-1)}A^{k-1}A^{\prime}e_{1}\quad\forall\,k\in\mathbb{N}.$

It holds that $A^{k-1}=A^{\prime(-1)}A^{k-1}A^{\prime}$ have the same
eigenvalues, so $A^{k-1}=-A^{\prime(-1)}A^{k-1}A^{\prime}$ must have
eigenvalues of opposite sign. That is, if $\lambda_{i}$ and
$\lambda_{i}^{\prime}$ are the eigenvalues of $A$ and $\pm
A^{\prime(-1)}AA^{\prime}$ respectively, then
$\lambda_{i}=\pm\lambda_{i}^{\prime}$. Following the same steps as in the
proof of Theorem 2 we get

$\sum_{i}\lambda_{i}^{k-1}v_{i}=(-1)^{q(k)}\sum_{i}\lambda_{i}^{k-1}v_{i}^{\prime}\eta_{i1}^{\prime
T}\quad\forall\,k\in\mathbb{N}.$

Subtracting the right hand side from both sides gives us

$\sum_{i}\lambda_{i}^{k-1}(v_{i}-(-1)^{q(k)}v_{i}^{\prime}\eta_{i1}^{\prime
T})=0\quad\forall\,k\in\mathbb{N}.$ (8)

We now show that if $(A,B)\in(\mathbb{R}^{2\times 2},\mathbb{R}^{2})$, then
equation (8) implies $A=\pm A^{\prime(-1)}AA^{\prime}$. Recall that
$q(k)\in\\{0,1\\}$ and so $(q(1),q(2))\in\\{(0,0),\,(0,1),\,(1,0),\,(1,1)\\}$.
When $(q(1),q(2))=(0,0)$, equation (8) is the same as equation (6) for $k=1$
and $k=2$. If we write these equations in matrix form as in Theorem 2 we again
have the Vandermonde matrix on the left-hand side. Following the same steps as
Theorem 2, we see that $v_{i}=v_{i}^{\prime}$ for all $i$. Since
$\lambda_{i}=\pm\lambda_{i}^{\prime}$, then $A=\pm A^{\prime(-1)}AA^{\prime}$.
Similarly, if $(q(1),q(2))=(1,1)$, if we follow the same procedure to find
$A=\pm A^{\prime(-1)}AA^{\prime}$.

The most interesting cases are when $(q(1),q(2))\in\\{(0,1),(1,0)\\}$. Let us
first consider $(q(1),q(2))=(1,0)$. Recall $q(k)\in\\{0,1\\}$, so if $q(3)=0$,
then $(q(2),q(3))=(0,0)$. If $(q(k),q(k+1))=(0,0)$ for some $k$, we then have

$\Lambda S_{j}=\begin{bmatrix}\lambda_{1}^{k}&\lambda_{2}^{k}\\\
\lambda_{1}^{k+1}&\lambda_{2}^{k+1}\end{bmatrix}\begin{bmatrix}v_{1j}-v^{\prime}_{1j}\eta_{11}^{\prime
T}\\\ v_{2j}-v^{\prime}_{2j}\eta_{21}^{\prime
T}\end{bmatrix}=\begin{bmatrix}0\\\ 0\end{bmatrix}.$ (9)

We note

$\mathrm{det}(\Lambda)=\lambda_{1}^{k}\lambda_{2}^{k}\begin{vmatrix}1&1\\\
\lambda_{1}&\lambda_{2}\end{vmatrix}\neq 0$

since we have two non-zero scalars multiplied by the non-zero Vandermonde
determinant in the case of distinct eigenvalues. Hence, $\Lambda$ as defined
in (9) is invertible and we again conclude that $A=\pm
A^{\prime(-1)}AA^{\prime}$.

We lastly consider cases where $q(k)$ is alternating, namely
$\\{q(k)\\}_{k=1}^{3}=(0,1,0)$ and $\\{q(k)\\}_{k=1}^{3}=(1,0,1)$. In the
former case, we have

$\Lambda S_{j}=\begin{bmatrix}1&1\\\
\lambda_{1}^{2}&\lambda_{2}^{2}\end{bmatrix}\begin{bmatrix}v_{1j}-v^{\prime}_{1j}\eta_{11}^{\prime
T}\\\ v_{2j}-v^{\prime}_{2j}\eta_{21}^{\prime
T}\end{bmatrix}=\begin{bmatrix}0\\\ 0\end{bmatrix}.$

The generic assumption that all eigenvalues are distinct modulo a sign implies
$\Lambda$ is invertible, thus we again find $v_{i}=v_{i}^{\prime}$ and thus
$A=\pm A^{\prime(-1)}AA^{\prime}$. By following the same steps, we arrive at
the same conclusion when $\\{q(k)\\}_{k=1}^{3}=(1,0,1)$.

We now have that $A=(-1)^{q(2)}A^{\prime(-1)}AA^{\prime}$. If $q(2)=0$, then
$A$ and $A^{\prime}$ commute. Using assumptions of the theorem statement, we
can conclude that $A$ and $A^{\prime}$ have the same eigenvectors [28]. If
$q(2)=1$, then $A=-A^{\prime(-1)}AA^{\prime}$ and so
$A^{2}=A^{\prime(-1)}A^{2}A^{\prime}$. Clearly, $A^{2}$ and $A^{\prime}$
commute, and by the theorem statement, $A^{2}$ has distinct eigenvalues, which
again implies that $A$ and $A^{\prime}$ share the same eigenvectors.

We now follow the same steps as in the latter part of the proof of Theorem 2.
Namely, we can diagonalize $A$ and $A^{\prime}$; taking the eigenvalue
expansion of equation (7) and multiplying both sides on the left by the matrix
of left eigenvectors gives us the series of equations

$\lambda_{i}^{k-1}=(-1)^{q(k)}\lambda_{i}^{\prime
k-1}\quad\forall~{}k\in\mathbb{N}.$

Since $q(2)=0$ or $q(2)=1$, then $\lambda_{i}=\lambda_{i}^{\prime}$ or
$\lambda_{i}=-\lambda_{i}^{\prime}$ for all $i$. Since both $A$ and
$A^{\prime}$ have the same eigenvectors and eigenvalues same to a sign, then
$A=\pm A^{\prime}$. ∎

Theorem 2 solves Problem 1 in the generic case where $\mathcal{U}\neq[-c,c]$
while Theorem 3 proves there exists a $\pm$-unique solution to Problem 1 in
the two-dimensional case where $\mathcal{U}=[-c,c]$. The proof of Theorem 3
drives our intuition for Conjecture 1 in general: intuitively, adding
dimensions to the system should not make it more likely that multiple generic
systems can produce the same reachable sets for all time, especially
considering no two such systems exist when the input set is asymmetric.
Formalizing this statement is left for future work.

We remark that if the system dynamics do not satisfy the assumptions of
Theorem 2 or Theorem 3, they might not be (uniquely or $\pm$-uniquely)
recoverable. However, using a slight perturbation of the reachable sets might
recover a generic approximation of the true dynamics. Doing so, however,
introduces challenges on the method of perturbing these sets. We leave such a
discussion for future work.

##  IV Solving for the System Dynamics

We ultimately want to use reachable sets to solve for the system dynamics.
Equation (2) of Theorem 1 already gives us a formula for calculating
$A^{i-1}b\mathcal{U}$ for all $i\in\mathbb{N}$, namely

$A^{i-1}b\mathcal{U}=\mathcal{R}(i,0)\ominus\mathcal{R}(i-1,0).$

In Theorem 2, we proved that the answer to Problem 1 is affirmative for
generic, single-input linear systems, meaning that for cases where the linear
system dynamics satisfy the generic assumptions of Lemma 1, we can uniquely
determine the true dynamics from the system’s reachable sets. This motivates
us to devise a procedure to calculate $(A,b)$.

We will determine $(A,b)$ from reachable sets through a two step procedure.
First, we calculate $A^{i-1}b\mathcal{U}$ for $i~{}=~{}\\{1,\ldots,n+1\\}$. In
the case where $\mathcal{U}\neq[-c,c]$, the sequence of sets
$A^{i-1}b\mathcal{U}$ can be used to calculate $(A,b)$ directly. If
$\mathcal{U}=[-c,c]$, these same sets can be utilized to compute a number of
candidate dynamics $(A,b)$ which satisfy $\mathcal{R}(i,0)$ for all $i$. To
determine which candidate solutions are correct, we compute the forward
reachable sets of (1) using all candidate $(A,b)$. By Theorem 3, in the two-
dimensional case, only two solutions $(A,b)$ and $(A^{\prime},b^{\prime})$
such that $(A,b)=-(A^{\prime},b^{\prime})$ will satisfy $\mathcal{R}(i,0)$ for
all $i$.

We begin our method by first using an algorithm that takes reachable sets of
(1) and solves for $A^{i-1}b\mathcal{U}$. By equation (2), we can utilize
existing methods [20, 22, 30] to compute the Minkowski difference between two
polygons to calculate $A^{i-1}b\mathcal{U}$ given $\mathcal{R}(i,0)$ for all
$i\in\mathbb{N}$. For this narrative, we adopt the method in [20]. By Lemma 1
of [20], if we let $v^{(i)}\in\mathcal{V}$ be the vertices of
$\mathcal{R}(i-1,0)$, then the Minkowski difference
$\mathcal{R}(i,0)\ominus\mathcal{R}(i-1,0)$ may be computed by taking the
intersection of the translation of the set $\mathcal{R}(i,0)$ by vertices
$v^{(i)}\in\mathcal{V}$ of $\mathcal{R}(i-1,0)$:

$\mathcal{R}(i,0)\ominus\mathcal{R}(i-1,0)=\bigcap_{v^{(i)}\in\mathcal{V}}(\mathcal{R}(i,0)-v^{(i)}).$
(10)

While computing the intersection in (10) is generally computationally
difficult, calculations are made significantly easier as $A^{i-1}b\mathcal{U}$
is a line segment; see [20] for details.

We now move to recover $A^{i-1}b$ from $A^{i-1}b\mathcal{U}$. We consider two
cases: $\mathcal{U}\neq[-c,c]$ and $\mathcal{U}=[-c,c]$ for some
$c\in\mathbb{R}$. In the former case, taking the mean of the vertices of
$A^{i-1}b\mathcal{U}$ will provide $A^{i-1}b\frac{c+d}{2}$. Multiplying this
vector by $\frac{2}{c+d}$ recovers $A^{i-1}b$.

###### Theorem 4

Let us assume the $n$-dimensional system (1) is controllable. Let
$C_{A,b}=\begin{bmatrix}b&Ab&\ldots&A^{n-1}b\end{bmatrix}$. For the single-
input case, $A=AC_{A,b}C_{A,b}^{-1}$.

The proof of Theorem 4 is trivial, noting that $C_{A,b}$ is full rank for
controllable systems. We note that the assumption of controllability is
generic [29].

In the case where $\mathcal{U}=[-c,c]$, by multiplying the vertices of
$A^{i-1}b\mathcal{U}$ by $c$, we can only recover $A^{i-1}b$ up to a sign,
generating two candidates for each $i$. Substituting all possible candidates
for $A^{i-1}b$ into the columns of $C_{A,b}$ and $AC_{A,b}$ generates
$2^{n+1}$ candidate matrices $A$.

To determine which candidate solutions yield the correct $\pm$-unique matrix
pair $(A,b)$, we can plot the reachable sets of all $2^{n+1}$ candidate
solutions to solve for the desired unknown $\pm$-unique system dynamics. In
the next section, we use the CORA toolkit [31] and adopt methods of computing
the Minkowski difference detailed in [20] to numerically calculate the
dynamics $(A,b)$ for an unknown band-pass filter circuit system and a two-
dimensional unknown system with $\mathcal{U}=[-1,1]$, validating the developed
theory.

##  V Numerical Examples

To validate the developed theory and demonstrate how to apply the proposed
method, we first consider a scenario of reverse engineering an electric
circuit from manufacturer specifications. At times, manufacturers will only
release partial information about a system. For example, instead of providing
a dynamic model of a manufactured part, manufacturers might convey the set of
all voltages a circuit may output within a set amount of time given the set of
all viable input frequencies. Such information can be interpreted as the
minimum time in which a state can be reached, providing a picture of the
system’s reachable sets. Motivated by such an example, in this section, we
provide an example of identifying the matrices $(A,b)$ of a band-pass filter
circuit from its reachable sets. In a subsequent example, we identify the
$\pm$-unique dynamics of an unknown two-dimensional system with an input set
symmetric around zero. Both examples utilize the CORA toolkit [31] for set
computations, namely to calculate convex hulls and Minkowski differences.

###  V-A Band-Pass Filter Circuit

We present the linear dynamic model of a band-pass filter circuit [32]. Let us
assume $x[0]=0$. The state-space controllable canonical representation [29] of
this circuit is

$\begin{gathered}x[j+1]=Ax[j]+bv_{c}[j]\\\ =\begin{bmatrix}0&1&0&0\\\
0&0&1&0\\\ 0&0&0&1\\\
-a_{0}&-a_{1}&-a_{2}&-a_{3}\end{bmatrix}x[j]+\begin{bmatrix}0\\\ 0\\\ 0\\\
1\end{bmatrix}v_{c}[j]\end{gathered}$ (11)

such that $v_{c}[j]\in[0,1]$ for all $j\in\mathbb{Z}_{\geq 0}$.

Assume the reachable sets $\\{\mathcal{R}(j+1,0)\\}_{j=0}^{\infty}$ of the
controllable, dynamical system (11) are known. From controllability, we know
the form of system (11), but not parameters $a_{0},\,a_{1},\,a_{2},\,a_{3}$.
From this information, we want to recover the true parameters: $a_{0}=3$,
$a_{1}=2$, $a_{2}=3$, and $a_{3}=6$. It can be easily shown that if $a_{0}\neq
0$, the assumptions of Theorem 2 are satisfied. Clearly, $A$ from (11)
satisfies said assumptions. Clearly, the matrix $M$ for which $Mb=e_{1}$ is a
simple row permutation; the assumptions of Theorem 2 are invariant under
permutations, hence all assumptions are satisfied and the results of Theorem 2
apply when solving for the matrix pair $(A,b)$. That is, there exists a unique
matrix pair which satisfies $\\{\mathcal{R}(j+1,0)\\}_{j=0}^{\infty}$. Since
the system is four-dimensional, Theorem 4 shows we need only consider the sets
$\\{\mathcal{R}(j+1,0)\\}_{j=0}^{4}$ to calculate $(A,b)$.

Assume, that $\\{\mathcal{R}(j+1,0)\\}_{j=0}^{4}$ are known to equal

$\mathcal{R}(1,0)=\mathrm{conv}\left(\begin{bmatrix}0\\\ 0\\\ 0\\\
0\end{bmatrix},\begin{bmatrix}0\\\ 0\\\ 0\\\ 1\end{bmatrix}\right),$
$\mathcal{R}(2,0)=\mathrm{conv}\left(\begin{bmatrix}0\\\ 0\\\ 1\\\
-5\end{bmatrix},\begin{bmatrix}0\\\ 0\\\ 1\\\
-6\end{bmatrix},\begin{bmatrix}0\\\ 0\\\ 0\\\
0\end{bmatrix},\begin{bmatrix}0\\\ 0\\\ 0\\\ 1\end{bmatrix}\right),$
$\mathcal{R}(3,0)=\mathrm{conv}\left(\begin{bmatrix}0\\\ 0.86\\\ -6.02\\\
33.00\end{bmatrix},\begin{bmatrix}0\\\ -0.14\\\ 0.98\\\
-6.00\end{bmatrix},\begin{bmatrix}0\\\ -0.14\\\ 0.98\\\
-5.00\end{bmatrix},\begin{bmatrix}0\\\ 0.86\\\ -6.02\\\
34.00\end{bmatrix}\right),$
$\mathcal{R}(4,0)=\mathrm{conv}\left(\begin{bmatrix}-0.15\\\ 1.05\\\ -5.99\\\
33.00\end{bmatrix},\begin{bmatrix}0.85\\\ -5.95\\\ 34.01\\\
-188\end{bmatrix},\begin{bmatrix}0.85\\\ -5.95\\\ 34.01\\\
-187\end{bmatrix},\begin{bmatrix}-0.15\\\ 1.05\\\ -5.99\\\
34.00\end{bmatrix}\right),$ $\mathcal{R}(5,0)=$
$\mathrm{conv}\left(\begin{bmatrix}-5.93\\\ 33.88\\\ -188.02\\\
1035.00\end{bmatrix},\begin{bmatrix}1.07\\\ -6.12\\\ 33.98\\\
-188.00\end{bmatrix},\begin{bmatrix}1.07\\\ -6.12\\\ 33.98\\\
-187.00\end{bmatrix},\begin{bmatrix}-5.93\\\ 33.88\\\ -188.02\\\
1036.00\end{bmatrix}\right).$

Based on Theorem 2, the knowledge of only these five sets is sufficient to
reconstruct the true values of parameters $a_{0}$, $a_{1}$, $a_{2}$, and
$a_{3}$.

Since $\mathcal{R}(1,0)=b\mathcal{U}$ and $\mathcal{U}=[0,1]$, $b$ can be
trivially computed to equal $b=\begin{bmatrix}0&0&0&1\end{bmatrix}^{T}$. Next,
by equation (2) of Theorem 1,
$Ab\mathcal{U}=\mathcal{R}(2,0)\ominus\mathcal{R}(1,0)$. Given
$\mathcal{U}=[0,1]$, by taking the Minkowski difference we get
$Ab=\begin{bmatrix}0&0&1&-6\end{bmatrix}^{T}$. Repeating the procedure, we
have $A^{2}b\mathcal{U}=\mathcal{R}(3,0)\ominus\mathcal{R}(2,0)$,
$A^{3}b\mathcal{U}=\mathcal{R}(4,0)\ominus\mathcal{R}(3,0)$,
$A^{4}b\mathcal{U}=\mathcal{R}(5,0)\ominus\mathcal{R}(4,0)$, and
$\mathcal{U}=[0,1]$. It follows that

$A^{2}b=\begin{bmatrix}0\\\ 1\\\ -6\\\
33\end{bmatrix},\,A^{3}b=\begin{bmatrix}1\\\ -6\\\ 33\\\
-182\end{bmatrix},\,A^{4}b=\begin{bmatrix}-6\\\ 33\\\ -182\\\
1002\end{bmatrix}.$

Recall we assume the system is controllable, and thus the controllability
matrix $C_{A,b}$ is invertible. Finally, by Theorem 4,

$A=AC_{A,b}C_{A,b}^{-1}$
$=\begin{bmatrix}A^{4}b&A^{3}b&A^{2}b&Ab\end{bmatrix}\begin{bmatrix}A^{3}b&A^{2}b&Ab&b\end{bmatrix}^{-1}$

which produces the correct matrix $A$ accurately reconstructing the parameters
$a_{0}=3$, $a_{1}=2$, $a_{2}=3$, and $a_{3}=6$.

###  V-B Numerical Example with a Symmetric Input Set

To validate Theorem 3, we present an example of a linear two-dimensional
dynamical system

$\begin{bmatrix}x_{1}[i+1]\\\ x_{2}[i+1]\end{bmatrix}=\begin{bmatrix}2&1\\\
2&3\end{bmatrix}\begin{bmatrix}x_{1}[i]\\\
x_{2}[i]\end{bmatrix}+\begin{bmatrix}0\\\ 1\end{bmatrix}u[i]$ (12)

with $\mathcal{U}=[-1,1]$. Such a system satisfies the assumptions of Lemma 1.
As in the previous example, we will show that we can reconstruct the values of
system matrices in (12) from reachable sets, albeit up to a sign. Assume,
thus, that we are given a sequence of reachable sets
$\\{\mathcal{R}(i,0)\\}_{i=1}^{4}$ as convex hulls of vertices:

$\mathcal{R}(1,0)=\mathrm{conv}\left(\pm\begin{bmatrix}0\\\
1\end{bmatrix}\right),$
$\mathcal{R}(2,0)=\mathrm{conv}\left(\pm\begin{bmatrix}-1\\\
-4\end{bmatrix},\pm\begin{bmatrix}1\\\ 2\end{bmatrix}\right),$
$\mathcal{R}(3,0)=\mathrm{conv}\left(\pm\begin{bmatrix}6\\\
15\end{bmatrix},\pm\begin{bmatrix}4\\\ 7\end{bmatrix},\pm\begin{bmatrix}6\\\
13\end{bmatrix}\right),$
$\mathcal{R}(4,0)=\mathrm{conv}\left(\pm\begin{bmatrix}27\\\
58\end{bmatrix},\pm\begin{bmatrix}15\\\
28\end{bmatrix},\pm\begin{bmatrix}25\\\
50\end{bmatrix},\pm\begin{bmatrix}27\\\ 56\end{bmatrix}\right).$

Clearly, $b\mathcal{U}=\mathcal{R}(1,0)$. Equation (2) of Theorem (1) also
shows that $Ab\mathcal{U}=\mathcal{R}(2,0)\ominus\mathcal{R}(1,0)$,
$A^{2}b\mathcal{U}=\mathcal{R}(3,0)\ominus\mathcal{R}(2,0)$, and
$A^{3}b\mathcal{U}=\mathcal{R}(4,0)\ominus\mathcal{R}(3,0)$. Since
$\mathcal{U}=[-1,1]$, through the same calculations in the previous example we
thus obtain

$b=\pm\begin{bmatrix}0\\\ 1\end{bmatrix},\,Ab=\pm\begin{bmatrix}1\\\
3\end{bmatrix},$ $A^{2}b=\pm\begin{bmatrix}5\\\
11\end{bmatrix},\,A^{3}b=\pm\begin{bmatrix}21\\\ 43\end{bmatrix}.$

Let us denote $b^{-}=\begin{bmatrix}0\\\ -1\end{bmatrix}$,
$b^{+}=\begin{bmatrix}0\\\ 1\end{bmatrix}$, $Ab^{-}=\begin{bmatrix}-1\\\
-3\end{bmatrix}$, etc. We now consider a set of $2^{3}$ possible candidate
pairs of $(C_{A,b},AC_{A,b})$ matrices:

$(C_{A,b},AC_{A,b})=\left\\{\begin{array}[]{cc}\left(\begin{bmatrix}b^{+}&Ab^{+}\end{bmatrix},\begin{bmatrix}Ab^{+}&A^{2}b^{+}\end{bmatrix}\right),\\\
\left(\begin{bmatrix}b^{-}&Ab^{+}\end{bmatrix},\begin{bmatrix}Ab^{+}&A^{2}b^{+}\end{bmatrix}\right),\\\
\left(\begin{bmatrix}b^{+}&Ab^{-}\end{bmatrix},\begin{bmatrix}Ab^{-}&A^{2}b^{+}\end{bmatrix}\right),\\\
\left(\begin{bmatrix}b^{-}&Ab^{-}\end{bmatrix},\begin{bmatrix}Ab^{-}&A^{2}b^{+}\end{bmatrix}\right),\\\
\left(\begin{bmatrix}b^{+}&Ab^{+}\end{bmatrix},\begin{bmatrix}Ab^{+}&A^{2}b^{-}\end{bmatrix}\right),\\\
\left(\begin{bmatrix}b^{-}&Ab^{+}\end{bmatrix},\begin{bmatrix}Ab^{+}&A^{2}b^{-}\end{bmatrix}\right),\\\
\left(\begin{bmatrix}b^{+}&Ab^{-}\end{bmatrix},\begin{bmatrix}Ab^{-}&A^{2}b^{-}\end{bmatrix}\right),\\\
\left(\begin{bmatrix}b^{-}&Ab^{-}\end{bmatrix},\begin{bmatrix}Ab^{-}&A^{2}b^{-}\end{bmatrix}\right).\end{array}\right.$
(13)

By Theorem 4, determining all candidate matrix pairs $(A,b)$ becomes a trivial
calculation using all possible pairs from (13). Doing so provides two
$\pm$-unique candidate pairs:

$(A,b)=\pm\left(\begin{bmatrix}2&1\\\ 2&3\end{bmatrix},\begin{bmatrix}0\\\
1\end{bmatrix}\right),\,(A^{\prime},b^{\prime})=\pm\left(\begin{bmatrix}8&-1\\\
20&-3\end{bmatrix},\begin{bmatrix}0\\\ -1\end{bmatrix}\right).$

(a)

(b)

Figure 1: Reachable sets for candidate dynamics at $i=4$.

While the calculations above used only $\mathcal{R}(1,0)$, $\mathcal{R}(2,0)$,
and $\mathcal{R}(3,0)$, to distinguish between these two final candidates we
need to employ $\mathcal{R}(4,0)$. Fig. 1 shows the plots of the forward
reachable sets for system (12) at time $i=4$ with matrix pairs $(A,b)$ and
$(A^{\prime},b^{\prime})$ on the left and right respectively.

Fig. 1(a) shows a reachable set that is identical to $\mathcal{R}(4,0)$, while
Fig. 1(b) illustrates the reachable set

$\mathcal{R}^{\prime}(4,0)=\mathrm{conv}\left(\pm\begin{bmatrix}35\\\
82\end{bmatrix},\,\pm\begin{bmatrix}25\\\
60\end{bmatrix},\,\pm\begin{bmatrix}33\\\
74\end{bmatrix},\,\pm\begin{bmatrix}35\\\ 80\end{bmatrix}\right),$

which is not the same as $\mathcal{R}(4,0)$. Therefore, we can identify the
matrix pair $(A,b)$, up to a sign, as the true dynamics of the unknown linear
system (12). As mentioned before, reachable sets in this case do not allow us
to distinguish any further: the dynamics that differ only in the sign generate
the same reachable sets.

##  VI Conclusion

This paper considers the problem of determining the dynamics of an unknown
discrete-time linear system using its reachable sets. The theory developed in
this paper proves that for input sets that are asymmetric around the origin,
the derived system dynamics are, given some technical assumptions, unique.
Thus, in such cases, we can determine the true dynamics of an unknown system
using the sequence of the system’s reachable sets. For the case where the
input set is symmetric, we prove that the derived dynamics are unique up to a
factor of $\pm 1$ for two-dimensional systems and provide a conjecture that
asserts the same result holds for the $n$-dimensional case. We then develop a
method for deriving the dynamics of a system given the sequence of the
system’s reachable sets using Minkowski differences and proceed to illustrate
by example how the method can be applied to identify the unknown linear model
of a band-pass filter. Novel identification methods are also applied to an
academic system with an input set symmetric around zero to detail how we can
adapt said methods to uniquely identify the model of a linear system.

A natural next step is to prove the stated conjecture to show $\pm$-uniqueness
for $n$-dimensional systems. Also, our current technical assumptions are
consistent with generic properties of matrices, but ideally we want to relax
these assumptions to identify necessary conditions for uniqueness. We also
want to consider cases when the state’s initial conditions are non-zero, when
there is only available knowledge of the system’s reachable sets at non-
consecutive time steps, and also when working with the more general framework
of multi-input systems.

## Acknowledgements

We thank Jeffrey Stuart from Pacific Lutheran University for providing
insights in combinatorial matrix theory that helped develop the scope this
project, namely addressing the question of uniqueness outlined in Problem 1.

## References

  * [1] G. Asaamoning, P. Mendes, D. Rosário, and E. Cerqueira, “Drone swarms as networked control systems by integration of networking and computing,” _Sensors_ , vol. 21, no. 8, p. 2642, 2021.
  * [2] I. Taub and T. Shima, “Intercept angle missile guidance under time varying acceleration bounds,” _Journal of Guidance, Control, and Dynamics_ , vol. 36, no. 3, pp. 686–699, 2013.
  * [3] S. R. Musse and D. Thalmann, “A model of human crowd behavior: Group inter-relationship and collision detection analysis,” in _Computer Animation and Simulation’97_. Springer, 1997, pp. 39–51.
  * [4] G. Schiebinger, J. Shu, M. Tabaka, B. Cleary, V. Subramanian, A. Solomon, S. Liu, S. Lin, P. Berube, L. Lee _et al._ , “Reconstruction of developmental landscapes by optimal-transport analysis of single-cell gene expression sheds light on cellular reprogramming,” _BioRxiv_ , 2017.
  * [5] S. Zhang, A. Afanassiev, L. Greenstreet, T. Matsumoto, and G. Schiebinger, “Optimal transport analysis reveals trajectories in steady-state systems,” _PLoS computational biology_ , vol. 17, no. 12, 2021.
  * [6] L. Ljung, “System identification,” in _Signal analysis and prediction_. Springer, 1998, pp. 163–173.
  * [7] R. Pintelon and J. Schoukens, _System identification: a frequency domain approach_. John Wiley & Sons, 2012.
  * [8] W. T. Miller, R. S. Sutton, and P. J. Werbos, _Neural networks for control_. MIT Press, 1995.
  * [9] S. Chen, S. A. Billings, and P. Grant, “Non-linear system identification using neural networks,” _International Journal of Control_ , vol. 51, no. 6, pp. 1191–1214, 1990.
  * [10] S. Lu and T. Basar, “Robust nonlinear system identification using neural network models,” in _34th IEEE Conference on Decision and Control_ , 1995, pp. 1840–1845.
  * [11] L. R. Medsker and L. Jain, “Recurrent neural networks,” _Design and Applications_ , vol. 5, pp. 64–67, 2001.
  * [12] W. Zaremba, I. Sutskever, and O. Vinyals, “Recurrent neural network regularization,” _arXiv preprint arXiv:1409.2329_ , 2014.
  * [13] C. Z. Aguilar, J. Gómez-Aguilar, V. Alvarado-Martínez, and H. Romero-Ugalde, “Fractional order neural networks for system identification,” _Chaos, Solitons & Fractals_, vol. 130, p. 109444, 2020\.
  * [14] C. Zúñiga-Aguilar, J. Gómez-Aguilar, H. Romero-Ugalde, H. Jahanshahi, and F. E. Alsaadi, “Fractal-fractional neuro-adaptive method for system identification,” _Engineering with Computers_ , pp. 1–24, 2021\.
  * [15] A. Ledezma, R. Aler, A. Sanchis, and D. Borrajo, “Predicting opponent actions by observation,” in _Robot Soccer World Cup_. Springer, 2004, pp. 286–296.
  * [16] V. Lisỳ, R. Píbil, J. Stiborek, B. Bošanskỳ, and M. Pěchouček, “Game-theoretic approach to adversarial plan recognition,” in _European Conference on Artificial Intelligence_ , 2012, pp. 546–551.
  * [17] N. Le Guillarme, A.-I. Mouaddib, X. Lerouvreur, and S. Gatepaille, “A generative game-theoretic framework for adversarial plan recognition,” in _10es Journées Francophones sur la Planification, la Décision et l’Apprentissage_ , 2015.
  * [18] J. A. Iglesias, P. Angelov, A. Ledezma, and A. Sanchis, “Evolving classification of agents’ behaviors: a general approach,” _Evolving Systems_ , vol. 1, no. 3, pp. 161–171, 2010.
  * [19] J. Zhang, M. Z. A. Bhuiyan, X. Yang, T. Wang, X. Xu, T. Hayajneh, and F. Khan, “Anticoncealer: reliable detection of adversary concealed behaviors in Edgeai Assisted Iot,” _IEEE Internet of Things Journal_ , 2021.
  * [20] M. Althoff, “On computing the Minkowski difference of zonotopes,” _arXiv preprint arXiv:1512.02794_ , 2015.
  * [21] M. Fradelizi, Z. Lángi, and A. Zvavitch, “Volume of the Minkowski sums of star-shaped sets,” _Proceedings of the American Mathematical Society, Series B_ , vol. 9, no. 34, pp. 358–372, 2022.
  * [22] Y. Yan and G. S. Chirikjian, “Closed-form characterization of the minkowski sum and difference of two ellipsoids,” _Geometriae Dedicata_ , vol. 177, no. 1, pp. 103–128, 2015.
  * [23] J. L. Kelley, _General topology_. Courier Dover Publications, 2017.
  * [24] T. Katō, _Perturbation theory for linear operators_. Springer, 1976, vol. 132.
  * [25] G. Strang, _Linear algebra and its applications._ Thomson, Brooks/Cole, 2006.
  * [26] F. Amato and M. Ariola, “Finite-time control of discrete-time linear systems,” _IEEE Transactions on Automatic control_ , vol. 50, no. 5, pp. 724–729, 2005.
  * [27] A. Klinger, “The Vandermonde matrix,” _American Mathematical Monthly_ , vol. 74, no. 5, pp. 571–574, 1967.
  * [28] G. Strang, _Introduction to Linear Algebra_. Wellesley-Cambridge Press, 2016.
  * [29] C.-T. Chen, _Linear system theory and design_. Saunders College Publishing, 1984.
  * [30] Y. Feng and Y. Tan, “On Minkowski difference-based contact detection in discrete/discontinuous modelling of convex polygons/polyhedra: Algorithms and implementation,” _Engineering Computations_ , vol. 37, no. 1, pp. 54 – 72, 2019.
  * [31] M. Althoff, D. Grebenyuk, and N. Kochdumper, “Implementation of Taylor models in CORA 2018,” in _5th International Workshop on Applied Verification for Continuous and Hybrid Systems_ , 2018, pp. 145 – 173.
  * [32] D. Denisenko, N. Prokopenko, N. Butyrlagin, and A. Titov, “Synthesis of a fourth-order bandpass filter circuit with considering the technological limitations on microelectronic realization,” in _2021 International Seminar on Electron Devices Design and Production_ , 2021.