Issue with Fitting Damping Model to Experimental Data

hello everyone,

I am trying to fit experimental data q(t) (response over time t) to a damping model. The model describes the system's behavior in three damping cases: overdamped, critically damped, and underdamped. Parameters R, L, and C are calculated from the data using a matrix-based approach.

Method for Calculating R, L, C:

1. From the data q(t), the following values are derived:

• q_i,
• dq_i/dt ,
• d²q_i/dt².

1. A matrix equation is formed:cssCopy codeA * X = V

Where:

• A = [ q_i, dq_i/dt, d²q_i/dt² ],
• X = [ 1/C, R, L ],
• V = the input voltage signal.

1. Solving this matrix equation provides R, L, and 1/C.

The damping behavior is determined based on the condition:

mathematcal

Damping Condition = R² - (4L/C)

1. Overdamped (R² > 4L/C):perlCopy codeq(t) = A1 * exp(s1 * t) + A2 * exp(s2 * t)

Where

s1 = -(R / 2L) + sqrt((R² / 4L²) - (1 / LC)),

s2 = -(R / 2L) - sqrt((R² / 4L²) - (1 / LC))

1. Critically Damped (R² = 4L/C):perlCopy codeq(t) = (A1 + A2 * t) * exp(-2L / R * t)

1. Underdamped (R² < 4L/C):perlCopy codeq(t) = exp(-2L / R * t) * (A1 * cos(omega_d * t) + A2 * sin(omega_d * t))

Where:

omega_d = sqrt((1 / LC) - (2L / R)²)

1. Complex Numbers in All Cases:Regardless of the damping condition, complex numbers sometimes appear during calculations or model evaluation. For example, in the underdamped case, if:scssCopy code(1 / LC) - (2L / R)² < 0

The value of omega_d becomes imaginary, making the model invalid for real data.

1. Poor Curve Fitting:I am using a least-squares method to fit the curve:scssCopy codeError = sum((q_obs(t) - q_pred(t))²)

However, the fit often fails to match the experimental data or converges to incorrect parameter values. This issue occurs across all damping cases.

[Torsten]

I don't understand what you try to fit against what.

My guess is that you want to fit your experimental data q(t) to one of the functional equations codeq(t), but I don't understand where you need the matrix equation A*X = V and where you have the V-data. Further what are A1 and A2 in codeq(t) ? Do they have to do with the matrix A from A*X = V ?

[Ts]

• What am I trying to fit?
• I am fitting the experimental data q(t) to one of the functional equations for q(t) based on the damping behavior: overdamped, critically damped, or underdamped. These equations depend on system parameters R, L, and C, as well as fitting parameters A1 and A2.
• Why do I use the matrix equation A * X = V?T
• he matrix equation is used to estimate the parameters R, L, and C from the input signal V(t) provided as a time series. Here's how it works:
• The cumulative sum of V(t) is used to compute q(t): q(t) = sum(V(t) * delta_t)
• The first and second numerical derivatives of q(t) are calculated as: dq/dt ≈ (q(t + delta_t) - q(t)) / delta_t
• d²q/dt² ≈ (dq/dt at t+delta_t - dq/dt at t) / delta_t
• Using the relationship between q(t), its derivatives, and V(t): V(t) = q(t)/C + R * (dq/dt) + L * (d²q/dt²)

This can be written in matrix form as:

A * X = V

Here, A contains q(t), dq/dt, and d²q/dt², X contains the unknown parameters [1/C, R, L], and V is the observed input voltage.

Solving this equation using least squares gives estimates for R, L, and C.

• What are A1 and A2 in q(t)?
• A1 and A2 are fitting parameters used in the functional form of q(t) that describe the damping behavior. These parameters are independent of the matrix A. Once R, L, and C are estimated, A1 and A2 are determined through a curve-fitting method to minimize the error between the observed and predicted q(t) values.

[Torsten]

Please include your code so far. The second column of the Excel file is q(t) ? Where is V(t) ?

[Ts]

% --- Step 1: Load and prepare data ---

[file, path] = uigetfile('*.csv', 'Select a CSV file');

data = readtable(fullfile(path, file));

column_name = input('Enter the column name for the price data: ', 's');

time_column = input('Enter the column name for the time data (or leave blank if no time column): ', 's');

V_data = data{:, column_name};

if ~isempty(time_column) && ismember(time_column, data.Properties.VariableNames)

time_data = datetime(data{:, time_column}, 'ConvertFrom', 'excel');

t_numeric = days(time_data - time_data(1));

else

t_numeric = (0:length(V_data)-1)';

end

delta_t = t_numeric(2) - t_numeric(1);

V_data = (V_data - mean(V_data)) / std(V_data);

% Apply smoothing to reduce noise

V_data_smoothed = smoothdata(V_data, 'sgolay', 11);

figure;

plot(t_numeric, V_data_smoothed, 'b');

xlabel('Time');

ylabel('Voltage (V)');

title('Smoothed Voltage Data');

grid on;

% --- Step 2: Perform wavelet decomposition of the signal ---

[coeffs, lengths] = wavedec(V_data_smoothed, 3, 'db4');

reconstructed_signals = {};

for i = 1:length(lengths)

coeff_temp = zeros(size(coeffs));

idx_start = sum(lengths(1:i-1)) + 1;

idx_end = sum(lengths(1:i));

coeff_temp(idx_start:idx_end) = coeffs(idx_start:idx_end);

reconstructed_signal = waverec(coeff_temp, lengths, 'db4');

reconstructed_signal = (reconstructed_signal - mean(reconstructed_signal)) / std(reconstructed_signal);

reconstructed_signals{i} = reconstructed_signal;

end

% --- Solve RLC ---

function [q_i, C_i, R_i, L_i] = solve_rlc_with_charge(V_signal, delta_t)

q_i = cumsum(V_signal) * delta_t;

dq_i_dt = gradient(q_i, delta_t);

d2q_i_dt2 = gradient(dq_i_dt, delta_t);

A = [q_i, dq_i_dt, d2q_i_dt2];

I = eye(size(A, 2));

X = (A' * A + I) \ (A' * V_signal);

C_i = max(1 / X(1), 1e-6);

R_i = max(X(2), 1e-6);

L_i = max(X(3), 1e-6);

end

% --- Define RLC Model ---

function q_t = rlc_charge_model(t, R, L, C, A1, A2)

R_squared = R^2;

LC_ratio = 4 * L / C;

if R_squared > LC_ratio

% Overdamped case

s1 = (-R + sqrt(R_squared - LC_ratio)) / (2 * L);

s2 = (-R - sqrt(R_squared - LC_ratio)) / (2 * L);

q_t = A1 * exp(s1 * t) + A2 * exp(s2 * t);

elseif R_squared == LC_ratio

% Critically damped case

s = -R / (2 * L);

q_t = (A1 + A2 * t) .* exp(s * t);

else

% Underdamped case

alpha = -R / (2 * L);

omega_d = sqrt(1 / (L * C) - (R / (2 * L))^2);

q_t = exp(alpha * t) .* (A1 * cos(omega_d * t) + A2 * sin(omega_d * t));

end

end

% --- Function to minimize (error function) ---

function error = error_function(params, t, q_actual, R, L, C)

A1 = params(1);

A2 = params(2);

q_t = rlc_charge_model(t, R, L, C, A1, A2);

error = sum((q_t - q_actual).^2);

end

% --- Compute R, L, C and reconstruct q_t for each sub-signal ---

for i = 1:length(reconstructed_signals)

signal = reconstructed_signals{i};

[q_i, C_i, R_i, L_i] = solve_rlc_with_charge(signal, delta_t);

% Optimize A1 and A2 for this sub-signal using q_i

opt_params = fminunc(@(params) error_function(params, t_numeric, q_i, R_i, L_i, C_i), [0, 0]);

A1_opt = opt_params(1);

A2_opt = opt_params(2);

% Reconstruct q_t for this sub-signal

q_t_signal = rlc_charge_model(t_numeric, R_i, L_i, C_i, A1_opt, A2_opt);

% Plot each sub-signal with its q_t

figure;

plot(t_numeric, q_i, 'b', 'DisplayName', sprintf('q_{actual} (Signal %d)', i));

hold on;

plot(t_numeric, q_t_signal, 'r--', 'DisplayName', sprintf('q_t (Model for Signal %d)', i));

xlabel('Time (s)');

ylabel('Charge (q)');

title(sprintf('Sub-Signal %d: q_{actual} vs q_t', i));

legend;

grid on;

% Print results for this sub-signal

fprintf('Signal %d: C = %.4f, R = %.4f, L = %.4f, A1 = %.4f, A2 = %.4f\\n', i, C_i, R_i, L_i, A1_opt, A2_opt);

end

• How can I adjust the parameters so that the two curves match?I am trying to fit the experimental data q(t) to the predicted model of q(t). However, the two curves do not align well. What steps or adjustments should I make to ensure the two curves align accurately?
• What is the explicit equation for q(t) at each i?Given that q(t) is derived from the input signal V(t), what would be the explicit equation for q(t) at each time step i, considering the relationships between q(t), its derivatives, and V(t)?

Answer 1

[William Rose]

@Ts,

You say q(t) is a response. Response to what? To v(t), I assume. But you have not told us v(t). If you apply the differential equations for linear second order system you your data, without a v(t) input, then it is as if v(t)=0. Then you can adjust the R,L, and C, and you can adjust the initial conditions. The result will be a decaying exponential or a sum of 2 decaying exponentials, a damped sinusoid, or an undamped sinusoid (if R=0). You will never get a response that looks much like like the plot you shared.

If you want to estimate v(t) from q(t), the problem is under-determined, or it needs more constraints to make it interesting. After all, you could create a v(t) waveform that generates the observed q(t) with a single capacitor: .

Good luck with your fitting.

[William Rose]

@Ts, thank you for posting the additional file which includes the voltage data.

I am currently unable to run the code below in the Matlab Answer window, because the data file won't load, even though I attach it with the paper clip. I don't know why it won't load. I have tried a few things. But I can run the code on my computer. Therefore I will insert figures showing the results I get. Sorry for the inconvenience.

data=load('combinedData.txt');
t=data(:,1); v=data(:,2); q=data(:,3);
% Plot data
figure; subplot(211), plot(t,q,'-r.')
grid on; ylabel('q')
subplot(212), plot(t,v,'b.')
grid on; ylabel('v'); xlabel('Time')

Code above produces figure below.

The signal v(t) has small fluctuations on top of a large mean value. v(t) shows weird quantization. Weird in the sense that the step size between levels is sometimes large and sometimes much smaller.

Let's assume this is voltage across a circuit which is an inductor L, resistor R, and capacitor C, in series. Then the total voltage across the circuit is

where v0 is a voltage offset. (Voltage can be have an arbitrary offset without affecting the way the circuit works.)

You make a very good point, which is that the system is linear in the unknown parameters L, R, C, and v0. Therefore can use the standard matrix method to find the least squares solution. We make an array qAll, with four columns for the four predictor variables: d2q/dt2, dq/dt, q, and a column of ones for the intercept term, v0. We cut 2 points off q, and 1 point off dqdt, so that all the columns of matrix qAll are the same length. We also trim the time vector and the vector v to the same length. The trimmed t and q and v are called tt, qt, and vt. I multiply q by 1e14, so it is of order 1, before doing the matrix analysis, because if I don't, the large difference between the scale of q and the column of ones in the matrix causes computational inaccuracy. More specifically, qAll'*qAll will be ill-conditioned, if the columns differ by many orders of magnitude.

I didn't need to compute dt in this example, since dt=1, but it is a good idea to do so, in case dt is not =1.

% Make matrix of predictor variables
dt=(t(end)-t(1))/(length(t)-1);     % delta t
q=q*1e14;                           % multiply q so it is of order 1
dqdt=diff(q)/dt;                    % first derivative
d2qdt2=(q(3:end)-2*q(2:end-1)+q(1:end-2))/dt^2; % centered second derivative
qt=q(2:end-1);                      % q, trimmed
qAll=[d2qdt2,dqdt(1:end-1),qt,ones(length(qt),1)];
vt=v(2:end-1);                      % v, trimmed
tt=t(2:end-1);                      % t, trimmed

Next, find the vector that minimizes the sum squared error between predicted and measured v. Use the standard matrix method.

fprintf('cond(qAll^T*qAll)=%.3g.\n',cond(qAll'*qAll)); % check condition number

cond(qAll^T*qAll)=19.2.

b=(qAll'*qAll)\qAll'*vt;
L=b(1); R=b(2); C=1/b(3); v0=b(4);

Display results

fprintf('Estimated values: L=%.3g, R=%.3g; C=%.3g, v0=%.5g.\n',L,R,C,v0);

Estimated values: L=0.656, R=0.313; C=2.36, v0=99512.

Predict v using the best-fit values for L,R, C:

vtPred=qAll*b;                      % predicted voltage

Plot measured and predicted v:

figure;
plot(tt,vt,'b.',tt,vtPred,'-b')
grid on; xlabel('Time'); ylabel('v')
title('v, measured and predicted')
legend('Measured','Prediction')

Code above yields figure below:

The prediction has the correct mean value, but the time-varying part of the prediction is not close to the time-varying part of the measured signal.

I don't know where your data came from. Based on the results above, I conclude that a series L-R-C circuit is not a good model for this voltage and charge data.