Someone please help to rectify this error: Error using barrier Objective function is undefined at initial point. Fmincon cannot continue. Optimization problem

I want to fit the mathematical model containing double numerical integration with the data. Please help me!

This is the driver code:

clc;
clear all;
u = load("thoracicdata11.mat");
xdata11 = u.Figure9S11(1:225,1);
ydata11 = u.Figure9S11(1:225,2);
x0=[5 2 1 0.01 0.1];%
lb=[0.001 0.001 0.001 1 0.0];
ub=[100.0 200.0 100.0 100 pi/2];
options = optimset('Display', 'iter', 'MaxIter', 100, 'TolFun', 1e-6);
obj = obj_GOHGR_planar(x0,xdata11,ydata11)

Warning: Non-finite result. The integration was unsuccessful. Singularity likely.

Warning: Non-finite result. The integration was unsuccessful. Singularity likely.

Index in position 2 exceeds array bounds. Index must not exceed 1.

Error in solution>obj_GOHGR_planar (line 19)
obj = sum((ydata11(:,2)-BiaxialResponse).^2);

%[x, fval]=fmincon(@(x)obj_GOHGR_planar(x,xdata11,ydata11),x0,[],[],[],[],lb,ub,[],options);

This is the objective function:

function obj=obj_GOHGR_planar(x,xdata11,ydata11)
%xdata11 = Figure9S11(1:225,1);
%ydata11 = Figure9S11(1:225,2);
BiaxialResponse=zeros(size(xdata11(:,1)));
for i=1:size(xdata11(:,1))
  BiaxialResponse(i)= BiaxialNewInvariantFunc11(x,xdata11(i,1));
end
obj = sum((ydata11(:,2)-BiaxialResponse).^2);
end

Calling function:

function K1 = BiaxialNewInvariantFunc11(par,lambda1)
  integrand=@(par,lambda1,theta,phi) (3/4).*exp(1).^((1/64).*lambda1.^(-12).*(1+(-2).*lambda1.^4+ ...
  lambda1.^6+(-1).*((-1)+lambda1.^6).*cos(2.*theta)+ ...
  lambda1.^4.*(lambda1.^(-8).*(2.*((-1)+lambda1.^2).^2.*(1+2.* ...
  lambda1.^2+lambda1.^4+lambda1.^8)+(-2).*((-1)+2.*lambda1.^4+ ...
  (-2).*lambda1.^10+lambda1.^12).*cos(2.*theta))).^(1/2)).^3.* ...
  par(4)+2.*cos(phi+(-1).*par(5)).^2.*par(2).*sin( ...
  theta).^2).*lambda1.^(-8).*(2.*pi).^(-1/2).*(1+(-2).* ...
  lambda1.^4+lambda1.^6+(-1).*((-1)+lambda1.^6).*cos(2.*theta) ...
  +lambda1.^4.*(lambda1.^(-8).*(2.*((-1)+lambda1.^2).^2.*(1+ ...
  2.*lambda1.^2+lambda1.^4+lambda1.^8)+(-2).*((-1)+2.* ...
  lambda1.^4+(-2).*lambda1.^10+lambda1.^12).*cos(2.*theta))) ...
  .^(1/2)).^2.*erfi(2.^(1/2).*par(2).^(1/2)).^(-1).*par( ...
  1).*par(2).^(1/2).*par(3).*sin(theta).*((1/2).* ...
  lambda1.^2.*cos(phi).^2.*sin(theta).^2.*(1+((-1)+lambda1.^2) ...
  .*(lambda1.^(-8).*((-1)+lambda1.^2).^2.*((1+lambda1.^2).^2.* ...
  cos(theta).^2+lambda1.^8.*sin(theta).^2)).^(-1/2))+(-1/2).* ...
  lambda1.^(-4).*cos(theta).^2.*(1+lambda1.^(-4).*(1+(-1).* ...
  lambda1.^4).*(lambda1.^(-8).*((-1)+lambda1.^2).^2.*((1+ ...
  lambda1.^2).^2.*cos(theta).^2+lambda1.^8.*sin(theta).^2)).^( ...
  -1/2)));
  K1 = 0.3E0.*((-1).*lambda1.^(-4)+lambda1.^2)+ integral2(@(phi,theta) integrand(par,lambda1,theta,phi),0,2*pi,0,pi);
 end 

Thanks in advance! :)

[Torsten]

Besides the unsuccessful integration one syntax error: ydata11(:,2) does not exist because ydata11 has only one column, not two.

Answer 1

[Steven Lord]

We can't run your code as we don't have all of the functions it uses, but what do you receive as an answer when you call obj_GOHGR_planar using the x0 you pass into fmincon as the input? I've commented out the call below because I want to run later code in this answer, but I wanted to show what I want you to run.

%{
y = obj_GOHGR_planar(x0)
%}

But I see another problem with your code that may be related. Does your initial point satisfy your lower and upper bounds?

x0=[5 2 1 0.01 0.1];%
lb=[0.001 0.001 0.001 1 0.0];
ub=[100.0 200.0 100.0 100 pi/2];
isx0InBounds = (lb <= x0) & (x0 <= ub)

isx0InBounds = 1×5 logical array
   1   1   1   0   1

Are you sure the fourth element of x0 should be 0.01?