In engineering we often develop parametrised mathematical models that describe a system and then make measurements on the real system in order to determine the model parameters. Parameter estimation or fitting is the process that we use to extract the model parameters from the system measurements.

This page gives a short introduction to parameter estimation (or curve fitting) using MATLAB^{®} by presenting two examples a linear fit (also known as linear regression) and a nonlinear fit to the transfer function of a low pass filter. I also provide a link to a short video tutorial showing how to do parameter estimation in Excel.

## Linear fit: Ohms Law

Consider the results obtained from measuring current (I) through a resistor as the voltage (V) across the resistor is varied. Here the system model is well described by Ohms law I=(1/R)V where the resistance R is the unknown parameter.

In order to set the problem up in MATLAB^{®} where first need to define the measured data in our main mfile.

xdata=[0 1 2 3 4 5 6 7 8 9 10]; % Input Voltage (V) ydata=[1.66 4.55 5.82 9.25 11.59 13.60 16.36 19.10 22.57 24.83 26.09]; % Current (mA)

We then need to define the linear model which is a function of the variable x and has parameters m and c, corresponding to the slope and y offset respectively.

We will write equation (1) as a function in a separate mfile.

function y = Linear(p,x) % General Linear Model e.g. Ohms Law m=p(1); c=p(2); y=m*x+c;

We then also need to define a cost function. The cost function defines the quality of the fit between the model and the measured data. By minimising the cost function with respect to the parameter we can obtain the optimum parameter value and hence best fit. The most commonly used cost function is the least squares cost function which is defined as the sum over all measurement points of the square of the difference between the measured value and the model evaluated at the measurement point.

(2)

Again we define the cost function equation (2) as a function in a separate mfile.

function error=LSCostFunction(p,model,xdata,ydata) % Fit cost function - Least squares % Evaluate model at x data points yfit=model(p,xdata); % Calculate Least squares error error=sum((ydata-yfit).^2);

We can now use use the builtin MATLAB^{®} function **fminsearch** to find the minima of the cost function and return the parameter that provides the best fit. To do this we first set the initial guess for the model parameters in the main mfile. Choosing an initial value that is close to the right value will speed up the minimisation and make the algorithm more robust to avoiding local minima. In this example, you might use the values that you obtained from independent measurements of the resistance using an ohm meter.

% Set initial guess for model parameters p_initial = [1,0];

We then call the** fminsearch** function. Note, that we pass it the function handle to the cost function **LSCostFunction** The handle (like a pointer in C) is indicated by the @. The (p) after the @ tells the **fminsearch** function that this is the parameter that it should minimise with respect to. In the case of a linear regression there are two parameters the slope and the constant so the search is a two dimensional search. We also need to pass the **LSCostFunction** the: pararmeters** p**; the model function, again this is passed as a function handle **@LinearModel**; and the data **xdata** and **ydata**.

% Find minima of cost function p_fit=fminsearch(@(p) LSCostFunction(p,@LinearModel,xdata,ydata),p_initial);

Finally we evaluate the model using the parameters obtained from the minimisation of the cost function and plot out the measurement results and the fitted model. You now see why it was nice to define the model as a function in a separate mfile.

% Evaluate the model using the optimum parameters obtained xfit = linspace(0,10,100); yfit = LinearModel(p_fit,xfit); % Plot the results figure(1) plot(xdata,ydata,'x',xfit,yfit) xlabel('Input Voltage (V)') ylabel('Current (mA)') title(['Linear fit m = ' num2str(p_fit(1)) ' c = ' num2str(p_fit(2))])

In MATLAB^{®} polynomial fits such as this can also be done using the **polyfit** function. In this case you do not need to define the model or the cost function. The x and y data is first defined as before and the order of the polynomial is also defined.

% Polynomial order order = 1; %Linear % Determine the parameters p_fit = polyfit(xdata,ydata,order);

We can then evaluate the polynomial using **polyval** and the parameters obtained from the **polyfit** function. Finally we plot out the measurement results and the fitted model.

% Evaluate the model using the optimum parameters obtained xfit = linspace(0,10,100); yfit = polyval(p_fit,xfit); % Plot the results figure(2) plot(xdata,ydata,'x',xfit,yfit) xlabel('Input Voltage (V)') ylabel('Current (mA)') title(['Linear fit m = ' num2str(p_fit(1)) ' c = ' num2str(p_fit(2))])

## Nonlinear fit: Low pass filter

The nonlinear power transfer function of a RC low pass filter in the frequency domain is given by

(3)

It is straightforward to fit this to a measured frequency response of a real filter in MATLAB^{®} using **fminsearch**. Firstly we define the measured data in our main mfile.

% Measured data from RC filter fdata=[10,14.3,20.6,29.7,42.8,6 1.5,88.5,127,183,263,379,545,784,1.12e3,... 1.62e3,2.33e3,3.35e3,4.83e3,6.95e3,1e4]; Gdata=[1.02,1.01,0.999,0.998,0.994,0.987,0.969,0.937,0.880,0.773,0.627,... 0.448,0.282,0.157,0.0838,0.0433,0.0239,0.0107,0.00849,0.00290];

We then define the model which in this case is the transfer function given in equation 3 as a function in a separate mfile. Note in this case the frequency f is the variable and the capacitance C is the unknown parameter. The resistance R is a known constant.

function y = LowPassFilter(p,f) % Low pass RC filter power transfer function R=990; % Resistance (ohms) C=p(1); y=1./(1+(R*C*2*pi*f).^2);

We can reuse the least squares cost function defined for the linear fit and just need to call it and the transfer function defined above from within **fminsearch**. Firstly we need to define the initial guess for the value of the capacitance.

% Set initial guess for the capacitor value. p_initial = 0.3e-6; % Capacitance(F) % Find minima of cost function. Note the LowPassFilter model function has been used p_fit=fminsearch(@(p) LSCostFunction(p,@LowPassFilter,fdata,Gdata),p_initial);

Finally we can evaluate the model using the value of the capacitance determined by the minimisation of the cost function and plot the results. Here we plot the gain in dB on the y axis and the frequency in (Hz) is plotted on a log scale to better show the dynamic range of the filter response.

% Evaluate the model using the optimum parameters obtained xfit = logspace(1,4,100); yfit = LowPassFilter(p_fit,xfit); % Plot the results in dB on a log scale figure(3) semilogx(fdata,10*log10(Gdata),'x',xfit,10*log10(yfit)) grid xlabel('Frequency (Hz)') ylabel('Gain (dB)') title(['Low pass filter C = ' num2str(p_fit(1)) ' F'])

Nonlinear Fit Problem: RLC Filter

The code snippet below provides a set of data of the current flowing through the a simple series RLC filter under sinusoidal excitation as a function of frequency. Use this data to determine the value of the resistance in the circuit.

% This data was obtained using a series RLC filter board from the first year experiment fdata=logspace(4,6,20); Idata=[0.000436,0.000563,0.000989,0.00152,0.00267,0.00524,0.0125,0.0439,... 0.999,0.109,0.0203,0.00755,0.00361,0.00194,0.001118,0.000645,0.000450,... 0.000246,0.000324,0.000119];

You can download the complete set of mfiles for these examples here as a zip file.

## Using Excel

Nonlinear curve fitting and parameter estimation can also be carried out in Excel. I have created multimedia tutorial that walks you through the steps required to do similar examples in Excel. The Excel spreadsheet used in the tutorial is available here.

Pingback: Using Energia to code in c » Hallo, here is Pane!

Hi! Thanks for this blog, it was really helpful for me. However can you please tell us how to estimate parameters in the first case you presented , subject to the constraint that errors have to be positive?

Use x = fminbnd(fun,x1,x2) where x1 is the lower bound and x2 is the upper bound, instead of fminsearch