Nonparametric Identification

1. Parametric vs Nonparametric Models

a) Parametric Models

• Definition: Models described by a finite set of parameters, e.g., transfer functions or state-space models.
• Examples:
  • Transfer function: $$G(s) = \frac{K}{(s+1)(s+3)(s+5)}$$
  • Discrete-time transfer function: $$G(z) = \frac{z^2 + 2z + 3}{z^3 + 3z^2 + 3z + 8}$$
• Characteristics:
  • Suitable for controller design, simulation, and prediction.
  • Different representations can be computed easily.
  • Require a priori assumptions about model structure.
  • Parameters are estimated via optimization.
  • Measurement noise can be filtered out during optimization.

b) Nonparametric Models

• Definition: Models based on data or graphical representations without assuming a fixed parametric form.
• Examples:
  • Step response $y(t)$
  • Impulse response
  • Frequency response $G(j\omega)$
• Characteristics:
  • Simple and fast; require minimal prior information.
  • Generally no optimization involved; noise affects the model directly.
  • Useful for controller design methods like PID, loop shaping, and robust control.
  • Provide general insights such as bandwidth and system order.
• Role: Nonparametric models are often a necessary preliminary step for parametric identification.

2. Nonparametric Identification Methods

3. Step Response Identification

• Step input: Apply a step of amplitude $\alpha$ to the system in steady state.
• Sampling: Choose sampling period as small as possible.
• Noise: Measure noise level before applying step.

a) First-Order Model (No Delay)

• Model: $$G(s) = \frac{\gamma}{\tau s + 1}$$
• Theoretical step response: $$y(t) = K \left(1 - e^{-t/\tau}\right)$$
• Parameter relation: $$\gamma = \frac{K}{\alpha}$$
• Key points:
  • At $t = \tau$, $y(t) = 0.63 K$
  • $K$ is steady-state gain.

b) First-Order Model with Delay (Ziegler-Nichols Approximation)

• Model: $$G(s) \approx \gamma e^{-\theta s} \frac{1}{\tau s + 1}$$
• Parameters:
  • $\theta$: delay time
  • $\tau$: time constant
  • $\gamma$: gain

c) Second-Order Model

• Model: $$G(s) = \frac{\gamma \omega_n^2}{s^2 + 2 \zeta \omega_n s + \omega_n^2}$$
• Step response: $$y(t) = K - \frac{K}{\beta} e^{-\zeta \omega_n t} \sin(\omega_n \beta t + \theta)$$
• Parameters:
  • $\beta = \sqrt{1 - \zeta^2}$
  • $\theta = \arccos \zeta$
  • $\zeta$: damping ratio
  • $\omega_n$: natural frequency
  • $K$: steady-state gain

d) Identification Procedure for Second-Order Model

1. Measure steady-state gain $K$ and compute $\gamma = \frac{K}{\alpha}$.
2. Measure peak time $t_p$ and peak value $y(t_p)$.
3. Compute overshoot $M_p = y(t_p) - K$ (used to estimate damping ratio $\zeta$).

To retain: Parametric models require structural assumptions and optimization, allowing noise filtering and precise controller design; nonparametric models are simpler, faster, and provide general system insights but do not filter noise and serve as a basis for parametric identification.

1. Damping Factor and Natural Frequency

• Damping factor $\zeta$ is computed from the maximum overshoot $M_p$:

  $$M_p = e^{-\frac{\zeta \pi}{\beta}} \implies \zeta = \sqrt{\frac{\ln^2(M_p)}{\pi^2 + \ln^2(M_p)}}$$

• Natural frequency $\omega_n$ is computed from the peak time $t_p$ and $\beta$:

  $$\omega_n = \frac{\pi}{t_p \beta}$$

2. Impulse Response in Discrete and Continuous Time

• The impulse response $g(t)$ is the derivative of the step response.

• In discrete time, a scaled Kronecker delta $\alpha \delta(k)$ can be applied.

• The discrete impulse response is:

  $$g(k) = \frac{1}{\alpha} g^*(t) \big|_{t=kT}$$

• If $y(k)$ is the step response to a step of amplitude $\alpha$, then:

  $$g(k) = \frac{1}{\alpha} [y(k) - y(k-1)]$$

• Approximation linking continuous and discrete impulse responses:

  $$g(t) \approx \frac{1}{\alpha T} g^*(t)$$

3. Numerical Deconvolution to Estimate Impulse Response

• Output $y(k)$ relates to input $u(k)$ and impulse response $g(k)$ by convolution:

  $$y(k) = \sum_{j=0}^k g(j) u(k-j), \quad k=0,1,\dots,N-1$$

• Recursive formula for $g(k)$ (assuming $u(0) \neq 0$):

  $$g(k) = \frac{1}{u(0)} \left[ y(k) - \sum_{j=0}^{k-1} g(j) u(k-j) \right]$$

• For unit step input $u(k) = 1$:

  $$g(k) = y(k) - y(k-1)$$

4. Matrix Formulation and Ill-Conditioning

• Convolution can be written as a matrix equation:

  $$\mathbf{Y} = \mathbf{U} \mathbf{\Theta}$$

  where

  $$\mathbf{Y} = \begin{bmatrix} y(0) \\ y(1) \\ \vdots \\ y(N-1) \end{bmatrix}, \quad \mathbf{\Theta} = \begin{bmatrix} g(0) \\ g(1) \\ \vdots \\ g(N-1) \end{bmatrix}$$

  and $\mathbf{U}$ is a lower-triangular Toeplitz matrix built from $u(k)$.

• Problem: $\mathbf{U}$ is often ill-conditioned, making direct inversion $\mathbf{\Theta} = \mathbf{U}^{-1} \mathbf{Y}$ numerically unstable.

5. Solutions to Ill-Conditioning

Key point: The impulse response can be estimated from step or arbitrary inputs using numerical deconvolution, but care must be taken to handle noise and ill-conditioning via FIR truncation or regularization.

1. Numerical Deconvolution for Impulse Response Estimation

• For stable systems without integrators, the impulse response $g(k)$ tends to zero as $k \to \infty$. We approximate:

  $$g(k) = 0 \quad \text{for} \quad k > K$$

• The output vector $Y$ and input matrix $U_K$ relate to the truncated impulse response vector $\Theta_K$ as:

  $$Y = U_K \Theta_K$$

  where

  $$Y = \begin{bmatrix} y(0) \\ y(1) \\ \vdots \\ y(N-1) \end{bmatrix}, \quad U_K = \begin{bmatrix} u(0) & 0 & \cdots & 0 \\ u(1) & u(0) & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ u(N-1) & u(N-2) & \cdots & u(N-K) \end{bmatrix}, \quad \Theta_K = \begin{bmatrix} g(0) \\ g(1) \\ \vdots \\ g(K-1) \end{bmatrix}$$

• Since $N > K$, the system is overdetermined and solved via least squares:

  $$\Theta_K = U_K^\dagger Y = (U_K^T U_K)^{-1} U_K^T Y$$

2. Regularization: Bias-Variance Trade-off

• To stabilize the solution and control noise amplification, use Tikhonov regularization:

  $$J(\Theta) = E^T E + \lambda \Theta^T \Theta, \quad \text{with} \quad E = Y - U \Theta$$

• Minimizing $J(\Theta)$ yields:

  $$\Theta = (U^T U + \lambda I)^{-1} U^T Y$$

• $\lambda > 0$ ensures invertibility and bounds the norm of $\Theta$.

• Choose $\lambda$ by trial and error, guided by expected impulse response length.

3. Correlation Approach

a) Cross-Correlation Function

• For stochastic signals $x(k)$ and $y(k)$, the cross-correlation at lag $h$ is:

  $$R_{xy}(h) = \mathbb{E}\{ x(k) y(k - h) \} = \lim_{N \to \infty} \frac{1}{N} \sum_{k=0}^{N-1} x(k) y(k - h)$$

• If $x$ and $y$ are independent and one has zero mean:

  $$R_{xy}(h) = 0 \quad \forall h$$

b) Properties of Correlation Functions

• These follow from the definitions and change of variables in expectation sums.

Key point: The impulse response can be estimated by solving $Y = U_K \Theta_K$ via least squares, optionally regularized to balance bias and variance, and correlation functions provide fundamental properties useful for stochastic signal analysis.

1. Correlation Approach in Time-Domain Numerical Deconvolution

The output of a noisy LTI system is modeled as:

$$y(k) = \sum_{j=0}^\infty g(j) u(k-j) + d(k) = g(k) * u(k) + d(k)$$

where:

• $u(k)$ is the input (random and independent of noise),
• $d(k)$ is zero-mean noise,
• $g(k)$ is the system impulse response.

Multiplying $y(k)$ by $u(k-h)$ and taking expectation yields the cross-correlation:

$$R_{yu}(h) = E\{y(k) u(k-h)\} = \sum_{j=0}^\infty g(j) R_{uu}(h-j) + R_{du}(h)$$

Since $d(k)$ and $u(k)$ are independent, $R_{du}(h) = 0$, so:

$$R_{yu}(h) = g(h) * R_{uu}(h)$$

This convolution relation allows numerical deconvolution to estimate $g(h)$ from $R_{yu}(h)$ and $R_{uu}(h)$.

2. Estimation of Correlation Functions

a) Autocorrelation Function Estimates

• Biased estimate:

$$\hat{R}_{uu}(h) = \frac{1}{N} \sum_{k=h}^{N-1} u(k) u(k-h), \quad \hat{R}_{uu}(-h) = \hat{R}_{uu}(h)$$

This estimate is biased because:

$$E\{\hat{R}_{uu}(h)\} = \frac{N-h}{N} R_{uu}(h) \neq R_{uu}(h)$$

• Unbiased estimate:

$$\hat{R}_{uu}(h) = \frac{1}{N - h} \sum_{k=h}^{N-1} u(k) u(k-h), \quad \hat{R}_{uu}(-h) = \hat{R}_{uu}(h)$$

3. Key Properties of Autocorrelation

For any lag $h$:

$$|R_{xx}(h)| \leq R_{xx}(0)$$

This follows from the non-negativity of variance and covariance expressions.

4. Cross-Correlation Function Estimate (Biased)

The biased estimate for cross-correlation between $y(k)$ and $u(k)$ is:

$$\hat{R}_{yu}(h) = \frac{1}{N} \sum_{k=0}^{N-1 + \min(h,0)} y(k) u(k - h)$$

5. Summary Table: Correlation Estimates

Key point: The system impulse response $g(h)$ can be estimated by numerically deconvolving the measured cross-correlation $R_{yu}(h)$ with the input autocorrelation $R_{uu}(h)$, assuming noise is independent and zero-mean.

1. Correlation Approach in Time-Domain Methods

a) Cross-Correlation Function (Unbiased Estimate)

• For signals $y(k)$ and $u(k)$, the unbiased cross-correlation estimate at lag $h$ is:

$$\hat{R}_{yu}(h) = \frac{1}{N - |h|} \sum_{k=\max(h,0)}^{N-1+\min(h,0)} y(k) u(k - h), \quad h = 0, \pm 1, \pm 2, \ldots$$

b) Impulse Response Estimation

• If input $u(k)$ is white noise:

  • Cross-correlation relates to impulse response $g(h)$ as:
  $$R_{yu}(h) = g(h) * \delta(h) R_{uu}(0) = g(h) R_{uu}(0)$$
  • Since $R_{uu}(0)$ is the input variance, estimate impulse response by:
  $$g(h) = \frac{\hat{R}_{yu}(h)}{\hat{R}_{uu}(0)}$$

• If input $u(k)$ is not white noise:

  • Use numerical deconvolution assuming $g(j) = 0$ for $j \geq K$.

  • Cross-correlation expressed as convolution sum:

  $$R_{yu}(h) = \sum_{j=0}^{K-1} g(j) R_{uu}(h - j), \quad h = 0, 1, 2, \ldots$$

• This leads to a linear system for $g$:

  $$\begin{bmatrix} \hat{R}_{yu}(0) \\ \hat{R}_{yu}(1) \\ \vdots \\ \hat{R}_{yu}(K-1) \end{bmatrix} = \begin{bmatrix} \hat{R}_{uu}(0) & \hat{R}_{uu}(-1) & \cdots & \hat{R}_{uu}(-K+1) \\ \hat{R}_{uu}(1) & \hat{R}_{uu}(0) & \cdots & \hat{R}_{uu}(-K+2) \\ \vdots & \vdots & \ddots & \vdots \\ \hat{R}_{uu}(K-1) & \hat{R}_{uu}(K-2) & \cdots & \hat{R}_{uu}(0) \end{bmatrix} \begin{bmatrix} g(0) \\ g(1) \\ \vdots \\ g(K-1) \end{bmatrix}$$

• Solve this system to estimate the impulse response vector $\mathbf{g}$.

2. Summary Table: Impulse Response Estimation Methods

Key point: When $u(k)$ is white noise, impulse response estimation reduces to dividing cross-correlation by input autocorrelation at zero lag; otherwise, solve a convolution linear system.

3. Additional Notes

• The cross-correlation function is central to nonparametric identification in the time domain.

• The assumption $g(j) = 0$ for $j \geq K$ truncates the impulse response to a finite length for practical computation.

• The unbiased estimate accounts for the reduced number of terms at large lags to avoid bias.

1. Frequency-Domain Methods: Frequency Analysis

a) Fourier Transform and System Response

• For continuous-time LTI systems: $y(t) = g(t) * u(t) = \int_{-\infty}^{\infty} g(\tau) u(t-\tau) d\tau$

• Taking Fourier transform:

  $$Y(j\omega) = G(j\omega) U(j\omega) \implies G(j\omega) = \frac{Y(j\omega)}{U(j\omega)}$$

• Where:

  $$U(j\omega) = \int_{-\infty}^\infty u(t) e^{-j\omega t} dt, \quad Y(j\omega) = \int_{-\infty}^\infty y(t) e^{-j\omega t} dt$$

• For discrete-time LTI systems:

  $$y(k) = \sum_{j=0}^\infty g(j) u(k-j) = g(k) * u(k)$$

• Fourier transform:

  $$Y(e^{j\omega}) = G(e^{j\omega}) U(e^{j\omega}) \implies G(e^{j\omega}) = \frac{Y(e^{j\omega})}{U(e^{j\omega})}$$

• Where:

  $$U(e^{j\omega}) = \sum_{k=-\infty}^\infty u(k) e^{-j k \omega T}, \quad Y(e^{j\omega}) = \sum_{k=-\infty}^\infty y(k) e^{-j k \omega T}$$

b) Three Types of Errors in Fourier Analysis

c) Sampling Error

• Occurs because discrete-time frequency response $G(e^{j\omega})$ is identified instead of continuous-time $G(j\omega)$.
• Shannon Sampling Theorem condition: if satisfied, $G(e^{j\omega}) = G(j\omega)$.
• In practice, this condition is never fully met, so sampling error always exists.

d) Measurement Error

• Model with noise: $$y(k) = \sum_{j=0}^\infty g(j) u(k-j) + d(k)$$
• Frequency domain: $$Y(e^{j\omega}) = G(e^{j\omega}) U(e^{j\omega}) + D(e^{j\omega})$$
• Frequency response estimate: $$\hat{G}(e^{j\omega}) = \frac{Y(e^{j\omega})}{U(e^{j\omega})} = G(e^{j\omega}) + \frac{D(e^{j\omega})}{U(e^{j\omega})}$$
• Unbiased estimate condition: noise $d(k)$ is zero-mean, so $$E\{\hat{G}(e^{j\omega}) - G(e^{j\omega})\} = 0$$
• Variance of estimate: $$\mathrm{var}[\hat{G}(e^{j\omega})] = \frac{E\{|D(e^{j\omega})|^2\}}{|U(e^{j\omega})|^2} = \frac{\Phi_{dd}(\omega)}{\Phi_{uu}(\omega)}$$
• Interpretation: Measurement error is large where noise power spectrum $\Phi_{dd}(\omega)$ is high and input power spectrum $\Phi_{uu}(\omega)$ is low.
• Increasing data length does not reduce variance caused by measurement noise.

e) Truncation Error

• Fourier transform for aperiodic discrete-time signals: $$X(e^{j\omega}) = \sum_{k=-\infty}^\infty x(k) e^{-j k \omega T}$$
• For periodic discrete-time signals (period $M$): $$X(e^{j\omega_n}) = \frac{1}{M} \sum_{k=0}^{M-1} x(k) e^{-j k \omega_n T}, \quad \omega_n = \frac{2\pi n}{M T}, \quad n=0,\ldots,M-1$$
• Source of truncation error: For aperiodic signals, the Fourier transform integral/sum is infinite but computed only on finite interval $[0, N-1]$.
• Mitigation techniques:
  • Use excitation signals with finite energy.
  • Use periodic excitation signals.
  • Apply windowing to reduce spectral leakage.

> Truncation error arises from computing Fourier transforms on finite data segments instead of infinite intervals, causing bias unless corrected by signal design or windowing.

1. Excitation with Time-Limited Signals

• Input signal $u(k)$ is zero outside $0 \leq k < M$.

• Fourier transform can be computed exactly without truncation error:

  $$U(e^{j\omega}) = \sum_{k=0}^{M-1} u(k) e^{-j k \omega T}$$

• For a stable system, output $y(k) \approx 0$ for large $k \geq N$.

• Wait until output settles before recording data to minimize truncation error on output (only noise truncated).

2. Excitation with Periodic Signals

• Input $u(k)$ is periodic with period $M$.

• Fourier series computed exactly (no truncation error):

  $$U(e^{j \omega_n}) = \frac{1}{M} \sum_{k=0}^{M-1} u(k) e^{-j k \omega_n T}$$

  where

  $$\omega_n = \frac{2\pi n}{M T}, \quad n=0,1,\ldots,M-1$$

• Output at steady-state is almost periodic (ignoring noise):

  $$Y(e^{j \omega_n}) = \frac{1}{M} \sum_{k=L}^{L+M-1} y(k) e^{-j k \omega_n T}$$

  with $L$ larger than settling time.

• Noise reduced by averaging over $p$ periods:

  $$G(e^{j \omega_n}) = \frac{1}{p} \sum_{i=1}^p \frac{Y_i(e^{j \omega_n})}{U(e^{j \omega_n})}$$

3. Windowing and Truncation Error

• Truncation in time domain = multiplication by rectangular window $f(t)$.

• In frequency domain, this corresponds to convolution:

  $$Z(\omega) = X(\omega) * F(\omega) = \int_{-\infty}^\infty X(\nu) F(\omega - \nu) d\nu$$

• Windowing introduces two main errors:

  Error Type	Cause	Effect
  Loss of resolution	Width of main lobe (MLW)	Smoothing, reduced frequency resolution
  Spectral leakage	Amplitude of secondary lobes (SLA)	Energy leaks to neighboring frequencies

• Trade-off: Increasing MLW reduces spectral leakage but lowers resolution.

• Usually prefer increasing MLW to reduce noise and leakage (e.g., Hann window).

4. Discrete Fourier Transform (DFT)

• For a sampled signal $x(k)$ of length $N$, DFT is:

  $$X(n) = \sum_{k=0}^{N-1} x(k) e^{-j \frac{2\pi}{N} n k}, \quad n=0,1,\ldots,N-1$$

• Frequency vector associated with $X$:

  $$0, \frac{\omega_s}{N}, \frac{2\omega_s}{N}, \ldots, \frac{(N-1)\omega_s}{N}$$

Key point: Using time-limited or periodic inputs allows exact computation of Fourier transforms without truncation error, while windowing controls the trade-off between spectral leakage and resolution in frequency analysis.

1. Discrete Fourier Transform (DFT)

• DFT approximates the Fourier transform of continuous-time signals by sampling uniformly $N$ points on the Fourier transform of an aperiodic discrete-time signal.
• For discrete-time periodic signals, DFT gives $N$ times the exact value of the Fourier transform.
• Sampling frequency is denoted by $\omega_s$.

2. Practical Aspects of Fourier Analysis

• The input signal must be rich in frequency content to excite the system over a wide band; typical choices are discrete white noise or PRBS (Pseudo-Random Binary Sequence) with uniform spectrum.
• Applying multiple periods ($m$) of a periodic PRBS reduces measurement error by a factor of $m$ and avoids truncation errors.
• To ensure the output is periodic (deterministic part), allow the transient response to decay before measurement.
• If a periodic input is not possible, use window functions (Hamming, Hann, etc.) to reduce truncation errors.
• Prefer fewer frequency points with higher accuracy: split data into parts, compute individual frequency responses $G_i$, then average.

3. Spectral Analysis for Linear Systems

For a linear system:

$$y(k) = g(k) * u(k) + d(k) \implies Y(e^{j\omega}) = G(e^{j\omega}) U(e^{j\omega}) + D(e^{j\omega})$$

Using correlation functions:

$$R_{yu}(h) = g(h) * R_{uu}(h) \implies \Phi_{yu}(\omega) = G(e^{j\omega}) \Phi_{uu}(\omega)$$

Hence, the frequency response function is:

$$G(e^{j\omega}) = \frac{\Phi_{yu}(\omega)}{\Phi_{uu}(\omega)} = \frac{\sum_{h=-\infty}^{\infty} R_{yu}(h) e^{-j \omega h T}}{\sum_{h=-\infty}^{\infty} R_{uu}(h) e^{-j \omega h T}}$$

In practice, estimates are computed using DFT:

$$\hat{G}(e^{j\omega_n}) = \frac{\hat{\Phi}_{yu}(\omega_n)}{\hat{\Phi}_{uu}(\omega_n)} = \frac{\sum_{h=0}^{N-1} \hat{R}_{yu}(h) e^{-j h \frac{2\pi n}{N}}}{\sum_{h=0}^{N-1} \hat{R}_{uu}(h) e^{-j h \frac{2\pi n}{N}}}$$

with

$$\omega_n = n \frac{\omega_s}{N}, \quad n = 0, 1, \ldots, N-1$$

4. Estimating Noise Spectrum

Given the system output:

$$Y(e^{j\omega}) = G(e^{j\omega}) U(e^{j\omega}) + D(e^{j\omega})$$

Taking the expected power spectrum:

$$E\{|Y(e^{j\omega})|^2\} = E\{|G(e^{j\omega})|^2 |U(e^{j\omega})|^2\} + E\{|D(e^{j\omega})|^2\}$$

Define the output power spectral density:

$$\Phi_{yy}(\omega) = E\{|Y(e^{j\omega})|^2\}$$

Then the noise power spectral density is:

$$\Phi_{dd}(\omega) = \Phi_{yy}(\omega) - \frac{|\Phi_{yu}(\omega)|^2}{\Phi_{uu}(\omega)}$$

5. Remark

• For periodic signals, Fourier analysis and spectral analysis are equivalent.
• In this case, the frequency response can be directly computed as:

$$G(e^{j\omega_n}) = \frac{\Phi_{yu}(\omega_n)}{\Phi_{uu}(\omega_n)} = \frac{Y(e^{j\omega_n}) U(e^{-j\omega_n})}{U(e^{j\omega_n}) U(e^{-j\omega_n})} = \frac{Y(e^{j\omega_n})}{U(e^{j\omega_n})}$$

Key point: Spectral analysis estimates the system frequency response by dividing the cross-spectrum of output and input by the input power spectrum, allowing nonparametric identification of linear systems.

1. Closed-Loop Identification

Motivation:

• Some systems are unstable in open-loop.
• Open-loop operation can cause output drift.
• Identification should focus on the frequency range relevant for control.

2. Problem in Closed-Loop Identification

• Input $u$ and disturbance $d$ are correlated via feedback, i.e., $R_{du}(h) \neq 0$.
• This correlation invalidates direct spectral analysis on input/output data.
• The classical spectral estimate of the system is biased:

$$\hat{G}(e^{j\omega}) = \frac{\Phi_{yu}(\omega)}{\Phi_{uu}(\omega)} = \frac{\Phi_{ypu}(\omega) + \Phi_{du}(\omega)}{\Phi_{uu}(\omega)} = G(e^{j\omega}) + \frac{\Phi_{du}(\omega)}{\Phi_{uu}(\omega)}$$

where $\Phi_{du}(\omega) \neq 0$ introduces an error term.

3. Extreme Cases in Closed-Loop Identification

Key point: Closed-loop feedback induces correlation between input and disturbance, causing bias in spectral estimates; special methods are required to obtain unbiased system identification.

4. Summary

• Closed-loop identification is necessary when open-loop operation is not feasible.
• Direct spectral analysis fails due to feedback-induced correlation.
• Understanding the correlation structure is essential to correct or compensate the bias in estimates.

1. Parametric Model Estimation from Frequency Data

a) Identification of the Controller Inverse

• Given spectral densities: $$\Phi_{yu}(\omega) = -\Phi_{eu}(\omega), \quad \Phi_{uu}(\omega) = K(e^{j\omega}) \Phi_{eu}(\omega)$$
• The estimated transfer function (TF) is: $$\hat{G}(e^{j\omega}) = \frac{\Phi_{yu}(\omega)}{\Phi_{uu}(\omega)} = -\frac{\Phi_{eu}(\omega)}{\Phi_{uu}(\omega)} = -\frac{1}{K(e^{j\omega})}$$
• Interpretation: The inverse of the controller $K(e^{j\omega})$ is identified.

b) Identification Procedure Using Excitation Signal $s$

1. Identify TF between $s$ and $y$:

   $$\hat{G}(e^{j\omega}) \left(1 + K(e^{j\omega}) \hat{G}(e^{j\omega})\right) = \frac{\Phi_{ys}(\omega)}{\Phi_{ss}(\omega)}$$

2. Identify TF between $s$ and $u$:

   $$\frac{1}{1 + K(e^{j\omega}) \hat{G}(e^{j\omega})} = \frac{\Phi_{us}(\omega)}{\Phi_{ss}(\omega)}$$

3. Divide the two models to get:

   $$\hat{G}(e^{j\omega}) = \frac{\Phi_{ys}(\omega)}{\Phi_{us}(\omega)}$$

2. Parametric Model from Frequency Data

• Given a nonparametric frequency domain model: $$G(j\omega_k), \quad k=1,\ldots,N$$
• Goal: Find a parametric rational model: $$\hat{G}(s) = \frac{B(s)}{A(s)} = \frac{b_m s^m + \ldots + b_1 s + b_0}{a_n s^n + \ldots + a_1 s + 1}, \quad m \leq n$$
• Minimize the cost function: $$J = \sum_{k=1}^N |\varepsilon(\omega_k)|^2, \quad \varepsilon(\omega_k) = G(j\omega_k) - \frac{\hat{B}(j\omega_k)}{\hat{A}(j\omega_k)}$$
• This is a nonlinear optimization problem in the complex plane.

3. Linearization of the Error for Least Squares

• Define a linear error with respect to parameters: $$\varepsilon_\ell(\omega_k) = \hat{A}(j\omega_k) G(j\omega_k) - \hat{B}(j\omega_k), \quad k=1,\ldots,N$$
• Form the system of linear equations: $$\hat{A}(j\omega_k) G(j\omega_k) - \hat{B}(j\omega_k) = 0, \quad k=1,\ldots,N$$
• Solve in the least squares sense by separating real and imaginary parts: $$\begin{cases} \Re\{\hat{A}(j\omega_k) G(j\omega_k)\} - \Re\{\hat{B}(j\omega_k)\} = 0 \\ \Im\{\hat{A}(j\omega_k) G(j\omega_k)\} - \Im\{\hat{B}(j\omega_k)\} = 0 \end{cases} \quad k=1,\ldots,N$$
• Note: The solution parameters may be complex; separating real and imaginary parts ensures a real-valued solution.

4. Problem with Frequency Weighting in Least Squares

• The linear error relates to the original error by: $$|\varepsilon_\ell(j\omega_k)| = |a_n (j\omega_k)^n + \ldots + a_1 (j\omega_k) + 1| \cdot |\varepsilon(j\omega_k)|$$
• The weighting term $|a_n (j\omega_k)^n + \ldots + a_1 (j\omega_k) + 1|$ increases with frequency, causing:
  • Good fit at high frequencies (HF)
  • Poor fit at low frequencies (LF)

5. Iterative Solution to Address Frequency Weighting

• Use an iterative approach minimizing at iteration $i$: $$J(\theta_i) = \sum_{k=1}^N \left|\varepsilon_\ell(\omega_k, \theta_i) \frac{1}{\hat{A}(\theta_{i-1})}\right|^2$$
• Here, $\hat{A}(\theta_{i-1})$ is computed from the previous iteration.
• This balances the frequency weighting and improves fit across all frequencies.

Key point: Parametric model estimation from frequency data involves transforming a nonlinear problem into a linear least squares problem via error linearization, then correcting frequency-dependent weighting through iterative refinement.