## Chapter 10


```python
import sympy
import numpy as np
sympy.init_printing(use_latex='mathjax')
s, w = sympy.symbols('s w', real=True)
import control
import matplotlib.pyplot as plt
from scipy.optimize import minimize
```

### Problem 10.1
####  a)
##### i)
The spectral density $S(\omega)$ is:
\begin{eqnarray}
S(\omega) &=& H(-j \omega) H^T(j \omega) \\
&=& \frac{1}{(\omega^2 + \alpha^2)(\omega^2 + \beta^2)} \\
&=& \frac{1}{\beta^2 - \alpha^2} \biggl [ \frac{1}{\omega^2 + \alpha^2} - \frac{1}{\omega^2 + \beta^2} \biggr ]
\end{eqnarray}
where the last result comes from partial fraction decomposition.

##### ii)
\begin{eqnarray}
\lim_{T \to \infty} \int_{-T/2}^{T/2} x^2 (\tau) d \tau &=& \rho(0) \\
&\triangleq& \frac{1}{2 \pi} \int_{-\infty}^\infty S(\omega) d \omega \\
&=& \frac{1}{2 \pi} \frac{1}{\beta^2 - \alpha^2} \int_{-\infty}^\infty \biggl [ \frac{1}{\omega^2 + \alpha^2} - \frac{1}{\omega^2 + \beta^2} \biggr ] \\
&=& \frac{1}{2 \pi} \frac{1}{\beta^2 - \alpha^2} \biggl [ \frac{\pi}{\alpha} - \frac{\pi}{\beta} \biggr ] \\
&=& \frac{1}{2} \frac{1}{(\beta + \alpha)(\beta - \alpha)} \biggl [ \frac{\beta - \alpha}{\alpha \beta} \biggr ] \\
&=& \frac{1}{2} \frac{1}{\alpha \beta(\beta + \alpha)}
\end{eqnarray}

##### iii)
Put the transfer function into the "first-companion" form:


```python
p1, p2, p3, a, b = sympy.symbols('p1 p2 p3 a b')
# use first companion form (Eqn. 3.88)
A = sympy.Matrix([[0, 1], [-a*b, -(a+b)]])
F = sympy.Matrix([[0], [1]])
Q = 1
P = sympy.Matrix([[p1, p2], [p2, p3]])
Z = A*P + P*A.T + F*Q*F.T
# Z == [[0, 0], [0, 0]], so 3 equations and 3 unknowns
sympy.pprint('{} = 0'.format(Z[0, 0]))
sympy.pprint('{} = 0'.format(Z[0, 1]))
sympy.pprint('{} = 0'.format(Z[1, 1]))
```

    2*p2 = 0
    -a*b*p1 + p2*(-a - b) + p3 = 0
    -2*a*b*p2 + 2*p3*(-a - b) + 1 = 0


The above linear equation leads to the covariance matrix:
\begin{eqnarray}
\bar P = \begin{pmatrix}
\frac{1}{2 \alpha \beta (\alpha+\beta)} & 0 \\
0 & \frac{1}{2 \alpha \beta}
\end{pmatrix},
\end{eqnarray}
and, for the mean-squared value, $\rho(0) = C^T \bar P e^{A \cdot 0} C = C^T \bar P C$.


```python
P = sympy.Matrix([[1/(2*a*b*(a+b)), 0], [0, 1/(2*a*b)]])
C = sympy.Matrix([[1], [0]])  # from the first companion form definition
qdrtc = sympy.simplify((C.T*P*C)[0])
sympy.pprint(qdrtc)
```

          1      
    ─────────────
    2⋅a⋅b⋅(a + b)


####  b)
_omitted because it's basically the same problem as above_

### Problem 10.2
_this is how I would solve the problem, but the kernel keeps crashing so I'm moving on_


```python
# data to use for regression
frequencies = np.array([0.1, 0.35, 1., 5., 100.])
Sw = np.array([90., 19000., 10., 1.])
# objective function - solve an optimization problem for the coefficients in feature
def obj(feature, freqs, powSpecDens):
    K, a, icsi1, icsi2, icsi3, w1, w2, w3 = feature
    def S(w):
        num = K**2 * w**4 * ( (w3-w)**2 + 4.*icsi3**2*w3**2*w**2 )
        den = ( (w**2 + a**2) * ( (w1-w)**2 + 4.*icsi1**2*w1**2*w**2 )*( (w2-w)**2 + 4.*icsi2**2*w2**2*w**2 ) )
        return num / den
    return np.sum([(S(f) - psd)**2 for f,psd in zip(frequencies, powSpecDens)])
# initial guess
feature = np.array([100., 1., 100., 100., 100., 0.5, 0.5, 0.5])
# bounds
bnds = ((0, 1000),(0, 10),(0, 1000),(0, 1000),(0, 1000),(0, 1),(0, 1),(0, 1))
# solve optimization problem
sol = minimize(obj, feature, args=(frequencies, Sw), bounds=bnds)
print(sol)
```

### Problem 10.3
#### a)
The covariance equation can be written as:
\begin{eqnarray}
\dot P = AP + PA^T + FQ_vF^T,
\end{eqnarray}
where $Q_v$ is the diagonal matrix representing the expected value of the white-noise and
\begin{eqnarray}
A = \begin{pmatrix}
0 & 1 & 0 \\
0 & 0 & -g \\
0 & \frac{1}{R} & 0
\end{pmatrix}, F = \begin{pmatrix}
0 & 0 \\
1 & 0 \\
0 & 1
\end{pmatrix}.
\end{eqnarray}

#### b)
_There are no numbers for the Q matrix, so just setup the problem_


```python
# numbers from Example 3E
Omega = 1.235e-3
g = 9.8
R = g / Omega**2
A = np.array([[0, 1, 0], [0, 0, -g], [0, 1./R, 0]]).astype('float')
F = np.array([[0, 0], [1, 0], [0, 1]]).astype('float')
p1d, p2d, p3d, p4d, p5d, p6d = sympy.symbols('p1d p2d p3d p4d p5d p6d')
p1, p2, p3, p4, p5, p6 = sympy.symbols('p1 p2 p3 p4 p5 p6')
lhs = sympy.Matrix([[p1d, p2d, p3d], [p2d, p4d, p5d], [p3d, p5d, p6d]])
P = sympy.Matrix([[p1, p2, p3], [p2, p4, p5], [p3, p5, p6]])
w1, w2 = sympy.symbols('w1 w2')
Q = sympy.Matrix([[w1, 0], [0, w2]])
rhs = sympy.Matrix(A)*P + P*sympy.Matrix(A.T) + sympy.Matrix(F)*Q*sympy.Matrix(F.transpose())
print('{} = {}'.format(lhs, rhs))
```

    Matrix([[p1d, p2d, p3d], [p2d, p4d, p5d], [p3d, p5d, p6d]]) = Matrix([[2.0*p2, -9.8*p3 + 1.0*p4, 1.55635204081633e-7*p2 + 1.0*p5], [-9.8*p3 + 1.0*p4, -19.6*p5 + 1.0*w1, 1.55635204081633e-7*p4 - 9.8*p6], [1.55635204081633e-7*p2 + 1.0*p5, 1.55635204081633e-7*p4 - 9.8*p6, 3.11270408163265e-7*p5 + 1.0*w2]])


57ac95e9b99db6b24af8fabdf641cd20e9bc1888
ipynb
Jupyter Notebook
Ch10/FriedlandCh10.ipynb
jwdinius/friedland-csd-solutions
2d7fec1b6055ea9777dc931bc1cf0f0df6c8ed7b
MIT
2019-03-24T12:06:18.000Z
2021-09-21T17:56:51.000Z
Ch10/FriedlandCh10.ipynb
jwdinius/friedland-csd-solutions
2d7fec1b6055ea9777dc931bc1cf0f0df6c8ed7b
MIT
2021-03-31T04:19:23.000Z
2022-03-08T21:16:28.000Z
Ch10/FriedlandCh10.ipynb
jwdinius/friedland-csd-solutions
2d7fec1b6055ea9777dc931bc1cf0f0df6c8ed7b
MIT
2019-10-05T04:23:29.000Z
2022-03-08T21:17:17.000Z
Qwen/Qwen-72B
1. YES
2. YES