[10.8 MATLAB: TOOLBOXES AND](#page-15-0) STATE-SPACE ANALYSIS

← Back to LINEAR SYSTEMS AND SIGNALS Overview

10.8 MATLAB: TOOLBOXES AND STATE-SPACE ANALYSIS

The preceding MATLAB sections provide a comprehensive introduction to the basic MATLAB environment. However, MATLAB also offers a wide range of toolboxes that perform specialized tasks. Once installed, toolbox functions operate no differently from ordinary MATLAB functions. Although toolboxes are purchased at extra cost, they save time and offer the convenience of predefined functions. It would take significant effort to duplicate a toolbox’s functionality by using custom user-defined programs.

Three toolboxes are particularly appropriate in the study of signals and systems: the control system toolbox, the signal-processing toolbox, and the symbolic math toolbox. Functions from these toolboxes have been utilized throughout the text in the MATLAB examples as well as certain end-of-chapter problems. This section provides a more formal introduction to a selection of functions, both standard and toolbox, that are appropriate for state-space problems.

10.8-1 z-Transform Solutions to Discrete-Time, State-Space Systems

As with continuous-time systems, it is often more convenient to solve discrete-time systems in the transform domain rather than in the time domain. As given in Ex. 10.13, consider the state-space description of the system shown in Fig. 10.12.

\begin{bmatrix} q_1[n+1] \\ q_2[n+1] \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ -\frac{1}{6} & \frac{5}{6} \end{bmatrix} \begin{bmatrix} q_1[n] \\ q_2[n] \end{bmatrix} + \begin{bmatrix} 0 \\ 1 \end{bmatrix} x[n]

and

y[n] = [-1 \quad 5] \begin{bmatrix} q_1[n] \\ q_2[n] \end{bmatrix}

We are interested in the output y[n] in response to the input x[n] = u[n] with initial conditions q1[0] = 2 and q2[0] = 3.

To describe this system, the state matrices A, B, C, and D are first defined.

A = [0 1;-1/6 5/6]; B = [0; 1]; C = [-1 5]; D = 0;

Additionally, the vector of initial conditions is defined.

q_0 = [2;3];

In the transform domain, the solution to the state equation is

\mathbf{Q}[z] = (\mathbf{I} - z^{-1}\mathbf{A})^{-1}\mathbf{q}[0] + (z\mathbf{I} - \mathbf{A})^{-1}\mathbf{B}\mathbf{X}[z] \tag{10.74}

The solution is separated into two parts: the zero-input response and the zero-state response.

MATLAB’s symbolic toolbox makes possible a symbolic representation of Eq. (10.74). First, a symbolic variable z needs to be defined.

z = sym(‘z’);

The sym command is used to construct symbolic variables, objects, and numbers. Typing whos confirms that z is indeed a symbolic object. The syms command is a shorthand command for constructing symbolic objects. For example, syms z s is equivalent to the two instructions z = sym(‘z’); and s = sym(‘s’);.

Next, a symbolic expression for X[z] needs to be constructed for the unit step input, x[n] = u[n]. The z-transform is computed by means of the ztrans command.

X = ztrans(sym(‘1’)) X = z/(z-1)

Several comments are in order. First, the ztrans command assumes a causal signal. For n ≥ 0, u[n] has a constant value of 1. Second, the argument of ztrans needs to be a symbolic expression, even if the expression is a constant. Thus, a symbolic one sym(‘1’) is required. Also note that continuous-time systems use Laplace transforms rather than z-transforms. In such cases, the laplace command replaces the ztrans command.

Construction of Q[z] is now trivial.

\begin{aligned}\n&\text{Q} = \text{inv}\left(\text{eye}(2) - z^(-1) * A\right) * q_0 + \text{inv}\left(z * \text{eye}(2) - A\right) * B * X \\ &\text{Q} = \\ &\left(18 * z\right) / \left(6 * z^2 - 5 * z + 1\right) + \left(2 * z * \left(6 * z - 5\right)\right) / \left(6 * z^2 - 5 * z + 1\right) + \left(6 * z\right) / \left((z - 1) * \left(6 * z^2 - 5 * z + 1\right)\right) \\ &\left(18 * z^2\right) / \left(6 * z^2 - 5 * z + 1\right) - \left(2 * z\right) / \left(6 * z^2 - 5 * z + 1\right) + \left(6 * z^2\right) / \left((z - 1) * \left(6 * z^2 - 5 * z + 1\right)\right)\n\end{aligned}

Unfortunately, not all MATLAB functions work with symbolic objects. Still, the symbolic toolbox overloads many standard MATLAB functions, such as inv, to work with symbolic objects. Recall that overloaded functions have identical names but different behavior; proper function selection is typically determined by context.

The expression Q is somewhat unwieldy. The simplify command uses various algebraic techniques to simplify the result.

>> Q = simplify(Q)
    Q = -(2*z*(- 6*z^2 + 2*z + 1))/(6*z^3 - 11*z^2 + 6*z - 1)
           (2*z*(9*z^2 - 7*z + 1))/(6*z^3 - 11*z^2 + 6*z - 1)

The resulting expression is mathematically equivalent to the original but notationally more compact.

Since D = 0, the output Y[z] is given by Y[z] = CQ[z].

>> Y = simplify(C*Q)
    Y = (6*z*(13*z^2 - 11*z + 2))/(6*z^3 - 11*z^2 + 6*z - 1)

The corresponding time-domain expression is obtained by using the inverse z-transform command iztrans.

y = iztrans(Y) y = 3*(1/2)^n - 2*(1/3)^n + 12

Like ztrans, the iztrans command assumes a causal signal, so the result implies multiplication by a unit step. That is, the system output is y[n] =(3(1/2)ⁿ−2(1/3)ⁿ+12)u[n], which is equivalent to Eq. (10.69) derived in Ex. 10.13. Continuous-time systems use inverse Laplace transforms rather than inverse z-transforms. In such cases, the ilaplace command therefore replaces the iztrans command.

Following a similar procedure, it is a simple matter to compute the zero-input response yzir[n]:

y_zir = iztrans(simplify(C*inv(eye(2)-z^(-1)*A)*q_0)) y_zir = 21*(1/2)^n - 8*(1/3)^n

The zero-state response is given by

y_zsr=y- y_zir y_zsr = 6*(1/3)^n - 18*(1/2)^n + 12

Typing iztrans(simplify(C*inv(z*eye(2)-A)*B*X)) produces the same result.

MATLAB plotting functions, such as plot and stem, do not directly support symbolic expressions. By using the subs command, however, it is easy to replace a symbolic variable with a vector of desired values.

Figure 10.14 Output y[n] computed by using the symbolic math toolbox.

>> n = [0:25]; stem(n,subs(y,n),'k.');
>> xlabel('n'); ylabel('y[n]'); axis([-.5 25.5 11.5 13.5]);

Figure 10.14 shows the results, which are equivalent to the results obtained in Ex. 10.13. Although there are plotting commands in the symbolic math toolbox such as ezplot that plot symbolic expression, these plotting routines lack the flexibility needed to satisfactorily plot discrete-time functions.

10.8-2 Transfer Functions from State-Space Representations

A system’s transfer function provides a wealth of useful information. From Eq. (10.73), the transfer function for the system described in Ex. 10.13 is

>> H = collect(simplify(C*inv(z*eye(2)-A)*B+D))
    H = (30*z - 6)/(6*z^2 - 5*z + 1)

It is also possible to determine the numerator and denominator transfer function coefficients from a state-space model by using the signal-processing toolbox function ss2tf.

[num,den] = ss2tf(A,B,C,D) num = 0 5.0000 -1.0000 den = 1.0000 -0.8333 0.1667

The denominator of H[z] provides the characteristic polynomial

^γ ² ⁻ ⁵ ⁶^γ ⁺ ¹ 6

Equivalently, the characteristic polynomial is the determinant of (zI−A).

syms gamma; char_poly = subs(det(z*eye(2)-A),z,gamma) char_poly = gamma^2 - (5*gamma)/6 + 1/6

Here, the subs command replaces the symbolic variable z with the desired symbolic variable gamma.

The roots command does not accommodate symbolic expressions. Thus, the sym2poly command converts the symbolic expression into a polynomial coefficient vector suitable for the roots command.

roots(sym2poly(char_poly)) ans = 0.5000 0.3333

Taking the inverse z-transform of H[z] yields the impulse response h[n].

h = iztrans(H) h = 18*(1/2)^n - 12*(1/3)^n - 6*kroneckerDelta(n, 0)

As suggested by the characteristic roots, the characteristic modes of the system are (1/2)ⁿ and (1/3)ⁿ. Notice that the symbolic math toolbox represents δ[n] as kroneckerDelta(n, 0). In general, δ[n − a] is represented as kroneckerDelta(n-a, 0). This notation is frequently encountered. Consider, for example, delaying the input x[n] = u[n] by 2, x[n − 2] = u[n − 2]. In the transform domain, this is equivalent to z−2X[z]. Taking the inverse z-transform of z−2X[z] yields

\Rightarrow \text{ 27 (--2)*X} = \text{ 27 (--2)*X}

\text{ans} = 1 - \text{ 27 (--2)*X} = \text{ 27 (--2)*X}

That is, MATLAB represents the delayed unit step u[n−2] as (−δ[n−1] −δ[n−0] +1)u[n]. The transfer function also permits convenient calculation of the zero-state response.

y\_zsr = iztrans(H*X)

y_ $zsr = 6*(1/3)^n - 18*(1/2)^n + 12$

The result agrees with previous calculations.

10.8-3 Controllability and Observability of Discrete-Time Systems

In their controllability and observability, discrete-time systems are analogous to continuous-time systems. For example, consider the LTID system described by the constant coefficient difference equation

y[n] + \frac{5}{6}y[n-1] + \frac{1}{6}y[n-2] = x[n] + \frac{1}{2}x[n-1]

Figure 10.15 illustrates the direct form II (DFII) realization of this system. The system input is x[n], the system output is y[n], and the outputs of the delay blocks are designated as state variables q1[n] and q2[n].

The corresponding state and output equations (see Prob. 10.7-1) are

\mathbf{Q}[n+1] = \begin{bmatrix} q_1[n+1] \\ q_2[n+1] \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ -\frac{1}{6} & -\frac{5}{6} \end{bmatrix} \begin{bmatrix} q_1[n] \\ q_2[n] \end{bmatrix} + \begin{bmatrix} 0 \\ 1 \end{bmatrix} x[n] = \mathbf{A}\mathbf{Q}[n] + \mathbf{B}x[n]

and

y[n] = \begin{bmatrix} -\frac{1}{6} & -\frac{1}{3} \end{bmatrix} \begin{bmatrix} q_1[n] \\ q_2[n] \end{bmatrix} + 1x[n] = \mathbf{CQ}[n] + \mathbf{D}x[n]

To describe this system in MATLAB, the state matrices A, B, C, and D are first defined.

A =

[0 \ 1; -1/6 \ -5/6]

; B = $[0; 1]$ ; C = $[-1/6 \ -1/3]$ ; D = 1;

Figure 10.15 Direct form II realization of y[n] + (5/6)y[n − 1] + (1/6)y[n − 2] = x[n] + (1/2)x[n−1].

To assess the controllability and observability of this system, the state matrix A needs to be diagonalized.† As shown in Eq. (10.55), this requires a transformation matrix P such that

PA = AP

(10.75)

where is a diagonal matrix containing the unique eigenvalues of A. Recall, the transformation matrix P is not unique.

To determine a matrix P, it is helpful to review the eigenvalue problem. Mathematically, an eigendecomposition of A is expressed as

AV = V\Lambda

where V is a matrix of eigenvectors and is a diagonal matrix of eigenvalues. Pre- and post-multiplying both sides of this equation by V−¹ yields

\mathbf{V}^{-1}\mathbf{A}\mathbf{V}\mathbf{V}^{-1} = \mathbf{V}^{-1}\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{-1}

Simplification yields

\mathbf{V}^{-1}\mathbf{A} = \mathbf{\Lambda}\mathbf{V}^{-1} \tag{10.76}

Comparing Eqs. (10.75) and (10.76), we see that a suitable transformation matrix P is given by an inverse eigenvector matrix V−¹ .

The eig command is used to verify that A has the required distinct eigenvalues as well as compute the needed eigenvector matrix V.

>> [V,Lambda] = eig(A)
   V =
       0.9487 -0.8944
      -0.3162 0.4472
   Lambda =
      -0.3333 0
           0 -0.5000

Since the diagonal elements of Lambda are all unique, a transformation matrix P is given by

P = inv(V);

The transformed state matrices ^A^ˆ ⁼ PAP−¹ , ^B^ˆ ⁼ PB, and ^C^ˆ ⁼ CP−¹ are easily computed by using transformation matrix P. Notice that matrix D is unaffected by state variable transformations.

>> Ahat = P*A*inv(P), Bhat = P*B, Chat = C*inv(P)
   Ahat =
      -0.3333 -0.0000
       0.0000 -0.5000

^† This approach requires that the state matrix A have unique eigenvalues. Systems with repeated roots require that state matrix A be transformed into a modified diagonal form, also called the Jordan form. The MATLAB function jordan is used in these cases.

Bhat =
    6.3246
    6.7082
Chat =
   -0.0527 -0.0000

The proper operation of P is verified by the correct diagonalization of A, Aˆ = . Since no row of Bˆ is zero, the system is controllable. Since, however, at least one column of Cˆ is zero, the system is not observable. These characteristics are no coincidence. The DFII realization, which is more descriptively called the controller canonical form, is always controllable but not always observable.

As a second example, consider the same system realized using the transposed direct form II structure (TDFII), as shown in Fig. 10.16. The system input is x[n], the system output is y[n], and the outputs of the delay blocks are designated as state variables v1[n] and v2[n].

The corresponding state and output equations (see Prob. 10.7-2) are

\mathbf{V}[n+1] = \begin{bmatrix} v_1[n+1] \\ v_2[n+1] \end{bmatrix} = \begin{bmatrix} 0 & -\frac{1}{6} \\ 1 & -\frac{5}{6} \end{bmatrix} \begin{bmatrix} v_1[n] \\ v_2[n] \end{bmatrix} + \begin{bmatrix} -\frac{1}{6} \\ -\frac{1}{3} \end{bmatrix} x[n] = \mathbf{A}\mathbf{V}[n] + \mathbf{B}x[n]

and

y[n] = [0 \quad 1] \begin{bmatrix} v_1[n] \\ v_2[n] \end{bmatrix} + 1x[n] = \mathbf{CV}[n] + \mathbf{D}x[n]

To describe this system in MATLAB, the state matrices A, B, C, and D are defined.

A =

[0 -1/6; 1 -5/6]

; B = $[-1/6; -1/3]$ ; C = $[0 1]$ ; D = 1;

To diagonalize A, a transformation matrix P is created.

[V, Lambda] = eig(A)
\n

V =

0.4472 0.3162
\n0.8944 0.9487
\nLambda =
\n-0.3333 0
\n0 -0.5000

⁶ – Figure 10.16 Transposed direct form II realization of y[n]+(5/6)y[n−1]+(1/6)y[n−2] = x[n]+(1/2)x[n−1].

The characteristic modes of a system do not depend on implementation, so the eigenvalues of the DFII and TDFII realizations are the same. However, the eigenvectors of the two realizations are quite different. Since the transformation matrix P depends on the eigenvectors, different realizations can possess different observability and controllability characteristics.

Using transformation matrix ^P, the transformed state matrices ^A^ˆ ⁼ PAP−¹ , Bˆ = PB, and ^C^ˆ ⁼ CP−¹ are computed.

>> P = inv(V);
>> Ahat = P*A*inv(P), Bhat = P*B, Chat = C*inv(P)
   Ahat =
      -0.3333 0
       0.0000 -0.5000
   Bhat =
      -0.3727
      -0.0000
   Chat =
       0.8944 0.9487

Again, the proper operation of P is verified by the correct diagonalization of A, Aˆ = . Since no column of Cˆ is zero, the system is observable. However, at least one row of Bˆ is zero, and therefore the system is not controllable. The TDFII realization, which is more descriptively called the observer canonical form, is always observable but not always controllable. It is interesting to note that the properties of controllability and observability are influenced by the particular realization of a system.

10.8-4 Matrix Exponentiation and the Matrix Exponential

Matrix exponentiation is important to many problems, including the solution of discrete-time state-space equations. Equation (10.64), for example, shows that the state response requires matrix exponentiation, A*ⁿ* . For a square A and specific n, MATLAB happily returns A*ⁿ* by using the ^ operator. From the system in Ex. 10.13 and n = 3, we have

>> A = [0 1;-1/6 5/6]; n = 3; A^n
   ans =
      -0.1389 0.5278
      -0.0880 0.3009

The same result is also obtained by typing A*A*A.

Often, it is useful to solve ^Aⁿ symbolically. Noting ^Aⁿ ⁼ ^Z⁻¹[(^I ⁻ ^z−¹A)−¹], the symbolic toolbox can produce a symbolic expression for A*ⁿ* .

>> syms z n; An = simplify(iztrans(inv(eye(2)-z^(-1)*A)))
   An =
   [ 3*(1/3)^n - 2*(1/2)^n, 6*(1/2)^n - 6*(1/3)^n]
   [ 1/3^n - 1/2^n, 3/2^n - 2/3^n]

Notice that this result is identical to Eq. (10.68), derived earlier. Substituting the case n = 3 into An provides a result that is identical to the one elicited by the previous A^n command.

>> double(subs(An,n,3)) ans = -0.1389 0.5278 -0.0880 0.3009

For continuous-time systems, the matrix exponential eA*^t* is commonly encountered. The expm command can compute the matrix exponential symbolically. Using the system from Ex. 10.8 yields

>> syms t; A = [-12 2/3;-36 -1]; eAt = simplify(expm(A*t))
    eAt =
    [ -(exp(-9*t)*(3*exp(5*t) - 8))/5, (2*exp(-9*t)*(exp(5*t) - 1))/15]
    [ -(36*exp(-9*t)*(exp(5*t) - 1))/5, (exp(-9*t)*(8*exp(5*t) - 3))/5]

This result is identical to the result computed in Ex. 10.8. Similar to the discrete-time case, an identical result is obtained by typing syms s; simplify(ilaplace(inv(s*eye(2)-A))).

For a specific t, the matrix exponential is also easy to compute, either through substitution or direct computation. Consider the case t = 3.

double(subs(eAt,t,3)) ans = 1.0e-004 * -0.0369 0.0082 -0.4424 0.0983

The command expm(A*3) produces the same result.