Can I use multiple values for lb and ub in fmincon?

I'm trying to minimize a function, and part of that function involves a sum (lsqcurvefit would not apply in this case). It looks as follows:

function [out] = sm_signal_min(params)

[F, par, perp, signal, bvalues] = deal(params(1,1),params(1,2),params(1,3),params(:,4),params(:,5));
D = 3.00e-3;

out = 0;
for b = 1:length(bvalues)
    out = out + (log((signal(b) - (1-F).*exp(-bvalues(b).*D))./F) + bvalues(b).*perp + ...
        log((sqrt(4*bvalues(b).*(par - perp)./pi))./erf(sqrt(bvalues(b).*(par - perp))))).^2;
end
end

function [ineq, eq] = constraints(params)

[perp, par] = deal(params(2),params(2)+params(3));
ineq = perp - par;
eq = [];
end

And the call to fmincon looks like:

params0 = zeros(b_num_unique, 5);
params0(1, 1) = 1; % non-free water fraction
params0(1, 2) = 2.1e-4; % initial lambda_perp guess
params0(1, 3) = 2.1e-3; % initial deltaLambda guess
for b = 1:b_num_unique
    params0(b, 4) = S_avg(x, y, z, b);
    params0(b, 5) = b_unique(b);
end
                
params = fmincon(@sm_signal_min, params0, [], [], [], [],...
         [0,0,0, params0(b, 4), params0(b, 5)], [1, 3.0e-3, 3.0e-3, params0(b, 4), params0(b, 5)],...
         @constraints, opts);

When using this function as a part of fmincon, I receive the error "Error using erf. Input must be real and full", which is because bvalues(b) is also being fit when it meant to be held constant. How do I hold params(:,4) and params(:,5) constant within my fit? I tried including that as a part of the bounds of fmincon, but it appears that lb and ub must include only constants.

Thank you!

Warren

[Matt J]

ineq = perp - par;

This constraint is equivalent to params(3)>=0, which is both not nonlinear and also already expressed in your lb, ub bound constraint arguments. So there is no need to hae the constraint() function.

Answer 1

[John D'Errico]

If you want to have one of the parameters as a bound, then it is NOT a bound constraint!!!!!

For example, suppose I wanted to solve the problem where I would minimize a function over a triangular region? That is, consider the region:

P = polyshape([0 0;1 1;0 1]);
plot(P)

You might think that you could set up the bounds as for x as [0,1], and then the lower bound for y is x, with the upper bound being constant at 1.

And while that works if you are doing an integration, that is not what lower and upper bounds means in an optimization. These codes need CONSTANT values for the bounds.

As such, the "bound" on y must be implemented as a linear inequality constraint, that

y >= x

I'm sorry. But you cannot do what you want without using an inequality constraint here. (Could you do it using a mathematical technique called conformal mapping? Well, technically yes, by employing a mathematical trick that would subtly change your problem, while also potentially creating a subtle problem in parts of the domain. Is it probably simpler to just use an inequality constraint? Yes.)

Answer 2

[Walter Roberson]

erf(sqrt(bvalues(b).*(par - perp)))

erf requires its inputs to be real-valued. The result of a sqrt() will not be real-valued if the expression is negative. That expression could be negative if bvalues(b) or par-perp are negative.

You construct the upper and lower bounds for your fourth and fifth inputs to be the same. Are they negative? We do not know what S_avg or b_unique do, so we cannot tell.

How many parameters are there? There are b_num_unique by 5 of them. They are numbered in linear indexing for ub and lb purposes, and we do not know b_num_unique so we do not know which entries are being constrained. You are providing 5 lb and ub but those only correspond to a row of 5 values if b_num_unique happens to be 1. Entries beyond the 5th would be given lb -inf and ub +inf

For the purpose of further analysis let us imagine that b_num_unique is 1 after all, and trace the deal() to match the calculated constraints to par and perp.

[F, par, perp, signal, bvalues] = deal(params(1,1),params(1,2),params(1,3),params(:,4),params(:,5));

so par is coming from params(1,2) and perp is coming from params(1,3) . Your bounds are

[0,0,0, params0(b, 4), params0(b, 5)], [1, 3.0e-3, 3.0e-3, params0(b, 4), params0(b, 5)]

the second and third of which are constrained between 0 and 3.0e-3 . However you are only initializing 5 bound constraints which is not the number you need if b_num_unique does not happen to be 1.

You have

function [ineq, eq] = constraints(params)
[perp, par] = deal(params(2),params(2)+params(3));

which pulls out the 2nd and 3rd linear entries from params -- which is the same as params(1,2) and params(1,3) if b_num_unique happens to be 1, but is otherwise different.

But are perp and par in that function the same as par and perl in the sm_signal_min function? NO.

1. perp comes from params(2) in the constraints function but from params(3) in sm_signal_min function
2. par comes from params(2)+params(3) in the constraints function but from params(2) alone in sm_signal_min function. Even if you had not flipped the order, here you make one of the two a sum rather than the plain parameter
3. ineq = perp - par; -- why bother to form the sum and then subtract off one of the two values that went into the sum ??

Because of the way you form a sum and then subtract, the ineq is constraining only one of the two param values, not constraining par - perp for the terms of sqrt(bvalues(b).*(par - perp))

If you need constrain par - perp to be positive, then do not use nonlinear constraints for that purpose: use linear inequality constraints, such as

A = [0 -1 1 0 0]; b = 0;

This means that -par + perp <= 0 which is par - perp >= 0

Note that if your b_num_unique is guaranteed to be 1, then do not use it in the code, hard-code it instead.

If it might be something other than 1, then you need to construct your constraints properly and make sure you are extracting the correct values from the arrays.