An Improvement on Data-Driven Pole Placement for State Feedback Control and Model Identification

Show more

1. Introduction

In state feedback pole placement, the state feedback gain must be determined for a given system such that the closed-loop poles coincide with the desired locations. This is a well-known problem, and the pole placement methods have been extensively discussed in the literature [1] [2] [3] [4] . In standard pole placement methods, a state space model is assumed to be given by a system identification technique using data from past experiments. Whereas the traditional approach combines the identification of the state space model with the standard pole placement method; an alternative approach called “data-driven pole placement” has recently been proposed [5] . In this approach, the state space model and pole placement feedback gain are identified simultaneously from the set of state measurements and control input sequences. The method proposed in [5] is based on the data-driven control framework ( [6] and references therein) such as unfalsified control [7] , virtual reference feedback tuning (VRFT) [8] [9] , or fictitious reference iterative tuning (FRIT) [10] [11] [12] [13] . In the data-driven control framework, where no explicit mathematical plant model is used, a feedback controller must be derived that satisfies the prescribed closed-loop performance and fits to known experimental data. In contrast with traditional model-based controller designs, techniques such as controller identification [14] or a combination of plant model and controller identification must be applied [15] [16] .

Many studies of data-driven control have focused on output feedback control and data-driven state feedback control [11] [12] [13] , in which the prescribed closed-loop performance is achieved by applying a closed-loop reference transfer function. Such methods can be applied to the data-driven pole placement problem by choosing a reference transfer function with the desired poles. However, the zeros of the reference transfer function cannot normally be specified, because the zeros of the plant are unknown. In contrast, the data-driven pole placement method presented in [5] requires only a state space representation of the closed- loop system to specify the prescribed closed-loop performance, as shown in Section 2. This avoids the zero assignment issue that arises in the transfer function approach used in [5] .

This data-driven pole placement method can, therefore, be applied to linear and time-invariant systems with measurable states. The method is briefly reviewed in Section 2. However, the capacity of the data-driven pole placement method to handle noise remains an open issue, though in [5] , the total least square (TLS) method [17] was claimed to be eﬀective. Measurement noise is one of the issues which may surely face in practical applications. Therefore, to resolve this, we introduced a prefiltering technique that reduces the eﬀect of measurement noise in Section 3. More specifically, a finite impulse response (FIR) filter is used to prefilter the data, as this makes them easier to manipulate. In Section 4, by using the numerical example of a self-balancing robot, we discuss the eﬀect of applying this prefiltering technique, together with the least square (LS) and TLS methods, to a self-balancing robot model. We investigate the ability of the data-driven pole placement method to produce a linearized model using numerical simulations as in [18] . A nonlinear diﬀerential equation was used to represent the dynamics of a self-balancing robot there. Moreover, we evaluate the effects by two different exciting signals, the random and the chirp exciting signal, along with TLS and prefiltering. Finally, we compare all the results for the pole placement error and identification error when two exciting signals are applied.

Notation: Let $A$ and $B$ be $m\times n$ and $p\times q$ matrices, respectively. Then, the Kronecker product of $A$ and $B$ is a $mp\times nq$ matrix, defined as follow:

$A\otimes B=\left[\begin{array}{ccc}{a}_{11}B& \cdots & {a}_{1n}B\\ \vdots & & \vdots \\ {a}_{m1}B& \cdots & {a}_{mn}B\end{array}\right],$ (1)

where ${a}_{ij}(i=1,\cdots ,m,\text{}j=1,\cdots ,n)$ is the $i{j}^{th}$ element of $A$ . The vectorization of then stacks the columns into a vector:

$\text{vec}\left(A\right)=\left[\begin{array}{c}{a}_{1}\\ \vdots \\ {a}_{n}\end{array}\right],$ (2)

in which ${a}_{j}$ is the ${j}^{th}$ column of $A$ . The Frobenius norm of matrix $A\in {R}^{m\times n}$ is defined as

${\Vert A\Vert}_{F}=\sqrt{{\displaystyle \underset{i=1}{\overset{m}{\sum}}{\displaystyle \underset{j=1}{\overset{n}{\sum}}{\left|{a}_{ij}\right|}^{2}}}}.$ (3)

2. Data-Driven Pole Placement

In this section, we briefly review the data-driven pole placement method formulated in [5] .

Consider a discrete-time linear time-invariant system and static state feedback

$x\left(k+1\right)=Ax\left(k\right)+Bu\left(k\right)$ (4)

$u\left(k\right)=Fx\left(k\right)+v\left(k\right)$ (5)

where $A\in {\mathbb{R}}^{n\times n},B\in {\mathbb{R}}^{n\times m},x\in {\mathbb{R}}^{n}$ is the state vector, $u\in {\mathbb{R}}^{m}$ is the input vector, $F\in {\mathbb{R}}^{m\times n}$ is the feedback gain, and $v\in {\mathbb{R}}^{m}$ is the external input to the closed loop system.

The data-driven pole placement problem was formulated in [5] as follows:

Problem 1. We assume that the order of the plant n is known, state n is measurable, pair $\left(A,B\right)$ is controllable but the exact value is unknown and $B$ is of full rank. Let $\text{\Lambda}=\left\{{p}_{1},\cdots ,{p}_{n}\right\}$ be a self-conjugate set of n complex numbers in the unit circle. Given the input and output measurement data sequence $\left({x}_{0}\left(k\right),{u}_{0}\left(k\right)\right)$ of (4), find a state feedback gain $F$ from the observed data $\left({x}_{0}\left(k\right),{u}_{0}\left(k\right)\right)$ such that $\left\{{\lambda}_{i}\left(A+BF\right)\right\}=\text{\Lambda}.$

In a conventional approach, this problem is solved in two steps: $A$ and $B$ are identified from $\left({x}_{0}\left(k\right),{u}_{0}\left(k\right)\right)$ , then $F$ is derived using the standard pole placement algorithms. In contrast, the data-driven pole placement method solves the two steps simultaneously. To achieve this, the method uses the equivalency between the closed-loop system

$x\left(k+1\right)=\left(A+BF\right)x\left(k\right)+Bv\left(k\right),$ (6)

with the desired pole placement gain $F$ and

${x}_{\text{d}}\left(k+1\right)={A}_{\text{d}}{x}_{\text{d}}\left(k\right)+{B}_{\text{d}}v\left(k\right),$ (7)

${x}_{\text{d}}\left(k\right)=Tx\left(k\right),$ (8)

where $\left({A}_{\text{d}},{B}_{\text{d}}\right)$ with ${\lambda}_{i}\left({A}_{\text{d}}\right)={p}_{i}$ is an appropriate controllable pair. This equivalency requires the nonsingular matrix $T$ to exist. Then, we remove $v$ from (7) by using (5), to obtain

${x}_{\text{d}}\left(k+1\right)={A}_{\text{d}}{x}_{\text{d}}\left(k\right)+{B}_{\text{d}}u\left(k\right)-{B}_{\text{d}}Fx\left(k\right).$ (9)

Then, using (8), we obtain

$Tx\left(k+1\right)={A}_{\text{d}}Tx\left(k\right)+{B}_{\text{d}}u\left(k\right)-{B}_{\text{d}}Fx\left(k\right).$ (10)

If $\left({x}_{0}\left(k\right),{u}_{0}\left(k\right)\right)\text{\hspace{0.17em}}\left(k=i,\cdots ,i+N\right)$ satisfies (10),

$T{X}_{0}{P}_{1}={A}_{\text{d}}T{X}_{0}{P}_{2}+{B}_{\text{d}}{U}_{0}-{B}_{\text{d}}F{X}_{0},$ (11)

where

${X}_{0}=\left[{x}_{0}\left(i\right)\text{}{x}_{0}\left(i+1\right)\text{}\cdots \text{}{x}_{0}\left(i+N\right)\right],$ (12)

${U}_{0}=\left[{u}_{0}\left(i\right)\text{}{u}_{0}\left(i+1\right)\text{}\cdots \text{}{u}_{0}\left(i+N-1\right)\right],$ (13)

${P}_{1}=\left[\begin{array}{c}{0}_{1\times N}\\ {I}_{N}\end{array}\right],\text{}{P}_{2}=\left[\begin{array}{c}{I}_{N}\\ {0}_{1\times N}\end{array}\right].$ (14)

In [5] , Equation (11) is cast into

${S}_{1}\left[\begin{array}{c}T\\ F\end{array}\right]{X}_{0}{P}_{1}+{S}_{2}\left[\begin{array}{c}T\\ F\end{array}\right]{X}_{0}{P}_{2}={B}_{\text{d}}{U}_{0},$ (15)

${S}_{1}=\left[\begin{array}{cc}{I}_{n}& {0}_{n\times m}\end{array}\right],\text{}{S}_{2}=\left[\begin{array}{cc}-{A}_{\text{d}}& {B}_{\text{d}}\end{array}\right],$ (16)

and

$F=\left[\begin{array}{c}{f}_{1}\\ \vdots \\ {f}_{m}\end{array}\right]\in {\mathbb{R}}^{m\times n},\text{}T=\left[\begin{array}{c}{t}_{1}\\ \vdots \\ {t}_{m}\end{array}\right]\in {\mathbb{R}}^{n\times n}.$ (17)

Remark 1. The system in (7) can be interpreted as a reference model within VRFT (e.g., [8] [9] ) and FRIT (e.g., [10] [11] [12] [13] ). The idea of eliminating $v$ in (9) is also based on FRIT. In [10] [11] [12] , a similar state feedback control problem has been discussed within the FRIT framework. To apply these FRIT techniques to the data-driven pole placement problem, the desired transfer function must be specified from $u$ to $x$ , rather than ${x}_{\text{d}}$ . When precise values for $\left(A,B\right)$ are not available, it becomes impossible to specify the zeros of the desired transfer function.

Remark 2. To obtain the datasets in (12) by applying state feedback in (5) to the system in (4), the initial feedback gain $F$ should be based on $\left(A,B\right)$ . Hence, in Problem 1, the exact value of $\left(A,B\right)$ is assumed to be unknown.

When applying the property of Kronecker product $\text{vec}\left(MND\right)=\left({N}^{T}\otimes M\right)\text{vec}D$ (see for example Th.2.13 in [19] ) to the transpose of (15) to solve (15) for $F$ and $T$ , a further linear equation is derived, as follows:

${\rm X}\eta =U,$ (18)

where

$\eta ={\left[{t}_{1}\text{}\cdots \text{}{t}_{n}\text{}{f}_{1}\text{}\cdots \text{}{f}_{m}\right]}^{{\rm T}}\in {\mathbb{R}}^{\left(n+m\right)n},$ (19)

${\rm X}={S}_{1}\otimes {\left({X}_{0}{P}_{1}\right)}^{{\rm T}}+{S}_{2}\otimes {\left({X}_{0}{P}_{2}\right)}^{{\rm T}}\in {\mathbb{R}}^{nN\times \left(n+m\right)n},$ (20)

$U=\left({B}_{\text{d}}\otimes {U}_{0}^{{\rm T}}\right)\left(\text{vec}{I}_{m}\right)\in {\mathbb{R}}^{nN}.$ (21)

If $T$ is nonsingular, the model coefficients can be obtained

$A={T}^{-1}{A}_{\text{d}}T-{T}^{-1}{B}_{\text{d}}F,\text{}B={T}^{-1}{B}_{\text{d}}.$ (22)

3. Prefiltering Noisy Measurement

When the measurement of $x$ is contaminated by noise $\epsilon $ ,

${x}_{0}\left(k\right)=x\left(k\right)+\epsilon \left(k\right).$ (23)

Then, (10) becomes

$\begin{array}{l}T\left({x}_{0}\left(k+1\right)-\epsilon \left(k+1\right)\right)\\ ={A}_{\text{d}}T\left({x}_{0}\left(k\right)-\epsilon \left(k\right)\right)+{B}_{\text{d}}{u}_{0}\left(k\right)-{B}_{\text{d}}F\left({x}_{0}\left(k\right)-\epsilon \left(k\right)\right).\end{array}$ (24)

Hence, if $\left({x}_{0}\left(k\right),{u}_{0}\left(k\right)\right)\text{\hspace{0.17em}}\left(k=i,\cdots ,i+N\right)$ satisfies the above equation,

$\begin{array}{l}T\left({X}_{0}-E\right){P}_{1}\\ ={A}_{\text{d}}T\left({X}_{0}-E\right){P}_{2}+{B}_{\text{d}}{U}_{0}-{B}_{\text{d}}F\left({X}_{0}-E\right){P}_{2},\end{array}$ (25)

where

$E=\left[\epsilon \left(i\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\epsilon \left(i+1\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\cdots \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\epsilon \left(i+N\right)\right].$ (26)

Then, the resulting linear equation is given as

$\left({\rm X}+\Delta {\rm X}\right)\eta =U+\Delta U,$ (27)

where the effect of noise $\Delta {\rm X}$ has the same structure as ${\rm X}$ in (20), then

$\Delta {\rm X}=-{S}_{1}\otimes {\left(E{P}_{1}\right)}^{\text{T}}-{S}_{2}\otimes {\left(E{P}_{2}\right)}^{\text{T}},$ (28)

and $\Delta U$ is the equation error. Following [5] , we can solve $\eta \in {R}^{\left(n+m\right)n}$ to (27) as a TLS problem [17] , by minimizing the Frobenius norm ${\Vert \left[\begin{array}{cc}\Delta {\rm X}& \Delta U\end{array}\right]\Vert}_{F}$ . It is known that the TLS solution is given as

$\eta =-\frac{1}{{V}_{22}}{V}_{12}$ (29)

based on the singular value decomposition

$\left[{\rm X}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}U\right]=\left[{U}_{1}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{U}_{2}\right]\left[\begin{array}{cc}{\Sigma}_{1}& 0\\ 0& {\Sigma}_{2}\end{array}\right]{\left[\begin{array}{cc}{V}_{11}& {V}_{12}\\ {V}_{21}& {V}_{22}\end{array}\right]}^{\text{T}},$ (30)

where these matrices are partitioned into blocks corresponding to ${\rm X}$ and $U$ .

Here, we assume that there exists $M>0$ such that

$\frac{1}{M}{\displaystyle \underset{j=1}{\overset{M}{\sum}}\epsilon \left(i+j\right)}\approx 0$ (31)

for all $i$ . This means that when $N>M,$

$E{P}_{1}\Phi \approx 0,\text{}E{P}_{2}\Phi \approx 0$ (32)

for the matrix

$\Phi =\frac{1}{M}\left[\begin{array}{cccc}1& 0& \cdots & 0\\ \vdots & \ddots & \ddots & \vdots \\ 1& \ddots & \ddots & 0\\ 0& \ddots & \ddots & 1\\ \vdots & \ddots & \ddots & \vdots \\ 0& \cdots & 0& 1\end{array}\right]\in {\mathbb{R}}^{N\times \left(N-M+1\right)},$ (33)

where each column has $M$ elements of 1. Therefore,

$T{\stackrel{\u02dc}{X}}_{0}{P}_{1}={A}_{\text{d}}T{\stackrel{\u02dc}{X}}_{0}{P}_{2}+{B}_{\text{d}}{\stackrel{\u02dc}{U}}_{0}-{B}_{\text{d}}F{\stackrel{\u02dc}{X}}_{0}{P}_{2}$ (34)

where

${\stackrel{\u02dc}{X}}_{0}={X}_{0}\Phi ,\text{}{\stackrel{\u02dc}{U}}_{0}={U}_{0}\Phi .$ (35)

This multiplication by $\Phi $ represents the prefiltering of signals via an $M\text{th}$ order FIR filter.

When the systems (4) and (7) are driven by the exciting signal, we have

$\left({X}_{0}-E\right){P}_{1}=A\left({X}_{0}-E\right){P}_{2}+B{U}_{0},$ (36)

${U}_{0}=F\left({X}_{0}-E\right){P}_{2}+V,$ (37)

${X}_{\text{d}}=T\left({X}_{0}-E\right){P}_{2},$ (38)

${X}_{\text{d}}{P}_{1}={A}_{\text{d}}{X}_{\text{d}}{P}_{2}+{B}_{\text{d}}V,$ (39)

where

${X}_{\text{d}}=\left[{x}_{\text{d}}\left(i\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{x}_{\text{d}}\left(i+1\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\cdots \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{x}_{\text{d}}\left(i+N\right)\right],$ (40)

$V=\left[v\left(i\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}v\left(i+1\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\cdots \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}v\left(i+N-1\right)\right].$ (41)

By applying $\Phi $ to these systems, we obtain

${X}_{0}{P}_{1}\Phi =A{X}_{0}{P}_{2}\Phi +B{U}_{0}\Phi ,$ (42)

${U}_{0}\Phi =F{X}_{0}{P}_{2}\Phi +V\Phi ,$ (43)

${X}_{\text{d}}{P}_{2}\Phi =T{X}_{0}{P}_{2}\Phi ,$ (44)

${X}_{\text{d}}{P}_{1}\Phi ={A}_{\text{d}}{X}_{\text{d}}{P}_{2}\Phi +{B}_{\text{d}}V\Phi .$ (45)

Here, if $V\Phi \approx 0$ , (34) cannot be satisfied. Hence, for all i, $V\Phi \ne 0$ , that is

$\frac{1}{M}{\displaystyle \underset{j=1}{\overset{M}{\sum}}v\left(i+j\right)}\ne 0$ (46)

must be satisfied.

4. Numerical Example: Self-Balancing Robot

We next applied the data-driven pole placement method described above to the model of a self-balancing robot [21] [20] as shown in Figure 1. The robot is equipped with right and left wheels driven by direct current (DC) motors whose voltages ${v}_{r}$ and ${v}_{1}$ can be controlled. Because the motion dynamics can be decomposed by the input $u$ , the control input to the robot was represented as

$u=\left[\begin{array}{c}{u}_{1}\\ {u}_{2}\end{array}\right]=\left[\begin{array}{c}{v}_{r}+{v}_{1}\\ {v}_{r}-{v}_{1}\end{array}\right].$ (47)

Figure 1. (a) Coordinates of the self-balancing robot; (b) photo.

We assume that the pitch angle ${\theta}_{\text{b}}$ and the pitch angular velocity ${\stackrel{\dot{}}{\theta}}_{\text{b}}$ of the body could be measured, as well as the angles ${\theta}_{\text{r}}$ and ${\theta}_{\text{1}}$ of the right and left wheels, and their angular velocities ${\stackrel{\dot{}}{\theta}}_{\text{r}}$ and ${\stackrel{\dot{}}{\theta}}_{1}$ , respectively. We define the mean values of the right and left wheel angles ${\theta}_{\text{r}}$ and ${\theta}_{1}$ , and the yaw angle of the body as follows:

${\theta}_{\text{w}}=\frac{1}{2}\left({\theta}_{\text{r}}+{\theta}_{1}\right),$ (48)

$\varphi =\frac{r}{w}\left({\theta}_{\text{r}}-{\theta}_{1}\right),$ (49)

where $r$ is the radius of the wheel and $w=2d$ is the distance between the two wheels.

4.1. Equation of Motion

The equation of motion for the self-balancing robot can be derived as

$\{\begin{array}{l}{J}_{1}\left({\theta}_{\text{b}}\left(t\right)\right)\left[\begin{array}{c}{\stackrel{\xa8}{\theta}}_{\text{w}}\left(t\right)\\ {\stackrel{\xa8}{\theta}}_{\text{b}}\left(t\right)\end{array}\right]+{D}_{1}\left[\begin{array}{c}{\stackrel{\dot{}}{\theta}}_{\text{w}}\left(t\right)\\ {\stackrel{\dot{}}{\theta}}_{\text{b}}\left(t\right)\end{array}\right]-{M}_{\text{b}}l\mathrm{sin}{\theta}_{\text{b}}\left(t\right)\left[\begin{array}{c}r{\stackrel{\dot{}}{\theta}}_{\text{b}}^{2}\left(t\right)\\ g+l\mathrm{cos}{\theta}_{\text{b}}\left(t\right){\stackrel{\dot{}}{\varphi}}^{2}\left(t\right)\end{array}\right]={H}_{1}u\left(t\right)\\ {J}_{2}\left({\theta}_{\text{b}}\left(t\right)\right)\stackrel{\xa8}{\varphi}\left(t\right)+{D}_{2}\stackrel{\dot{}}{\varphi}\left(t\right)+\left(2{M}_{\text{b}}{l}^{2}\mathrm{sin}{\theta}_{\text{b}}\left(t\right)\mathrm{cos}{\theta}_{\text{b}}\left(t\right)\right){\stackrel{\dot{}}{\theta}}_{\text{b}}\left(t\right)\stackrel{\dot{}}{\varphi}\left(t\right)={H}_{2}u\left(t\right),\end{array}$ (50)

where

${J}_{1}\left({\theta}_{\text{b}}\right)=\left[\begin{array}{cc}2\left({J}_{\text{w}}+{g}_{\text{r}}^{2}{J}_{\text{m}}\right)+\left(2{M}_{\text{w}}+{M}_{\text{b}}\right){r}^{2}& -2{g}_{\text{r}}^{2}{J}_{\text{m}}+{M}_{\text{b}}lr\mathrm{cos}{\theta}_{\text{b}}\\ -2{g}_{\text{r}}^{2}{J}_{\text{m}}+{M}_{\text{b}}lr\mathrm{cos}{\theta}_{\text{b}}& {J}_{\text{b}}+2{g}_{\text{r}}^{2}{J}_{\text{m}}+{M}_{\text{b}}{l}^{2}\end{array}\right],$

${J}_{2}\left({\theta}_{\text{b}}\right)={J}_{\varphi}+2\frac{{d}^{2}}{{r}^{2}}\left({J}_{\text{w}}+{g}_{\text{r}}^{2}{J}_{\text{m}}\right)+2{M}_{\text{w}}{d}^{2}+{M}_{\text{b}}{l}^{2}{\mathrm{sin}}^{2}{\theta}_{\text{b}},$

${D}_{1}=2\left[\begin{array}{cc}{d}_{\text{b}}+{d}_{\text{w}}& -{d}_{\text{b}}\\ -{d}_{\text{b}}-{d}_{\text{w}}& {d}_{\text{b}}\end{array}\right],$ ${D}_{2}=2\frac{{d}^{2}}{{r}^{2}}\left({d}_{\text{b}}+{d}_{\text{w}}\right),$

${H}_{1}={b}_{\text{v}}\left[\begin{array}{cc}1& 0\\ -1& 0\end{array}\right],$ ${H}_{2}={b}_{\text{v}}\frac{d}{r}\left[\begin{array}{cc}0& 1\end{array}\right],$

${d}_{\text{b}}:=\frac{{g}_{\text{r}}^{2}{K}_{\text{t}}{K}_{\text{e}}}{{R}_{\text{m}}}+{g}_{\text{r}}{d}_{\text{m}},$ ${b}_{\text{v}}:=\frac{{g}_{\text{r}}{K}_{\text{t}}}{{R}_{\text{m}}}.$

The symbols are explained in Table 1. The parameters used in the simulations were taken from [20] [21] .

Table 1. Parameters of the self-balancing robot [20] [21] .

4.2. Linear Model and Feedback Gain

We linearized the equations of motion (50) around equilibrium states ${\theta}_{\text{w}}=0$ , ${\theta}_{\text{b}}=0$ , $\varphi =0$ , ${\stackrel{\dot{}}{\theta}}_{\text{w}}=0$ , ${\stackrel{\dot{}}{\theta}}_{\text{b}}=0$ , $\stackrel{\dot{}}{\varphi}=0$ , and $u=0$ . Then, under the assumption that $\mathrm{sin}{\theta}_{\text{b}}\left(t\right)\approx {\theta}_{\text{b}}\left(t\right)$ , $\mathrm{cos}{\theta}_{\text{b}}\left(t\right)\approx 1$ , ${\mathrm{sin}}^{2}{\theta}_{\text{b}}\left(t\right)\approx 0$ , ${\theta}_{\text{b}}^{2}\left(t\right)\approx 0$ , ${\varphi}^{2}\left(t\right)\approx 0$ , and $\mathrm{sin}{\theta}_{\text{b}}\left(t\right)\mathrm{cos}{\stackrel{\dot{}}{\theta}}_{\text{b}}\left(t\right)\stackrel{\dot{}}{\varphi}\left(t\right)\approx 0$ , the linearized equations of motion can be derived as

$\{\begin{array}{l}{J}_{1}{\stackrel{\xa8}{x}}_{a}\left(t\right)+{D}_{1}{\stackrel{\dot{}}{x}}_{a}\left(t\right)+{K}_{1}{x}_{a}\left(t\right)={H}_{1}u\left(t\right)\\ {J}_{2}{\stackrel{\xa8}{x}}_{b}\left(t\right)+{D}_{2}{\stackrel{\dot{}}{x}}_{b}\left(t\right)={H}_{2}u\left(t\right),\end{array}$ (51)

where

${x}_{a}\left(t\right):=\left[\begin{array}{c}{\theta}_{\text{w}}\left(t\right)\\ {\theta}_{\text{b}}\left(t\right)\end{array}\right],\text{\hspace{0.17em}}\text{\hspace{0.17em}}{x}_{\text{b}}\left(t\right):=\varphi \left(t\right),$ (52)

${J}_{1}={J}_{1}\left(0\right),\text{\hspace{0.17em}}{J}_{2}={J}_{2}\left(0\right),\text{\hspace{0.17em}}{K}_{1}={M}_{\text{b}}\mathrm{lg}\left[\begin{array}{cc}0& 0\\ 0& -1\end{array}\right].$ (53)

By defining the state vector

${x}_{1}\left(t\right)=\left[\begin{array}{c}{x}_{a}\left(t\right)\\ {\stackrel{\dot{}}{x}}_{a}\left(t\right)\end{array}\right]=\left[\begin{array}{c}{\theta}_{\text{w}}\left(t\right)\\ {\theta}_{\text{b}}\left(t\right)\\ {\stackrel{\dot{}}{\theta}}_{\text{w}}\left(t\right)\\ {\stackrel{\dot{}}{\theta}}_{\text{b}}\left(t\right)\end{array}\right],$ ${x}_{2}\left(t\right)=\left[\begin{array}{c}{x}_{b}\left(t\right)\\ {\stackrel{\dot{}}{x}}_{b}\left(t\right)\end{array}\right]=\left[\begin{array}{c}\varphi \left(t\right)\\ \stackrel{\dot{}}{\varphi}\left(t\right)\end{array}\right],$ (54)

the linear state space model can be derived as

$\{\begin{array}{l}{\stackrel{\dot{}}{x}}_{1}\left(t\right)={A}_{c1}{x}_{1}\left(t\right)+{B}_{c}{}_{1}{u}_{1}\left(t\right),\\ {\stackrel{\dot{}}{x}}_{2}\left(t\right)={A}_{c2}{x}_{2}\left(t\right)+{B}_{c2}{u}_{2}\left(t\right),\end{array}$ (55)

where

${A}_{c1}=\left[\begin{array}{cc}{0}_{2\times 2}& {I}_{2}\\ -{J}^{-1}{K}_{1}& -{J}^{-1}{D}_{1}\end{array}\right],$ ${B}_{c1}=\left[\begin{array}{c}{0}_{2\times 1}\\ {b}_{\text{v}}{J}_{1}{}^{-1}\left[\begin{array}{c}1\\ -1\end{array}\right]\end{array}\right],$

${A}_{c2}=\left[\begin{array}{cc}0& 1\\ 0& -{J}_{2}{}^{-1}{D}_{2}\end{array}\right],$ ${B}_{c2}=\left[\begin{array}{c}0\\ {b}_{\text{v}}\frac{d}{r}{J}_{2}{}^{-1}\end{array}\right].$

Then, the feedback can be independently designed as

${u}_{1}={F}_{1}{x}_{1}+{v}_{1},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{u}_{2}={F}_{2}{x}_{2}+{v}_{2}.$ (56)

Note that this can be more succinctly represented as

$\stackrel{\dot{}}{x}\left(t\right)={A}_{c}x\left(t\right)+{B}_{c}u\left(t\right),\text{}u\left(t\right)=Fx\left(t\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}x\left(t\right)=\left[\begin{array}{c}{x}_{1}\left(t\right)\\ {x}_{2}\left(t\right)\end{array}\right],$ (57)

${A}_{c}=\left[\begin{array}{cc}{A}_{c1}& {0}_{4\times 2}\\ {0}_{2\times 4}& {A}_{c2}\end{array}\right],\text{}{B}_{c}=\left[\begin{array}{cc}{B}_{c1}& {0}_{4\times 1}\\ {0}_{2\times 1}& {B}_{c2}\end{array}\right],\text{}F=\left[\begin{array}{cc}{F}_{1}& {0}_{1\times 2}\\ {0}_{1\times 4}& {F}_{2}\end{array}\right].$ (58)

When the parameters in Table 1 are used and the sampling period is $h=0.1\text{\hspace{0.17em}}\text{s}$ , the discrete-time model after discretizing (55) is

$\{\begin{array}{l}{x}_{1}\left(k+1\right)={A}_{1}{x}_{1}\left(k\right)+{B}_{1}{u}_{1}\left(k\right),\\ {x}_{2}\left(k+1\right)={A}_{2}{x}_{2}\left(k\right)+{B}_{2}{u}_{2}\left(k\right),\end{array}$ (59)

where

${A}_{1}=\left[\begin{array}{cccc}1& \text{0}\text{.1719}& \text{0}\text{.0226}& \text{0}\text{.0830}\\ 0& \text{1}\text{.1722}& \text{0}\text{.0113}& \text{0}\text{.0944}\\ 0& \text{3}\text{.5363}& \text{0}\text{.1388}& \text{1}\text{.0332}\\ 0& \text{3}\text{.5299}& \text{0}\text{.1386}& \text{1}\text{.0336}\end{array}\right],\text{\hspace{0.17em}}{B}_{1}=\left[\begin{array}{c}\text{0}\text{.0641}\\ -\text{0}\text{.0094}\\ \text{0}\text{.7131}\\ -\text{0}\text{.1148}\end{array}\right],$

${A}_{2}=\left[\begin{array}{cc}1& \text{0}\text{.0113}\\ 0& \text{0}\text{.0001}\end{array}\right],\text{\hspace{0.17em}}{B}_{2}=\left[\begin{array}{c}\text{0}\text{.0306}\\ \text{0}\text{.3450}\end{array}\right].$ (60)

Here, we assume that the exact values of (60) are not available, but that uncertain values are available:

${A}_{1}=\left[\begin{array}{cccc}1& \text{0}\text{.1897}& \text{0}\text{.0218}& \text{0}\text{.0844}\\ 0& \text{1}\text{.1900}& \text{0}\text{.0115}& \text{0}\text{.0947}\\ 0& \text{3}\text{.9115}& \text{0}\text{.1408}& \text{1}\text{.0489}\\ 0& \text{3}\text{.9151}& \text{0}\text{.1407}& \text{1}\text{.0492}\end{array}\right],\text{\hspace{0.17em}}{B}_{1}=\left[\begin{array}{c}\text{0}\text{.0648}\\ -\text{0}\text{.0095}\\ \text{0}\text{.7115}\\ -\text{0}\text{.1165}\end{array}\right],$

${A}_{2}=\left[\begin{array}{cc}1& \text{0}\text{.0103}\\ 0& \text{0}\text{.0001}\end{array}\right],\text{\hspace{0.17em}}{B}_{2}=\left[\begin{array}{c}\text{0}\text{.0310}\\ \text{0}\text{.3450}\end{array}\right].$ (61)

The coefficients can be derived from ${J}_{1},{J}_{2}$ , with an assumed uncertainty of 10%. By applying linear quadratic optimal control theory to (61), the desired closed-loop pole locations can be chosen as

$\lambda \left({A}_{1}+{B}_{1}{F}_{1}\right)\in {\Lambda}_{1}=\left\{6.0355\times {10}^{-5},\text{0}\text{.5253,0}\text{.5745,0}\text{.7630}\right\}$ , (62)

$\lambda \left({A}_{2}+{B}_{2}{F}_{2}\right)\in {\Lambda}_{2}=\left\{6.0426\times {10}^{-5},\text{0}\text{.7835}\right\}$ , (63)

and the initial feedback gains needed to obtain datasets for the data-driven pole placement as

${F}_{1}=\left[1.5216\text{}124.181\text{}2.3915\text{}18.3089\right],\text{\hspace{0.17em}}{F}_{2}=\left[-6.2764-0.0646\right].$ (64)

4.3. Comparison of Methods

Next, simulations were conducted and comparisons were made from the obtained results when using different methods and exciting signals.

Measurement noise was prepared with the Gaussian distribution $N\left(0,{\sigma}^{2}\right)$ , where ${\sigma}^{2}=1.0\times {10}^{-3}$ , $1.0\times {10}^{-4}$ and $1.0\times {10}^{-4}$ in ${\theta}_{\text{w}}$ , ${\stackrel{\dot{}}{\theta}}_{\text{b}}$ , and $\varphi $ , respectively. This is shown in Figure 2(a). We used the random exciting signal $v$ shown in Figure 3(a) and the linear chirp signal $v\left(k\right)$ shown in Figure 3(b) with the uniform distribution ${v}_{1}\left(k\right)\sim U\left(-0.5,0.5\right)$ and ${v}_{2}\left(k\right)\sim U\left(-0.1,0.1\right)$ . We set the order of the prefilter $\Phi $ (33) as $M=6$ . After prefiltering, the measurement noise in ${\theta}_{\text{w}}$ , ${\stackrel{\dot{}}{\theta}}_{\text{b}}$ , and $\varphi $ was reduced, as shown in Figure 2(b). The prefiltered exciting signals were shown in Figure 3(b) and Figure 3(d). It can be seen that the exciting signals $v$ were not eliminated by prefilter $\Phi $ , but that the high-frequency elements were reduced.

A closed-loop response in the presence of measurement noise by state feedback (56), with initial gain (64), is shown in Figure 4. The response to the random exciting signal and the chirp exciting signal are shown in Figure 4(a) and Figure 4(b), respectively. Of particular note is that the responses of ${\theta}_{\text{b}}$ , ${\stackrel{\dot{}}{\theta}}_{\text{w}}$ , and ${\stackrel{\dot{}}{\theta}}_{\text{b}}$ in Figure 4(b) show the high-pass filter-like gain characteristics of the transfer function from $v$ to $x$ .

For comparison, the dataset for the data-driven pole placement was chosen as ${\left\{\left({x}_{0}\left(k\right),{u}_{0}\left(k\right)\right)\right\}}_{k=50,\cdots ,450}$ where $i=50$ and $N=400$ .

To evaluate the obtained pole placement gain $\stackrel{\u02dc}{F}$ , we introduced an accuracy measurement that takes the largest absolute difference in value between each eigenvalue of ${A}_{i}+{B}_{i}{\stackrel{\u02dc}{F}}_{i}$ and the corresponding ${p}_{j}\in {\Lambda}_{i}$ ,

$\delta \lambda \left({A}_{\text{d}i}\right):=\mathrm{max}\left\{\left|{\lambda}_{j}\left({A}_{i}+{B}_{i}{\stackrel{\u02dc}{F}}_{i}\right)-{p}_{j}\right|{p}_{j}\in {\Lambda}_{i}\right\}.$ (65)

Figure 2. (a) Measurement noise; (b) Prefiltered measurement noise.

Figure 3. Exciting signal $v$ (a) random, (b) chirp, (c) prefiltered random, (d) prefiltered chirp.

Figure 4. Closed-loop response by an initial state feedback via (a) random exciting signal $v,$ (b) chirp exciting signal $v$ .

To evaluate the obtained model $\left(\stackrel{\u02dc}{A},\stackrel{\u02dc}{B}\right)$ , the following identification errors were used:

$\Delta {A}_{i}:=\Vert {\stackrel{\u02dc}{A}}_{i}-{A}_{i}\Vert ,\text{\hspace{0.17em}}\Delta {B}_{i}=\Vert {\stackrel{\u02dc}{B}}_{i}-{B}_{i}\Vert ,$ (66)

$\delta \lambda \left({A}_{i}\right):=\mathrm{max}\left\{\left|{\lambda}_{j}\left({A}_{i}\right)-{\lambda}_{j}\left({\stackrel{\u02dc}{A}}_{i}\right)\right|\right\}.$ (67)

The eigenvalues ${\lambda}_{j}$ were sorted by magnitude using the MATLAB command “sort”. This further sorts elements of equal magnitude by the phase angle on the interval $\left(-\text{\pi},\text{\pi}\right]$ . The impulse response $G\left(z\right)={\left(zI-\stackrel{\u02dc}{A}\right)}^{-1}\stackrel{\u02dc}{B}$ was used to evaluate the model obtained, as follows:

$\Delta {G}_{i}:=\sqrt{{\displaystyle \underset{k=0}{\overset{10}{\sum}}{\Vert {x}_{{\stackrel{\u02dc}{A}}_{i},{\stackrel{\u02dc}{B}}_{i}}\left(k\right)-{x}_{{A}_{i},{B}_{i}}\left(k\right)\Vert}^{2}}},$ (68)

where ${x}_{{\stackrel{\u02dc}{A}}_{i},{\stackrel{\u02dc}{B}}_{i}}$ and ${x}_{{A}_{i},{B}_{i}}$ are the impulse responses of ${G}_{i}\left(z\right):={\left(zI-{\stackrel{\u02dc}{A}}_{i}\right)}^{-1}{\stackrel{\u02dc}{B}}_{i}$ and $G\left(z\right):={\left(zI-A\right)}^{-1}B$ , respectively.

From the perspective of system control, smaller is better, particularly in the case of $\delta \lambda \left({A}_{\text{d}i}\right)$ , $\delta \lambda \left({A}_{i}\right)$ , and $\Delta {G}_{i}$ . The following key results were contrastively found in Table 2:

1) The initial model and feedback gain were affected by uncertainty: The model errors and pole placement errors are shown in Table 2 (initial).

2) The results when using the LS method to solve linear Equation (27) for noiseless data are shown in Table 2(a). All errors were reasonably small, confirming that the data-driven method performs well when the measurement data $\left({x}_{0}\left(k\right),{u}_{0}\left(k\right)\right)$ are noiseless.

3) The results when using the LS method to solve linear Equation (27) for noisy data are shown in Table 2(b). All errors became larger when noise was added, suggesting that LS analysis is inadequate when the measurement data are contaminated by noise.

4) The results when using the TLS method to solve linear Equation (27) are shown in Table 2(c). The errors were significantly smaller than those reported in [5] , using the LS method.

5) The results when applying prefiltering (PF) and using the TLS method to

Table 2. Comparison of errors.

Figure 5. Comparison of pole locations (“+” indicates the desired poles, “.” those obtained by the random exciting signal and “o” those obtained by the chirp exciting signal.).

solve linear Equation (27) are shown in Table 2(d). The prefilter further reduced the errors, in particular, the pole placement error $\delta \lambda \left({A}_{\text{d}1}\right)$ and the impulse response error $\Delta {G}_{1}$ .

6) The results when applying PF and using the TLS method to solve the linear Equation (27), but with $v$ as the chirp signal, are shown in Table 2(e). No significant improvement in error rates was found with respect to ${A}_{2}$ when using the chirp exciting signal. However, the errors with respect to ${A}_{1}$ became significantly worse than when a random exciting signal was used. This was assumed to be because ${A}_{1}$ has an unstable eigenvalue of 1.7838. We conclude that a random exciting signal is more appropriate than a chirp exciting signal when using data-driven methods.

Finally, we compare the pole locations obtained as shown in Figure 5. As can be seen, a better performance was achieved when using the random exciting signal.

5. Conclusion

In this study, we evaluated the different approaches reducing the effect of measurement noise in data-driven pole placement methods for deriving a state space model and pole placement state feedback. Using numerical simulations of a self-balancing robot, which is a nonlinear system, we demonstrated the important role that prefiltering can play in reducing the interference caused by noise. Again using numerical simulation, we compared the use of two exciting signals: a random signal and a chirp signal. The use of a random exciting signal was found to be more effective with our proposed method. Further developments are needed in the methods used to cope with noise. A method such as that used in [9] may be appropriate for use in practical applications where noise is present, and adaptive control based on real-time updating [22] is a future promising approach.

Acknowledgements

This work was partially supported by JSPS KAKENHI Grant Number 16H04385.

References

[1] Wonham, W.M. (1967) On Pole Assignment in Multi-input Controllable Linear Systems. IEEE Transaction on Automatic Control, 12, 660-665.

https://doi.org/10.1109/TAC.1967.1098739

[2] Ackermann, J.E. (1977) On the Synthesis of Linear Control Systems with Specified Characteristics. Automatica, 13, 89-94.

[3] Kimura, H. (1975) Pole Assignment by Gain Output Feedback. IEEE Transaction Automatic Control, 20, 509-516.

https://doi.org/10.1109/TAC.1975.1101028

[4] Hikita, H., Koyama, S. and Miura, R. (1975) The Redundancy of Feedback Gain Matrix and the Derivation of Low Feedback Gain Matrix in Pole Assignment. The Society of Instrument and Control Engineers, 11, 556-560. (In Japanese)

[5] Yamamoto, S., Okano, Y. and Kaneko, O. (2016) A Data-driven Pole Placement Method Simultaneously Identifying a State Space Model. Transaction of the Institute of Systems, Control and Information Engineers, 29, 275-284. (In Japanese)

[6] Hou, Z.S. and Wang, Z. (2013) From Model-Based Control to Data-Driven Control: Survey. Classification and Perspective, Information Sciences, 235, 3-35.

[7] Safonov, M.G. and Tsao, T.C. (1997) The Unfalsified Control Concept and Learning. IEEE Transactions on Automatic Control, 42, 843-847.

https://doi.org/10.1109/9.587340

[8] Campi, M.C., Lecchini, A. and Savaresi, S.M. (2002) Virtual Reference Feedback Tuning: A Direct Method for the Design of Feedback Controllers. Automatica, 38, 1337-1346.

[9] Sala, A. and Esparza, A. (2005) Extensions to “Virtual Reference Feedback Tuning: A Direct Method for the Design of Feedback Controllers”. Automatica, 41, 1473-1476.

[10] Souma, S., Kaneko, O. and Fujii, T. (2004) A New Method of a Controller Parameter Tuning Based on Input-output Data-Fictitious Reference Iterative Tuning. Proceedings of the 2nd IFAC Workshop on Adaptation and Learning in Control and Signal Processing, 37, 789-794.

[11] Matsui, Y., Akamatsu, S., Kimura, T., Nakano, K. and Sakurama, K. (2011) Fictitious Reference Iterative Tuning for State Feedback Control of Inverted Pendulum with Inertia Rotor. SICE Annual Conference, Tokyo, 13-18 September 2011, 1087-1092.

[12] Matsui, Y., Akamatsu, S., Kimura, T., Nakano, K. and Sakurama, K. (2014) An Application of Fictitious Reference Iterative Tuning to State Feedback, Electronics and Communications in Japan, 97, 1-11.

https://doi.org/10.1002/ecj.11506

[13] Kaneko, O. (2015) The Canonical Controller Approach to Data-Driven Update of State Feedback Gain. Proceedings of the 10th Asian Control Conference 2015 (ASCC 2015), Kota Kinabalu, 31 May-3 June 2015, 2980-2985.

https://doi.org/10.1109/ASCC.2015.7244744

[14] Van Heusden, K., Karimi, A. and Soderstrom, T. (2011) Extensions to “On Identification Methods for Direct Data-Driven Controller Tuning”. International Journal of Adaptive Control and Signal Processing, 25, 448-465.

https://doi.org/10.1002/acs.1213

[15] Kaneko, O., Miyachi, M. and Fujii, T. (2008) Simultaneous Updating of a Model and a Controller Based on the Data-Driven Fictitious Controller. 47th IEEE Conference on Decision and Control, Cancun, 9-11 December 2008, 1358-1363.

[16] Kaneko, O., Miyachi, M. and Fujii, T. (2011) Simultaneous Updating of Model and Controller Based on Fictitious Reference Iterative Tuning. SICE Journal of Control, Measurement, and System Integration, 4, 63-70.

https://doi.org/10.9746/jcmsi.4.63

[17] Markovsky, I. and Huffel, S.V. (2007) Overview of Total Least-Squares Methods. Signal Processing, 87, 2283-2302.

[18] Shwe, P.E.E. and Yamamoto, S. (2016) Data-Driven Method to Simultaneously Obtain a Linearized State Space Model and Pole Placement Gain. Proceedings of the 3rd Multi Symposium on Control Systems, Nagoya, 7-10 March 2016, 3B3-2.

[19] Brewer, J. (1978) Kronecker Products and Matrix Calculus in System Theory. IEEE Transactions on Circuits and Systems, 25, 772-781.

https://doi.org/10.1109/TCS.1978.1084534

[20] ZMP Inc. Stabilization and Control of Stable Running of the Wheeled Inverted Pendulum, Development of Educational Wheeled Inverted Robot e-nuvo WHEEL ver.1.0.

http://www.zmp.co.jp/e-nuvo/

[21] Nomura, T., Kitsuka, Y., Suemitsu, H. and Matsuo, T. (2009) Adaptive Back Stepping Control for a Two-Wheeled Autonomous Robot. Proceedings of ICCAS-SICE, Fukuoka, 18-21 August 2009, 4687-4692.

[22] Shwe, P.E.E. and Yamamoto, S. (2016) Real-Time Simultaneously Updating a Linearized State-Space Model and Pole Placement Gain. Proceedings of SICE Annual Conference, Tsukuba, 20-23 September 2016, 196-201.