Error: Custom function fitting in MATLAB

Hey guys, I would like to know why I got this error and how can I fix it. Can I please have your guys comments on this code:

% Load the experimental data and assgin the data to each parameter
data2fit2 = xlsread('data2fit2.xlsx');
ST = data2fit2(:,1);
theta_bar = data2fit2(:,2);
f_exp = data2fit2(:,3);

% Initial Guess
par0 = [0.417257,0.04132,0.001177,0.00961,1.946234];
pars = zeros(size(par0));

% Get the estimate parameters to have best fit to measurement data
par = fminsearch(@(par) err_fcn(pars,...
                                data2fit2,...
                                ones(size(data2fit2)),...
                                @(s,t,p) epsilon_f_of_sigma(s,t,p),...
                                ST,...
                                theta_bar),...
                 par0);

Arrays have incompatible sizes for this operation.

Error in solution>err_fcn (line 63)
err = sum((data2fit2(:) - model_data(:)).^2.*weights4scaling);

Error in solution>@(par)err_fcn(pars,data2fit2,ones(size(data2fit2)),@(s,t,p)epsilon_f_of_sigma(s,t,p),ST,theta_bar) (line 12)
par = fminsearch(@(par) err_fcn(pars,...

Error in fminsearch (line 216)
fv(:,1) = funfcn(x,varargin{:});

% Use optimized parameters to calculate the estimated epsilon_f by HC model
epsilon_f_modeled = epsilon_f_of_sigma(ST, theta_bar, par);

% Plot fracture criterion
figure;
plot(ST, f_exp, 'o', ST, epsilon_f_modeled, '-')
xlabel('Stress Triaxiality [-]')
% ylabel('Lode Angle Parameter [-]')
ylabel('Equivalent Plastic Strain at Fracture [-]')

function epsilon_f = epsilon_f_of_sigma(ST, theta_bar, pars)
    b0 = pars(1);
    gamma = pars(2);
    c = pars(3);
    epsilon_p = 0.0005;
    epsilon_0 = 0.0005;
    n = pars(4);
    a = pars(5); 
    
    % Determine strain rate effect on parameter b
    if epsilon_p < epsilon_0
      b = b0;
    else
      b = b0*(1 +gamma*log(epsilon_p/epsilon_0));
    end
    
    % Define Lode-dependent functions
    f1 = (2/3)*cos((pi/6)*(1-theta_bar));
    f2 = (2/3)*cos((pi/6)*(3+theta_bar));
    f3 = -(2/3)*cos((pi/6)*(1+theta_bar));

    % Define each term
    Strain_rate_term = b*(1+c).^(1./n);
    Lode_dependent_term = ((1/2).*((f1-f2).^a+(f2-f3).^a+(f1-f3).^a)).^(1./a);
    Triaxiality_term = c.*(2*ST+f1+f3);
    
    % Final expression for HC model
    epsilon_f = Strain_rate_term.* (Lode_dependent_term + Triaxiality_term).^(-1./n);
end

% Define error function
function err = err_fcn(pars,data2fit2,weights4scaling,model_fcn,ST,theta_bar)
    model_data = model_fcn(ST,theta_bar,pars);
    err = sum((data2fit2(:) - model_data(:)).^2.*weights4scaling);
end

Answer 1

[Cris LaPierre]

data2fit2(:) and model_data(:) and weights4scaling do not all have the same size.

A = rand(3)

A = 3×3

0.6170 0.3529 0.7404 0.5912 0.8156 0.0387 0.5109 0.5150 0.8625

B = rand(2,1)

B = 2×1

0.7242 0.0795

% Your error
A-B

Arrays have incompatible sizes for this operation.

[LM]

It works now! Thanks for your help, Cris!

Can I please have your comments:

% Load the experimental data and assgin the data to each parameter
data2fit2 = xlsread('data2fit2.xlsx');
ST = linspace(-0.3,2,24);
theta_bar = linspace(-1,1,24);
f_exp = data2fit2(:,3);

sigma = data2fit2(:,1);
theta_bar1 = data2fit2(:,2);

% Initial Guess
par0 = [0.417257,0.043132,0.001177,0.002961,1.946234];
pars = ones(size(par0));

% Get the estimate parameters to have best fit to measurement data
par = fminsearch(@(par) err_fcn(pars,...
                                f_exp,...
                                ones(size(f_exp)),...
                                @(s,t,p) epsilon_f_of_sigma(s,t,p),...
                                ST,...
                                theta_bar),...
                 par0);

% Use optimized parameters to calculate the estimated epsilon_f by HC model
epsilon_f_modeled = epsilon_f_of_sigma(ST, theta_bar, par);

% Plot fracture criterion
figure;
plot(sigma, f_exp, 'o', linspace(-0.3,2,24), epsilon_f_modeled, '-')
xlabel('Stress Triaxiality [-]')
% ylabel('Lode Angle Parameter [-]')
ylabel('Equivalent Plastic Strain at Fracture [-]')

function epsilon_f = epsilon_f_of_sigma(ST, theta_bar, pars)
    b0 = pars(1);
    gamma = pars(2);
    c = pars(3);
    epsilon_p = 0.0005;
    epsilon_0 = 0.0005;
    n = pars(4);
    a = pars(5); 
    
    % Determine strain rate effect on parameter b
    if epsilon_p < epsilon_0
      b = b0;
    else
      b = b0*(1 + gamma*log(epsilon_p/epsilon_0));
    end
    
    % Define Lode-dependent functions
    f1 = (2/3)*cos((pi/6)*(1-theta_bar));
    f2 = (2/3)*cos((pi/6)*(3+theta_bar));
    f3 = -(2/3)*cos((pi/6)*(1+theta_bar));

    % Define each term
    Strain_rate_term = b*(1+c).^(1./n);
    Lode_dependent_term = ((1/2).*((f1-f2).^a+(f2-f3).^a+(f1-f3).^a)).^(1./a);
    Triaxiality_term = c.*(2*ST+f1+f3);
    
    % Final expression for HC model
    epsilon_f = Strain_rate_term.* (Lode_dependent_term + Triaxiality_term).^(-1./n);
end

% Define error function
function err = err_fcn(pars,f_exp,weights4scaling,model_fcn,ST,theta_bar)
    model_data = model_fcn(ST,theta_bar,pars);
    err = sum((f_exp(:) - model_data(:)).^2.*weights4scaling);
end

[Cris LaPierre]

i see that you have asked this question here. I agree with what @Torsten said there - you should be fitting to your experimental data. Once you have that, you can use the results to then generate the fracture surface.

I think there is one thing to be aware of. After doing some digging in, it seems the equations you have shared are for generating the surface. However, the curves you are looking for seem to be generated under a specific condition, as illustrated here.

The curve you are currently plotting is the line along the surface that runs diagonally across the surface from the min Lode Angle and min Stress Triaxiality to the max values. However, the curve you are trying to mimic looks like those in the bottom plot shared above, which appear to follow the Plane Stress line shown in the surface plots.

That's about what I've been able to figure out tonight. Since this is your area of expertise, I leave it to you to provide us with the equation for solving the integral.

[LM]

Hi Cris, I literally confused right now. Would you mind to give me some tips on this?

- For you mentioned equation, which is for caculating the damage parameter and I think it will not be involved in the fitting equation.

And yes, I agree with you guys on fitting on the measurement data, and I just don't know which step I made wrong. My apologies to this.

[Cris LaPierre]

No, the new equation does not play a role in obtaining the fit parameters.

Perhaps you are trying to do too much at once. Simplify your code down to just fitting the equation to your data. Get rid of all your extra variables and unnecessary arrays of ones.