Hi all,

This originally shouldn't be a problem -- but it is now our most
annoying headache.

As many other people reported(we found many other posts talking about
this), Matlab's SQP based fmincon has a weird behavior, it needs to
first evaluate the functions slightly outside the boundaries, and then
get back in the boundaries.

But our model has log's and sqrt's, which generates complex, infinite
or NaN (not a number) if boundaries are violated.

As a simple example, let's say we have log(a(x, y, z)),

and we require a(x, y, z)>0,

and we put -a(x, y, z)+myeps <=0, because the solver can only handle
<= constraints.

we set myeps=1e-8.

Now recall that the solver will compute function values that slightly
violate the boundaries.

And the objective function f(x, y, z) is evaluated before the
constraint function c(x, y, z).

Q1.

When the objective function f(x, y, z) is evaluated at the points that
are slightly beyond the boundaries, it returns
complex, infinite or NaN (not a number) . What shall we do in this
case?

The Matlab Optimization Toolbox manual says: "you should check for
these cases, possible solutions are to (1) place bounds on variables;
(2) make a penalty function to give large positive number to both the
objective function f(x, y, z) and constraint function c(x, y, z). The
penalty function should be smooth and continuous."

But (1) we've already placed bounds on variables, it still gets
violated; (2) I  cannot think of a penalty function that is smooth and
continuous, when the complex, infinite, or NaN are detected, they are
already out of numeric domain, how to make a smooth and continuous
penalty function out of them?

Q2.

In the constraint function c(x, y, z):

The first constraint is:

-a(x, y, z)+myeps <=0;

Let's suppose the second constraint is:

b(log(a(x, y, z))) <=0;

As you can see, the second constraint relies on a(x, y, z),

and when the first constraint is violated,

the constraint function doesn't stop, it continues to evaluate
b(log(a(x,y, z))) and gives weird values.

This makes the constraints invalid themselves.

Again, the only way out might be impose a penalty function. But the
constraint function also has to be continuous and smooth, how to make
such nice a penalty function?

--------------------------
Any thoughts? We are stuck! Please shed some lights!

If you have better solvers that don't "slightly go beyond" the
boundaries like SQP does, please let me know.

Thanks a lot!

Luna Moon wrote:
> Hi all,
> 
> This originally shouldn't be a problem -- but it is now our most
> annoying headache.
> 
> As many other people reported(we found many other posts talking about
> this), Matlab's SQP based fmincon has a weird behavior, it needs to
> first evaluate the functions slightly outside the boundaries, and then
> get back in the boundaries.
> 

You're right, that's weird (and disconcerting).

> But our model has log's and sqrt's, which generates complex, infinite
> or NaN (not a number) if boundaries are violated.
> 
> As a simple example, let's say we have log(a(x, y, z)),
> 
> and we require a(x, y, z)>0,
> 
> and we put -a(x, y, z)+myeps <=0, because the solver can only handle
> <= constraints.
> 
> we set myeps=1e-8.
> 
> Now recall that the solver will compute function values that slightly
> violate the boundaries.
> 
> And the objective function f(x, y, z) is evaluated before the
> constraint function c(x, y, z).
> 
> Q1.
> 
> When the objective function f(x, y, z) is evaluated at the points that
> are slightly beyond the boundaries, it returns
> complex, infinite or NaN (not a number) . What shall we do in this
> case?

I'm not a Matlab user, but I assume that Matlab has an if-else construct 
of some kind.  Using your example of f() = log(a()), I'd try something 
like f(x, y, z) = log(a(x, y, z)) if a(x, y, z) > myeps and f(x, y, z) = 
10^20 if not.
> 
> The Matlab Optimization Toolbox manual says: "you should check for
> these cases, possible solutions are to (1) place bounds on variables;
> (2) make a penalty function to give large positive number to both the
> objective function f(x, y, z) and constraint function c(x, y, z). The
> penalty function should be smooth and continuous."
> 
> But (1) we've already placed bounds on variables, it still gets
> violated; (2) I  cannot think of a penalty function that is smooth and
> continuous, when the complex, infinite, or NaN are detected, they are
> already out of numeric domain, how to make a smooth and continuous
> penalty function out of them?

Sticking with the same example, you already basically have a penalty 
function -- log(a(x, y, z)) blows up smoothly as a(x, y, z) -> 0.  A 
penalty function won't work if the algorithm starts close enough to the 
boundary and then takes a large enough step to jump over the boundary, 
which I gather is happening here.  Possible cures include increasing 
myeps and/or shrinking the step size fmincon is using when it pokes 
around outside the search interval (if you can).
> 
> Q2.
> 
> In the constraint function c(x, y, z):
> 
> The first constraint is:
> 
> -a(x, y, z)+myeps <=0;
> 
> Let's suppose the second constraint is:
> 
> b(log(a(x, y, z))) <=0;
> 
> As you can see, the second constraint relies on a(x, y, z),
> 
> and when the first constraint is violated,
> 
> the constraint function doesn't stop, it continues to evaluate
> b(log(a(x,y, z))) and gives weird values.
> 
> This makes the constraints invalid themselves.
> 
> Again, the only way out might be impose a penalty function. But the
> constraint function also has to be continuous and smooth, how to make
> such nice a penalty function?

My suggestion would be to redefine a() so that it retains its original 
definition at all points where a(x, y, z) >= myeps, but for (x, y, z) 
outside that domain the new a() asymptotically approaches myeps/2.

For instance, define z(t) = (1 - cos(pi t))/2.  Note that z is smooth 
with z(0) = z'(0) = z'(1) = 0 and z(1) = 1.  Now define a*(x, y, z) as 
follows:

a*(x, y, z) = a(x, y, z) if a(x, y, z) >= myeps;

a*(x, y, z) = z(t)*a(x, y, z) + (1 - z(t))*myeps/2 if myeps/2 <= a(x, y, 
z) < myeps, where t = 2 a(x, y, z)/myeps - 1;

a*(x, y, z) = myeps/2 if a(x, y, z) < myeps.

a* is smooth, bounded below by myeps/2, and matches a in the feasible 
region (where a >= myeps).

/Paul

On Aug 5, 2:01 pm, "Paul A. Rubin" <ru...@msu.edu> wrote:
> Luna Moon wrote:
> > Hi all,
>
> > This originally shouldn't be a problem -- but it is now our most
> > annoying headache.
>
> > As many other people reported(we found many other posts talking about
> > this), Matlab's SQP based fmincon has a weird behavior, it needs to
> > first evaluate the functions slightly outside the boundaries, and then
> > get back in the boundaries.
>
> You're right, that's weird (and disconcerting).
>
>
>
> > But our model has log's and sqrt's, which generates complex, infinite
> > or NaN (not a number) if boundaries are violated.
>
> > As a simple example, let's say we have log(a(x, y, z)),
>
> > and we require a(x, y, z)>0,
>
> > and we put -a(x, y, z)+myeps <=0, because the solver can only handle
> > <= constraints.
>
> > we set myeps=1e-8.
>
> > Now recall that the solver will compute function values that slightly
> > violate the boundaries.
>
> > And the objective function f(x, y, z) is evaluated before the
> > constraint function c(x, y, z).
>
> > Q1.
>
> > When the objective function f(x, y, z) is evaluated at the points that
> > are slightly beyond the boundaries, it returns
> > complex, infinite or NaN (not a number) . What shall we do in this
> > case?
>
> I'm not a Matlab user, but I assume that Matlab has an if-else construct
> of some kind.  Using your example of f() = log(a()), I'd try something
> like f(x, y, z) = log(a(x, y, z)) if a(x, y, z) > myeps and f(x, y, z) =
> 10^20 if not.
>
>
>
> > The Matlab Optimization Toolbox manual says: "you should check for
> > these cases, possible solutions are to (1) place bounds on variables;
> > (2) make a penalty function to give large positive number to both the
> > objective function f(x, y, z) and constraint function c(x, y, z). The
> > penalty function should be smooth and continuous."
>
> > But (1) we've already placed bounds on variables, it still gets
> > violated; (2) I  cannot think of a penalty function that is smooth and
> > continuous, when the complex, infinite, or NaN are detected, they are
> > already out of numeric domain, how to make a smooth and continuous
> > penalty function out of them?
>
> Sticking with the same example, you already basically have a penalty
> function -- log(a(x, y, z)) blows up smoothly as a(x, y, z) -> 0.  A
> penalty function won't work if the algorithm starts close enough to the
> boundary and then takes a large enough step to jump over the boundary,
> which I gather is happening here.  Possible cures include increasing
> myeps and/or shrinking the step size fmincon is using when it pokes
> around outside the search interval (if you can).
>
>
>
>
>
> > Q2.
>
> > In the constraint function c(x, y, z):
>
> > The first constraint is:
>
> > -a(x, y, z)+myeps <=0;
>
> > Let's suppose the second constraint is:
>
> > b(log(a(x, y, z))) <=0;
>
> > As you can see, the second constraint relies on a(x, y, z),
>
> > and when the first constraint is violated,
>
> > the constraint function doesn't stop, it continues to evaluate
> > b(log(a(x,y, z))) and gives weird values.
>
> > This makes the constraints invalid themselves.
>
> > Again, the only way out might be impose a penalty function. But the
> > constraint function also has to be continuous and smooth, how to make
> > such nice a penalty function?
>
> My suggestion would be to redefine a() so that it retains its original
> definition at all points where a(x, y, z) >= myeps, but for (x, y, z)
> outside that domain the new a() asymptotically approaches myeps/2.
>
> For instance, define z(t) = (1 - cos(pi t))/2.  Note that z is smooth
> with z(0) = z'(0) = z'(1) = 0 and z(1) = 1.  Now define a*(x, y, z) as
> follows:
>
> a*(x, y, z) = a(x, y, z) if a(x, y, z) >= myeps;
>
> a*(x, y, z) = z(t)*a(x, y, z) + (1 - z(t))*myeps/2 if myeps/2 <= a(x, y,
> z) < myeps, where t = 2 a(x, y, z)/myeps - 1;
>
> a*(x, y, z) = myeps/2 if a(x, y, z) < myeps.
>
> a* is smooth, bounded below by myeps/2, and matches a in the feasible
> region (where a >= myeps).
>
> /Paul

Thanks Paul. I will have to digest and think about what you suggested
carefully and will keep you updated! Appreciate your help!