Problem in verifying zero blocking property of a zero of a multiple input multiple output system

Question

I have a multiple input multiple output control system, and in particular a two inputs two outputs system. I want to find the zero state direction and the zero input direction of the system, which has a zero at $s=0$. So, I have foud them with the following code:

s = tf('s');

P = 1/(s+5);
C = 6/s;

S = 1/(1+P*C);
T = P*C/(1+P*C);

G_2 = [T S; S -S]; %transfer matrix

S = ss(G_2)
Smr = minreal(S)
A = Smr.A
B = Smr.B
C = Smr.C
D = Smr.D
trz = tzero(A,B,C,D)
RSM_1 = [0*eye(2)-A B; -C D]
rRSM_1 = rank(RSM_1)

null = null(RSM_1)


%verifying zero blocking property
x0 = [-0.4388; -0.1482];
u0 = [-0.0000; 0.8863];
t = linspace(0,5);
u = exp(t);
[y,x] = lsim(A,B*u0,C,D*u0,u,t,x0);
plot(t,y)
axis([0 5 -1 1])
xlabel('time, seconds')
title('response to x_{0} and u(t)=u_{0}e^{t}')

and I want to verify the zero blocking property, which means that if I apply an input at $u(t)=u_0e^{t}$ starting from $x(0)=x_0$, where $u_0$ is the input direction and $x_0$ is the state direction, the output is zero.

I have found it analytically in this way:

I know that the since the system zeros of a MIMO transfer function are those values $s = z$ which make the system matrix $S(z)$ loose rank, there exists at least a direction $u_0$ (zero input direction) and a n-dimensional vector $x_0$ (zero state direction), not simultaneously zero, such that

$\begin{pmatrix} zI-A & -B\\ C&D \end{pmatrix}$ $\begin{bmatrix} x_0\\ u_0 \end{bmatrix}$ $=$ $\begin{bmatrix} 0\\ 0 \end{bmatrix}$

And by solving this system we obtain the zero state direction $x_0$ and the zero input direction $u_0$.

In order to make the question even clearer I can make an example using a different transfer function which has a simpler minimim realization. The example is taken from this : Here, and it can be found at the bottom of the paper.

The system in this case is:

$P(s)=\begin{bmatrix} \frac{2}{s^2+3s+2} & \frac{2s}{s^2+3s+2}\\ \frac{-2s}{s^2+3s+2}& \frac{-2}{s^2+3s+2} \end{bmatrix}$

Its minimum realization is:

$A=\begin{bmatrix} -1 &0 &0 \\ 0&-2 &0 \\ 0& 0 & -2 \end{bmatrix}$

$B=\begin{bmatrix} 2 &-2 \\ -2& 4\\ -4& 2 \end{bmatrix}$

$C=\begin{bmatrix} 1 &1 &0 \\ 1 &0 &1 \end{bmatrix}$

$D=\begin{bmatrix} 0 &0 \\ 0& 0 \end{bmatrix}$

and the after using the same code, which can also be found in the paper, the simulation for verifying the zero blocking property is:

s = tf('s');
P = [2/(s^2+3*s+2) 2*s/(s^2+3*s+2); -2*s/(s^2+3*s+2) -2/(s^2+3*s+2)]
S = ss(P)
Smr = minreal(S)
A = Smr.A
B = Smr.B
C = Smr.C
D = Smr.D
trz = tzero(A,B,C,D)
RSM_1 = [eye(3)-A B; -C D]
rRSM_1 = rank(RSM_1)

null(RSM_1)

x0 = [0.5345; -0.5345; -0.5345];
u0 = [0.2673; -0.2673];
t = linspace(0,5);
u = exp(t);
[y,x] = lsim(A,B*u0,C,D*u0,u,t,x0);
plot(t,y),grid
axis([0 5 -1 1])
xlabel('time, seconds')
title('response to x_{0} and u(t)=u_{0}e^{t}')

and the output of this is:

which is correct, since the output coes to zero.

But, what I obtain with my code is:

I am still working on this, and by using the Matlab command isstable(T) it tells me that the closed loop is not stable. But if I run the step response, step(T) ,I find:

what am I doing wrong?

Answer 1

The transmission zeros of a MIMO system are obtained by solving

\begin{align} \dot{x}^* = A\,x^* + B\,u^*, \\ y = C\,x^* + D\,u^*, \end{align}

for $x^*$ and $u^*$ such that $y = 0$. By using the Laplace variable $s$ for the time derivative that system of equation can also be written as

$$ \begin{bmatrix} 0 \\ 0 \end{bmatrix} = \begin{bmatrix} A - s\,I & B \\ C & D \end{bmatrix} \begin{bmatrix} x^* \\ u^* \end{bmatrix}, $$

which can be seen as the generalized eigenvalue problem $\mathcal{A}\,v = \lambda\,\mathcal{B}\,v$, with $\lambda$ equivalent to $s$,

$$ \mathcal{A} = \begin{bmatrix} A & B \\ C & D \end{bmatrix}, \quad \mathcal{B} = \begin{bmatrix} I & 0 \\ 0 & 0 \end{bmatrix}, \quad v = \begin{bmatrix} x^* \\ u^* \end{bmatrix}. $$

I believe that tzero in matlab also solves such generalized eigenvalue problem to obtain the transmission zeros (however by using eig yourself the associated eigenvectors are more accurate than the nullspace calculations). The associated solution in time for both $x(t)$ and $u(t)$ can be obtained using

$$ \begin{bmatrix} x(t) \\ u(t) \end{bmatrix} = e^{\lambda\,t}\, \begin{bmatrix} x^* \\ u^* \end{bmatrix}, $$

so you can't just pick any value for $\lambda$ when generating such solution. You could now simulate the system using $x(0)=x^*$ as initial condition and the expression for $u(t) = e^{\lambda\,t}u^*$ with lsim. However, I believe that slim uses zero-order-hold for the input so the actual applied input will deviate slightly causing the output $y$ to deviate away from zero using that approach. Instead, if you look at the analytical solution, both the output and the error in state dynamics should remain practically to zero (for non-positive $\lambda$):

n = length(A);
m = length(D);
[V,L] = eig([A B; C D], diag([ones(1,n) zeros(1,m)]));
ind = 1; % Index of associated eigenvalue (avoid infinite eigenvalues)
lambda = L(ind,ind);
x0 = V(1:n, ind);
u0 = V(1+n:n+m, ind);

N = 1e3;
t = linspace(0, 5, N);
x = @(t) x0 * exp(lambda* t);
u = @(t) u0 * exp(lambda* t);
y = @(t) C * x(t) + D * u(t);
e = @(t) A * x(t) + B * u(t) - lambda* x(t);

Y = zeros(n, N);
E = zeros(m, N);
for k = 1 : N
Y(:,k) = y(t(k));
E(:,k) = e(t(k));
end
figure, plot(t, E)
figure, plot(t, Y)