Estimating parameters from the step response

Often, our model has a couple parameters we don’t know well from the specifications, but must attempt to measure. For the system under consideration, perhaps the two parameters most interesting to measure are the dominant time constant and the transformer ratio \(TF\) (most important). In this section, we explore how one might estimate them from a measured step response. Other parameters in the system could be similarly estimated.

By way of the transfer function, the state-space model can be transformed into input-output differential equations.

syms B_ J_ TF_ L_ R_ Vs_ s % using underscore for syms

a_ = [-B_/J_,TF_/J_;-TF_/L_,-R_/L_];
b_ = [0;1/L_];

(s*eye(2)-a_)^-1*b_

ans =
         TF_/(TF_^2 + B_*R_ + B_*L_*s + J_*R_*s + J_*L_*s^2)
 (B_ + J_*s)/(TF_^2 + B_*R_ + B_*L_*s + J_*R_*s + J_*L_*s^2)

The differential equation for \(\Omega_J\) is \[\begin{aligned} \frac{d^2 \Omega_J}{d t^2} + \left(\frac{R}{L} + \frac{B}{J}\right) \frac{d \Omega_J}{d t} + \frac{TF^2 + B R}{J L} \Omega_J &= \frac{TF}{J L} V_s. \end{aligned}\] The corresponding characteristic equation is \[\begin{aligned} \lambda^2 + \left(\frac{R}{L} + \frac{B}{J}\right) \lambda + \frac{TF^2 + B R}{J L} = 0 \end{aligned}\] which has solution \[\begin{aligned} \lambda_{1,2} = -\frac{1}{2}\left(\frac{R}{L} + \frac{B}{J}\right) \pm \frac{1}{2} \sqrt{\left(\frac{R}{L} + \frac{B}{J}\right)^2 - 4\frac{TF^2 + B R}{J L}}. \end{aligned}\] For a step input \(V_s(t) = \overline{V}_s\), \(\Omega_J(0) = d\Omega_J(0)/dt = 0\), and distinct roots \(\lambda_1\) and \(\lambda_2\), the solution is \[\begin{aligned} \label{eq:forced2} \Omega_J(t) = \overline{V}_s \frac{TF}{TF^2 + B R} \left( 1 - \frac{1}{\lambda_2-\lambda_1} \left( \lambda_2 e^{\lambda_1 t} - \lambda_1 e^{\lambda_2 t} \right) \right) \end{aligned}\]

Let’s compute \(\lambda_1\) and \(\lambda_2\).

lambda12 = -1/2*(R/L+B/J) + ...
    [1,-1]*1/2*sqrt((R/L+B/J)^2 - 4*(TF^2+B*R)/(J*L))

lambda12 =
  -16.3467 -373.9941

Both values are real, so we expect not an oscillation, but a decay to a final value. However, that decay occurs with two different time constants: \(\tau_1 = -1/\lambda_1\) and \(\tau_2 = -1/\lambda_2\).

tau12 = -1./lambda12
disp(['ratio: ',num2str(tau12(1)/tau12(2))])

tau12 =
    0.0612    0.0027
ratio: 22.8788

So second decays much faster than the first. That’s good news for our estimation project because we can easily ignore the step response’s first \(5\tau_2 \approx 0.0134\) s and assume the rest is decaying at \(\tau_1\), which we call the dominant time constant and which we would like to estimate.

Let’s generate some fake response data to get the idea. We’ll layer on some Gaussian noise with randn to be more realistic. The data is plotted in figure 4.16.

t_data = linspace(0,-6/lambda12(1),200);

O_fun = @(t) TF/(TF^2+B*R)*...
    (1-1/(lambda12(2)-lambda12(1))*...
    (lambda12(2)*exp(lambda12(1)*t)-...
    lambda12(1)*exp(lambda12(2)*t)));
rng(2);
O_data = O_fun(t_data) + .5*randn(size(t_data));

Let’s trim the data to eliminate the time interval corresponding to the first five of the “fast” time constant \(\tau_2\).

[t_5,i_5] = min(abs(t_data-(-5/lambda12(2)))); % delete
t_data_trunc = t_data((i_5+1):end);
O_data_trunc = O_data((i_5+1):end);

We need want to take the natural logarithm of the data so we can perform a linear regression to estimate the “experimental” slow time constant \(\tilde{\tau}_1\). We must first estimate the steady-state value \(\Omega_{J\infty}\) (which we’ll also need). We don’t want to just take the last value in the array due to its noisiness. The data goes for six slow time constants, so averaging the data for the last time constant is a good estimate.

[t_ss,i_ss] = ...
    min(abs(t_data_trunc-(-5/lambda12(1)))); % start here
O_data_ss = O_data_trunc((i_ss+1):end);
mu_O_ss = mean(O_data_ss)
S_mu_O_ss = std(O_data_ss)/sqrt(length(O_data_ss))

mu_O_ss =
   10.1801
S_mu_O_ss =
    0.0763

Let’s use this result to transform the data into its linear form.

O_lin = log(-(O_data_trunc-mu_O_ss));
O_lin_complex = find(imag(O_lin)>0);
disp(['number of complex values: ',...
    num2str(length(O_lin_complex))])

number of complex values: 33

Now we have encountered a problem. The noisiness of the data makes some of our points wander into negative-land. Logarithms of negative numbers are complex. Naive approaches like just taking real parts, excluding complex values, or coercing complex values to \(-\infty\) all have the issue of biasing the data.

There are a lot of approaches we could take. The best approaches include nonlinear regression and discrete filtering to smooth the data (e.g. filtfilt).

We opt for an easier approach: we find the index at which the time series first transgresses the boundary and exclude the data beyond the previous index.

i_bad = O_lin_complex(1);
t_lin_trunc = t_data_trunc(1:i_bad-1);
O_lin_trunc = O_lin(1:i_bad-1);

This is plotted in figure 4.17 along with the linear regression least-squares fit, computed below.

pf = polyfit(t_lin_trunc,O_lin_trunc,1);
O_lin_fit = polyval(pf,t_lin_trunc);
tau_1_est = -1/pf(1)

tau_1_est =
    0.0603

So our estimate for \(\tau_1\) is \(\tilde{\tau}_1 = 60.3\) ms. Recall that our analytic expression for \(\tau_1\) is known in terms of other parameters. Similarly, the steady-state value of \(\Omega_J\), which has already been estimated to be \(\Omega_{J\infty} = 10.18\) (i.e. mu_O_ss). This occurs when the time-derivatives of \(\Omega_J\) are zero. From the solution for \(\Omega_J\) (or its differential equation), for constant \(V_s(t) = \overline{V_s}\), this occurs when \[\begin{aligned} \label{eq:omega_J_ss} \Omega_{J\infty} &= \frac{TF}{TF^2+B R} \overline{V_s}. \end{aligned}\] An analytic expression for \(TF\) can be found by solving , which yields \[\begin{aligned} TF = \overline{V_s} \pm \frac{1}{2\tilde{\Omega}_{J\infty}}\sqrt{V_s^2-4 B R \tilde{\Omega}_{J\infty}^2} \end{aligned}\]

We choose the solution closer to the a priori (spec) value of \(0.0974\).

TF_est = (1 + (- 4*B*R*mu_O_ss^2 + 1^2)^(1/2))/(2*mu_O_ss)

TF_est =
    0.0976

This estimate \(\tilde{TF} = 0.0976\) is very close to the value given in the specification sheet because we constructed it to be so. Real measurements would probably yield an estimate further from the specification, which is why we would estimate it.