# Fitting a transfer function with CurveFit and TransferFunctionModel
In this tutorial, we will use ``kontrol.curvefit.CurveFit`` to fit a measured transfer function.
As is mentioned in previous tutorial, ``kontrol.curvefit.CurveFit`` requires a few things, the independent variable data ``xdata``, the dependent variable ``ydata``, the model ``model``, the cost function ``cost``, and the optimizer.
They have the following signature:
```
xdata : array
ydata : array
model : func(x: array, *args, **kwargs) -> array
cost : func(args: array, model: func, xdata: array, ydata: array, model_kwargs: dict) -> float
optimizer : func(cost: func, **kwargs) -> scipy.optimize.OptimzeResult
```
To save you a lot of troubles, ``kontrol`` library provides various models classes ``kontrol.curvefit.model``, cost functions (error functions). And optimizers are readily available in ``scipy.optimize``. So we typically just need to prepare ``xdata`` and ``ydata``.
Of course, knowing what model, error function, and optimizer to use is crucial in a curve fitting task. 

This time, we will use the model ``kontrol.curvefit.model.TransferFunctionModel`` as our model this time.
This model can be defined by the number of zeros ``nzero`` and ``npole``. The parameters are simply concatenated coefficients of numerator and denominator arranging from high to low order.

Here, let's consider the transfer function

\begin{equation}
    H(s) = \frac{s^2 + 3s + 2}{s^3 + 12s^2 + 47s + 60}
\end{equation}

So the parameters we would like to recover are the coefficients ``[1, 3, 2, 1, 12, 47, 60]``.

For the sake of demonstration, let's assume that we know there are 2 zeros and 3 poles as this is required to define a model.

We will use ``kontrol.curvefit.error_func.tf_error`` as the error function and it is defined as
\begin{equation}
E_\mathrm{tf\_error}\left(H_1(f), H_2(f); w(f), \epsilon\right) = \frac{1}{N}\sum_{i=0}^{N} \log_{10}(\lvert(H_1(f_i) - H_2(f_i)+\epsilon)w(f_i)\rvert)\,,
\end{equation}
where $H_1(f)$ and $H_2(f)$ are the frequency responses values (complex array) of the measured system and the model, $\epsilon$ is a small number to prevent $\log_{10}$ from exploding, $w(f)$ is a weighting function, and $N$ is the total number of data points.


```python
import control
import numpy as np
import matplotlib.pyplot as plt

f = np.logspace(-3, 3, 100000)
s = 1j*2*np.pi*f
tf = (s**2 + 3*s + 2) / (s**3 + 12*s**2 + 47*s + 60)

plt.figure(figsize=(15, 5))
plt.subplot(121)
plt.loglog(f, abs(tf))
plt.grid(which="both")
plt.ylabel("Amplitude")
plt.xlabel("Frequency (Hz)")
plt.subplot(122)
plt.semilogx(f, np.angle(tf))
plt.grid(which="both")
plt.ylabel("Phase (rad)")
plt.xlabel("Frequency (Hz)")
plt.show()
```

Now, let's everything that ``kontrol.curvefit.CurveFit`` needs, namely ``xdata``, ``ydata``, ``model``, ``cost``, and ``optimizer``.


```python
import kontrol.curvefit
import scipy.optimize


xdata = f
ydata = tf
model = kontrol.curvefit.model.TransferFunctionModel(nzero=2, npole=3, log_args=False)
error_func = kontrol.curvefit.error_func.tf_error
cost = kontrol.curvefit.Cost(error_func=error_func)
optimizer = scipy.optimize.minimize
```

Since we're using ``scipy.optimize.minimize``, it requires an initial guess. Let's start with all ones.


```python
x0 = np.ones(7)  # There are 7 coefficients, 3 for numerator and 4 for denominator
optimizer_kwargs = {"x0": x0}
```

Now let's use ``kontrol.curvefit.CurveFit`` to fit the data.


```python
a = kontrol.curvefit.CurveFit()
a.xdata = xdata
a.ydata = ydata
a.model = model
a.cost = cost
a.optimizer = optimizer
a.optimizer_kwargs = optimizer_kwargs
res = a.fit()
```


```python
plt.figure(figsize=(15, 5))
plt.subplot(121)
plt.loglog(f, abs(tf), label="Data")
plt.loglog(f, abs(a.model(f, x0)), label="Initial guess")
plt.loglog(f, abs(a.yfit), label="Fit")
plt.legend(loc=0)
plt.grid(which="both")
plt.ylabel("Amplitude")
plt.xlabel("Frequency (Hz)")
plt.subplot(122)
plt.semilogx(f, np.angle(tf), label="Data")
plt.semilogx(f, np.angle(a.model(f, x0)), label="Initial guess")
plt.semilogx(f, np.angle(a.yfit), label="Fit")
plt.legend(loc=0)
plt.grid(which="both")
plt.ylabel("Phase (rad)")
plt.xlabel("Frequency (Hz)")
plt.show()
```


```python
res.x/res.x[0]
```




    array([ 1.        ,  1.78627675,  0.84976028,  1.00034299, 11.024541  ,
           35.31604716, 25.492815  ])



Looks great, right? But, this is not typically that easy. This above case was easy because I have set the frequency array to logspace whereas we typically get a linspace array from fourier transform of linearly spaced time domain data!

Now let's see what happens if we switch back to linspace.


```python
f = np.linspace(0.001, 1000, 100000)
s = 1j*2*np.pi*f
tf = (s**2 + 3*s + 2) / (s**3 + 12*s**2 + 47*s + 60)

xdata = f
ydata = tf
model = kontrol.curvefit.model.TransferFunctionModel(nzero=2, npole=3, log_args=False)
error_func = kontrol.curvefit.error_func.tf_error
cost = kontrol.curvefit.Cost(error_func=error_func)
optimizer = scipy.optimize.minimize

x0 = np.ones(7)  # There are 7 coefficients, 3 for numerator and 4 for denominator
optimizer_kwargs = {"x0": x0}

a = kontrol.curvefit.CurveFit()
a.xdata = xdata
a.ydata = ydata
a.model = model
a.cost = cost
a.optimizer = optimizer
a.optimizer_kwargs = optimizer_kwargs
res = a.fit()
```


```python
plt.figure(figsize=(15, 5))
plt.subplot(121)
plt.loglog(f, abs(tf), label="Data")
plt.loglog(f, abs(a.model(f, x0)), label="Initial guess")
plt.loglog(f, abs(a.yfit), label="Fit")
plt.legend(loc=0)
plt.grid(which="both")
plt.ylabel("Amplitude")
plt.xlabel("Frequency (Hz)")
plt.subplot(122)
plt.semilogx(f, np.angle(tf), label="Data")
plt.semilogx(f, np.angle(a.model(f, x0)), label="Initial guess")
plt.semilogx(f, np.angle(a.yfit), label="Fit")
plt.legend(loc=0)
plt.grid(which="both")
plt.ylabel("Phase (rad)")
plt.xlabel("Frequency (Hz)")
plt.show()
```

Not great! Terrible!

It turns out that there're a few tricks that can be used fit linspace transfer function data.

1. Weighing function (1/f).
2. Parameter scaling (log).
3. Use other optimizers (differential evolution, Nelder-Mead, Powell, etc).
4. Use better initial guess.
5. Tighter acceptable tolerance for convergence during minimization (e.g. ``ftol``, ``xtol``. Different optimizers use different convention)

We will use all of them!


```python
a = kontrol.curvefit.CurveFit()
a.xdata = f
a.ydata = tf
# Use 1/f weighting
weight = 1/f
error_func_kwargs = {"weight": weight}
a.cost = kontrol.curvefit.Cost(error_func=error_func, error_func_kwargs=error_func_kwargs)

# Use log scaling. Models have an argument for enabling that.
a.model = kontrol.curvefit.model.TransferFunctionModel(nzero=2, npole=3, log_args=True)

# Use better initial guess and Nelder-Mead optimizer.
np.random.seed(123)
true_args = np.array([1, 3, 2, 1, 12, 47, 60])  # These are the true parameters
noise = np.random.normal(loc=0, scale=true_args/10, size=7)
x0 = true_args + noise  # Now the initial guess is assumed to be some deviation from the true values.
# x0 = np.ones(7)
a.optimizer_kwargs = {"x0": np.log10(x0), "method": "powell", "options": {"xtol": 1e-12, 'ftol': 1e-12}}  # Note the inital guess is log.
a.optimizer = scipy.optimize.minimize
res = a.fit()
```


```python
10**a.optimize_result.x
```




    array([ 2.25124892,  3.66597884,  1.51466572,  2.25137845, 23.92566764,
           77.41975222, 45.43996825])




```python
plt.figure(figsize=(15, 5))
plt.subplot(121)
plt.loglog(f, abs(tf), label="Data")
plt.loglog(f, abs(a.model(f, np.log10(x0))), label="Initial guess")
plt.loglog(f, abs(a.yfit), label="Fit")
plt.legend(loc=0)
plt.grid(which="both")
plt.ylabel("Amplitude")
plt.xlabel("Frequency (Hz)")
plt.subplot(122)
plt.semilogx(f, np.angle(tf), label="Data")
plt.semilogx(f, np.angle(a.model(f, np.log10(x0))), label="Initial guess")
plt.semilogx(f, np.angle(a.yfit), label="Fit")
plt.legend(loc=0)
plt.grid(which="both")
plt.ylabel("Phase (rad)")
plt.xlabel("Frequency (Hz)")
plt.show()
```


```python
fitted_args = 10**a.optimized_args
fitted_args /= fitted_args[0]
fitted_args
```




    array([ 1.        ,  1.62842003,  0.6728113 ,  1.00005754, 10.62773089,
           34.38968995, 20.18433765])




```python
print("True arguments: ", true_args)
print("Initial arguments: ", x0/x0[0])
print("Arguments optained from fitting: ", fitted_args)
```

    True arguments:  [ 1  3  2  1 12 47 60]
    Initial arguments:  [ 1.          3.70099497  2.30705685  0.95281056 12.68253445 61.43087558
     50.97379581]
    Arguments optained from fitting:  [ 1.          1.62842003  0.6728113   1.00005754 10.62773089 34.38968995
     20.18433765]



```python
print("True transfer function")
s = control.tf("s")
tf = (s**2 + 3*s + 2) / (s**3 + 12*s**2 + 47*s + 60)
tf
```

    True transfer function





$$\frac{s^2 + 3 s + 2}{s^3 + 12 s^2 + 47 s + 60}$$




```python
print("Fitted transfer function")
a.model.args = a.optimized_args
a.model.tf.minreal()
```

    Fitted transfer function





$$\frac{0.9999 s^2 + 1.628 s + 0.6728}{s^3 + 10.63 s^2 + 34.39 s + 20.18}$$



The fit looks good.
But, apparently it has found a set of different parameters.
See further tutorials for alternative ways of fitting the same transfer function.

Nevertherless, we demonstrated the use of ``kontrol.curvefit.CurveFit`` for transfer function fitting.

52ead2c38df00c1f65095afaf5b6529b158ea757
ipynb
Jupyter Notebook
docs/source/tutorials/curve_fitting_transfer_function.ipynb
terrencetec/kontrol
ba6461784e38d01399efeb7a42911259f9254db0
MIT
2020-08-31T10:34:43.000Z
2021-08-10T20:48:59.000Z
docs/source/tutorials/curve_fitting_transfer_function.ipynb
terrencetec/kontrol
ba6461784e38d01399efeb7a42911259f9254db0
MIT
2020-06-16T18:38:25.000Z
2022-03-24T00:48:55.000Z
docs/source/tutorials/curve_fitting_transfer_function.ipynb
terrencetec/kontrol
ba6461784e38d01399efeb7a42911259f9254db0
MIT
Qwen/Qwen-72B
1. YES
2. YES