Problems

a. Eigendecomposition

The eigendecomposition is performed using the eig function.

[M,L] = eig(A)

M =

   0.6667 + 0.0000i   0.6667 + 0.0000i   0.8944 + 0.0000i
  -0.0000 + 0.6667i  -0.0000 - 0.6667i   0.0000 + 0.0000i
  -0.3333 - 0.0000i  -0.3333 + 0.0000i   0.4472 + 0.0000i


L =

  -2.0000 + 2.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i  -2.0000 - 2.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i  -2.0000 + 0.0000i

Let’s take a closer look at the eigenvalues.

disp(diag(L))

  -2.0000 + 2.0000i
  -2.0000 - 2.0000i
  -2.0000 + 0.0000i

b. Stability and character of response

All three eigenvalues have negative real parts; therefore, the system is stable. There is a single pair of complex conjugate eigenvalues, which will cause the system to oscillate while decaying. The real eigenvalue will contribute a real exponential decay with the same rate as the oscillating pair.

c. State transition matrices

First, we construct \(\Phi'\) symbolically.

So we can find this from the state transition matrix \(\Phi\), which is known from to be \(M \Phi'(t) M^{-1}\).

syms t
Phi_p = expm(L*t)

Phi_p =
 
[ exp(t*(- 2 + 2i)),                 0,         0]
[                 0, exp(t*(- 2 - 2i)),         0]
[                 0,                 0, exp(-2*t)]
 

Now we can apply our transformation.

Phi = M*Phi_p*inv(M); % messy to print

d. Free response

Let’s compute the free response to some initial conditions. The free state response is given by \[\begin{aligned} \bm{x}(t) = \Phi(t) \bm{x}(0). \end{aligned}\]

So our symbolic solution is to multiply the initial conditions by this matrix.

x_0 = [0;0;1]; % initial conditions
x = Phi*x_0; % free response, symbolically

e. Plotting a free response

Let’s make the symbolic solution into something we can evaluate numerically and plot.

x_fun = matlabFunction(x);

Now let’s set up our time array for the plot.

t_a = linspace(0,4,300);

figure
plot(t_a,real(x_fun(t_a)))
xlabel('time (s)')
ylabel('state free response')
legend('x_1','x_2','x_3')

a and b Eigenvalue and modal matrices

Use eig for the eigendecomposition. The problem asks for the eigenvalues in the order opposite that of Matlab’s default eig order, so we flip them in L and M.

[M,L] = eig(A)

M =
            1            0     -0.97014
            0            1            0
            0            0      0.24254
L =
    -1     0     0
     0    -2     0
     0     0    -3

c State transition matrix

Proceed symbolically with the definition of the diagonalized basis state-transition matrix \(\Phi'(t)\).

syms t
Phi_p = @(t) diag(diag(exp(L*t)));
Phi_p(t)

ans =
[ exp(-t),         0,         0]
[       0, exp(-2*t),         0]
[       0,         0, exp(-3*t)]

Now transform bases to the original basis using \[\begin{aligned} \Phi(t) = M \Phi'(t) M^{-1}. \end{aligned}\]

M_inv = M^-1;
Phi = @(t) M*Phi_p(t)*M_inv;
Phi(t)

ans =
[ exp(-t),         0, 4*exp(-t) - 4*exp(-3*t)]
[       0, exp(-2*t),                       0]
[       0,         0,               exp(-3*t)]

d Forced state response

The forced state response is the second term of the state response solution.

syms T
u = @(t) [4;sin(2*pi*t)]; % input
x = @(t) int(Phi(t - T)*B*u(T),T,0,t); % integrate!
pretty(x(t))

/ (sin(2 pi t) - pi cos(2 pi t) 2) 2   4 pi exp(-t) \
| ---------------------------------- + ------------ |
|                  2                         2      |
|              4 pi  + 1                 4 pi  + 1  |
|                                                   |
|                  6 - exp(-2 t) 6                  |
|                                                   |
\                         0                         /

Written out somewhat more accessibly, \[\begin{aligned} \bm{x}(t) &= \left(\begin{array}{c} \frac{2\,\left(\sin\left(2\,\pi \,t\right)-2\,\pi \,\cos\left(2\,\pi \,t\right)\right)}{4\,\pi ^2+1}+\frac{4\,\pi \,{\mathrm{e}}^{-t}}{4\,\pi ^2+1}\\ 6-6\,{\mathrm{e}}^{-2\,t}\\ 0 \end{array}\right). \end{aligned}\]

e Forced output response

The state forced response is known, so the output response is just \[\begin{aligned} \bm{y}(t) &= C \bm{x} + D \bm{u}. \end{aligned}\]

y = @(t) C*x(t) + D*u(t);
y(t)

ans =
(2*(sin(2*pi*t) - 2*pi*cos(2*pi*t)))/(4*pi^2 + 1) + (4*pi*exp(-t))/(4*pi^2 + 1)

More readably, \[\begin{aligned} \bm{y}(t) &= \frac{2\,\left(\sin\left(2\,\pi \,t\right)-2\,\pi \,\cos\left(2\,\pi \,t\right)\right)}{4\,\pi ^2+1}+\frac{4\,\pi \,{\mathrm{e}}^{-t}}{4\,\pi ^2+1}. \end{aligned}\]

We might consider combining the \(\sin\) and \(\cos\) terms, which are truly a single sinusoid with amplitude and phase depending on their coefficients.

f Plotting the forced output response

Plot the forced output response, shown in .

y_fun = matlabFunction(y(t));
t_a = linspace(0,7,200);
figure;
plot(t_a,y_fun(t_a),'linewidth',1.5);
xlabel('time (s)')
ylabel('output forced response')
grid on

a and b Eigenvalue and modal matrices

Use eig for the eigendecomposition. The problem asks for the eigenvalues in the order opposite that of Matlab’s default eig order, so we flip them in L and M.

[M,L] = eig(A);

disp('The eigenvalue matrix Lambda is:')
disp(L)
disp('The modal matrix M is:')
disp(M)

The eigenvalue matrix Lambda is:
    -4     0
     0   -11
The modal matrix M is:
      0.98639     -0.70711
       0.1644      0.70711

c State transition matrix

Proceed symbolically with the definition of the diagonalized basis state-transition matrix \(\Phi'(t)\).

syms t
Phip = @(t) diag(exp(diag(L)*t))

Now transform bases to the original basis using \[\begin{aligned} \Phi(t) = M \Phi'(t) M^{-1}. \end{aligned}\]

Phi = @(t) M*Phip(t)*M^-1

d Forced state response

The forced state response is the second term of the state response solution.

syms tau

u = @(t) [dirac(t)];
sympref('HeavisideAtOrigin',0); % consistent with RW

xfo = int(Phi(t-tau)*B*u(tau),tau,0,t);
assume(t > 0) % technically t >= 0, so we should check that x(0)=0
xfo = simplify(xfo,'steps',10);

Written out somewhat more accessibly, \[\begin{aligned} \bm{x}(t) &= \left(\begin{array}{c} \frac{2\,\left(\sin\left(2\,\pi \,t\right)-2\,\pi \,\cos\left(2\,\pi \,t\right)\right)}{4\,\pi ^2+1}+\frac{4\,\pi \,{\mathrm{e}}^{-t}}{4\,\pi ^2+1}\\ 6-6\,{\mathrm{e}}^{-2\,t}\\ 0 \end{array}\right). \end{aligned}\]

e Forced output response

The state forced response is known, so the output response is just \[\begin{aligned} \bm{y}(t) &= C \bm{x} + D \bm{u}. \end{aligned}\]

y_fo = C*xfo + D*u(t)

\(y_fo =\frac{2\,{\mathrm{e}}^{-11\,t} \,{\left({\mathrm{e}}^{7\,t} -1\right)}}{7}\)

f Plotting the forced output response

Plot the forced output response, shown in .

yfo_fun = matlabFunction(yfo);
tN = linspace(0,3,200);
yfoN = yfo_fun(tN);
figure
plot(tN,yfoN,'linewidth',2);
xlabel('time (s)')
ylabel('forced response y_{fo}(t)')
grid on

Define the System

Define a symbolic state-space matrix using the dysys package as follows:

A = [[-4, -3, 0], [0, -8, 4], [0, 0, -1]]
B = [[0], [1], [0]]
C = [[0, 1, 0]]
D = [[0]]
sys = dysys.sss(A, B, C, D)  # Create a symbolic state-space model

Eigenvalues and Stability

The eigenvalue matrix L can be found via the eig() method:

L, M = sys.eig()
print(f"Eigenvalues: {L.diagonal().tolist()}")

Eigenvalues: [[-8, -4, -1]]

The real parts of the eigenvalues are all negative; therefore, the system is asymptotically stable.

Eigenvectors and the Modal Matrix

The eigenvectors are stored in M:

print(M)

Matrix([[3/4, 1, -4/7], [1, 0, 4/7], [0, 0, 1]])

State Transition Matrix

The state transition matrix \(\Phi(t)\) can be found as follows:

t = sp.Symbol("t", nonnegative=True)  ### Solution valid for t >= 0
Phi = sys.state_transition_matrix(t)
pprint(Phi)

Matrix([
[exp(-4*t), -3*exp(-4*t)/4 + 3*exp(-8*t)/4, -4*exp(-t)/7 + exp(-4*t) - 3*exp(-8*t)/7],
[        0,                      exp(-8*t),              4*exp(-t)/7 - 4*exp(-8*t)/7],
[        0,                              0,                                  exp(-t)]])

Forced State Response

The forced state response for \(t \ge 0\) can be found as follows:

u_s = sp.Heaviside(t)
x_fo = sys.state_forced_response(t, u=u_s)
pprint(x_fo)

Matrix([
[-3/32 + 3*exp(-4*t)/16 - 3*exp(-8*t)/32],
[                      1/8 - exp(-8*t)/8],
[                                      0]])

Forced Output Response

y_fo = sys.output_forced_response(t, u=u_s)
pprint(y_fo)

Matrix([[1/8 - exp(-8*t)/8]])

Plot the Forced Output Response

First, convert the symbolic expression to a NumPy function:

y_fo_fun = sp.lambdify(t, y_fo, "numpy")

Now create arrays to plot:

t_ = np.linspace(0, 5, 101)
y_fo_ = y_fo_fun(t_).flatten()

Finally, plot:

fig, ax = plt.subplots()
ax.plot(t_, y_fo_)
ax.set_xlabel("Time (s)")
ax.set_ylabel("Forced Output Response $y_\\text{fo}$")
plt.show()