Problem with finding the global minimum with fmincon

I am currently trying to find the global minimum for a strain-energy-function ("Holzapfel-Model") and I am running into multiple problems:

The SEF has the form

With

We can calculate , where

We want to determine the minimum of the least-square function

My solution was to put all these equations into one long one:

fun = @(x) sum((sigma_11 - (lambda_1.^2 - lambda_2.^2 .* lambda_1.^2).*x(1) + 2 .*lambda_1.^2 .*cos(x(4))^2 .* (2.*x(2).*((lambda_1.^2 .* cos(x(4))^2 + lambda_2.^2 .* sin(x(4))^2) - 1) .* exp(x(2).*((lambda_1.^2 .* cos(x(4))^2 + lambda_2.^2 .* sin(x(4))^2)-1).^2))).^2 + (sigma_22 - (lambda_2.^2 - lambda_2.^2 .* lambda_1.^2).*x(1) + 2 .*lambda_2.^2 .*sin(x(4))^2 .* (2.*x(2).*((lambda_1.^2 .* cos(x(4))^2 + lambda_2.^2 .* sin(x(4))^2) - 1) .* exp(x(2).*((lambda_1.^2 .* cos(x(4))^2 + lambda_2.^2 .* sin(x(4))^2)-1).^2))).^2)

and then use the following parameters

x0 = [15,500,12,0.75*pi];
A = [];
b = [];
Aeq = [];
beq = [];
lb = [0,0,0,0];
ub = [inf, inf, inf, pi];
chi = fmincon(fun, x0,A,b,Aeq,beq,lb,ub,nonlcon)
function [c,ceq] = nonlcon(x,lambda_1,lambda_2)
c =1-(lambda_1.^2 .* cos(x(4)).^2 + lambda_2.^2 * sin(x(4)).^2)  ;
ceq = [];
end

With these parameters, I can somewhat get close to my data points.

Now my questions:

I don't think I understood c,ceq correctly. I used c to account for the constraint on I4, but I'm not sure if this was the right way to do it.

With the initial guess for x0, I can get close but it never seems to approach my curve nearly enough. How do I know if I have a good starting guess, and is fmincon even the right approach for this problem.

I have multiple data sets, for different stretch ratios (lambda_1:lambda_2: 1-1, 1-0.75, 0.75-1, 1-0.5,0.5-1) and since they are the same sample, I would like to use those datas to get one set of parameters for all of them. I tried to put all my data into a single vector, (1:30 would be the first data set, 31:^60 the second,...). This does not seem to work well. Should I find the solution for just one curve and than try to average over the parameters? As you guys can see, I am doing this parameter evaluation thing the first time ever and I would greatly appreciate help.

[Steven Lord]

My solution was to put all these equations into one long one:

It's tempting to write your expression as one big anonymous function, but in this case with a somewhat complicated equation I'd recommend writing up your objective function as a normal function instead. Doing so will have a few benefits:

1. You can define intermediate terms using good, descriptive names. When someone else reads your code (or when you read it six months from now) the hints about the purpose of those variables given by the names could prove very useful in understanding what the code is doing. As it stands right now your problem statement has a number of variables defined that don't appear in your objective function; where do c, k1, and k2 get used in your function fun?
2. Related to the first point, you can add comments indicating the purpose of each step. This doesn't necessarily mean explaining what the code is doing (ideally that should be obvious from what functions you're calling) but more importantly means explaining why the code does what it does. Something like "We tried using approach X here but need to use approach Y instead to avoid catastrophic cancellation" is an example of a good why comment IMO.
3. You can step through the code in the debugger. If you have an example of evaluating your function by hand that you've worked out step-by-step, you can step through the code with the debugger line-by-line and confirm that each step is doing what you expect it to do.

Your format for your nonlinear constraint function is incorrect. If you look at the description of the nonlcon input on the fmincon documentation page, you'll see that it states " nonlcon is a function that accepts a vector or array x and returns two arrays, c(x) and ceq(x)." There are two problems here.

chi = fmincon(fun, x0,A,b,Aeq,beq,lb,ub,nonlcon)
function [c,ceq] = nonlcon(x,lambda_1,lambda_2)

The first is that when you call fmincon you're calling nonlcon with 0 input arguments and passing whatever that function call returns as the nonlcon input argument. That won't work. Specify it as a function handle instead, @nonlcon.

The second is that your nonlcon function accepts three inputs (and requires all three to work), not one as the documentation states. If lambda_1 and lambda_2 are additional constant parameters you want to pass into nonlcon you can do that using one of the techniques shown on this documentation page.

Answer 1

[Torsten]

What are x(1) - x(4) in your mathematical description ? In other words: What are your optimization parameters ?

You don't define lambda_1 and lambda_2 in the code you posted.

There are many errors in the computation of sigma_11 and sigma_22 in your objective function.

Your constraint formulation for I4 in nonlcon is correct, but the call to "fmincon" must be changed to

chi = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,lambda_1,lambda_2))

(I assume lambda_1 and lambda_2 are given values defined somewhere before in your code).

[Harald]

x(1)-x(4) are my parameters that I want to find.

x(1) = c, x(2) = k_1, x(3) = k_2 and x(4) = phi

sigma_11, sigma_22, lambda_1 and lambda_2 are given in my data. They are vectors. Basically lambda is my x-axis and sigma is my y-axis.

I skipped over some steps:

Thanks for the help

[Torsten]

Are you sure that for your data vectors lambda_1 and lambda_2, there exists one parameter x(4) such that lambda_1.^2 .* cos(x(4)).^2 + lambda_2.^2 * sin(x(4)).^2 >=1 for all cases ?

I4 = @(x)lambda_1.^2*cos(x(4))^2 + lambda_2.^2*sin(x(4))^2;
psi1 = @(x)0.5*x(1);
psi4 = @(x)2*x(2)*(I4(x)-1).*exp(x(3)*(I4(x)-1).^2);
sigma_11_sim = @(x)(2*lambda_1.^2-lambda_1.^(-2).*lambda_2.^(-2)).*psi1(x)+2*lambda_1.^2.*cos(x(4))^2.*psi4(x);
sigma_22_sim = @(x)(2*lambda_2.^2-lambda_1.^(-2).*lambda_2.^(-2)).*psi1(x)+2*lambda_2.^2.*sin(x(4))^2.*psi4(x);
fun = @(x)sum((sigma_11-sigma_11_sim(x)).^2 + (sigma_22-sigma_22_sim(x)).^2);
x0 = [15,500,12,0.75*pi];
lb = [0,0,0,0];
ub = [inf, inf, inf, pi];
chi = fmincon(fun,x0,[],[],[],[],lb,ub,@(x)nonlcon(x,lambda_1,lambda_2))

function [c,ceq] = nonlcon(x,lambda_1,lambda_2)
    c = 1 - (lambda_1.^2*cos(x(4))^2 + lambda_2.^2*sin(x(4))^2)  ;
    ceq = [];
end

[Harald]

That is an excellent point. x(4) describes the angle between the fibers of the material. I used the paper of "Holzapfel GA. Nonlinear solid mechanics. A continuum approach for engineering. Chichester: John Wiley & Sons; 2000" and since I don't know the angle, I chose to set it as an phenomenological parameter. Theoretically, with higher strain, the fiber angle might change and I didn't account for that. I am not sure on how to change my equations now, but thank you for your help.