[SciPy-user] nonlinear optimisation with constraints
Ernest Adrogué
I am stuck in an obnoxious optimisation problem.
Essentially, I want to find the local maximum of a
multivariate nonlinear function subject to a linear
constraint.
x = (a1, a2, a2, ..., a_n, b1, b2, b3, ..., b_n)
Maximise: f(x)
Subject to: sum(a) - sum(b) = 0
No big deal, apparently. The problem is that f(x) is defined
only when x_n > 0 for all n, as it contains lots of
log(a[i] * b[j])
which are undefined when a[i] or b[j] < 0.
I have tried to specify a lower bound for x, but both
fmin_l_bfgs_b and fmin_tnc seem to evaluate the objective
function with elements of x < 0, regardless of the bounds
specified, making my programme to crash.
fmin_cobyla seems to ignore the constraint altogether.
I have run out of ideas on how to deal with this.
Any advice?
Thanks.
Ernest
_______________________________________________
SciPy-user mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user
Sebastian Walter
2009/6/22 Ernest Adrogué <eadr...@gmx.net>:
> Hi all,
>
> I am stuck in an obnoxious optimisation problem.
>
> Essentially, I want to find the local maximum of a
> multivariate nonlinear function subject to a linear
> constraint.
>
> x = (a1, a2, a2, ..., a_n, b1, b2, b3, ..., b_n)
>
> Maximise: f(x)
> Subject to: sum(a) - sum(b) = 0
>
> No big deal, apparently. The problem is that f(x) is defined
> only when x_n > 0 for all n, as it contains lots of
>
> log(a[i] * b[j])
>
> which are undefined when a[i] or b[j] < 0.
>
> I have tried to specify a lower bound for x, but both
> fmin_l_bfgs_b and fmin_tnc seem to evaluate the objective
> function with elements of x < 0, regardless of the bounds
> specified, making my programme to crash.
>
> fmin_cobyla seems to ignore the constraint altogether.
>
> I have run out of ideas on how to deal with this.
> Any advice?
are you sure you can't reformulate the problem?
maybe you should try an interior point method. By definition, all
iterates will be feasible.
There is a python wrapper for IPOPT out there. It's called pyipopt. It
worked reasonably well when I tried it.
OPENOPT also interfaces to IPOPT as far as I know, but I have never
used that interface.
Sebastian
Ernest Adrogué
22/06/09 @ 09:57 (+0200), thus spake Sebastian Walter:
>
> are you sure you can't reformulate the problem?
Another approach would be to try to solve the system of
equations resulting from equating the gradient to zero.
Such equations are defined for all x. I have already tried
that with fsolve(), but it only seems to find the obvious,
useless solution of x=0. I was going to try with a
Newton-Raphson alorithm, but since this would require the
hessian matrix to be calculated, I'm leaving this option
as a last resort :)
> maybe you should try an interior point method. By definition, all
> iterates will be feasible.
> There is a python wrapper for IPOPT out there. It's called pyipopt. It
> worked reasonably well when I tried it.
> OPENOPT also interfaces to IPOPT as far as I know, but I have never
> used that interface.
Thanks, this looks interesting. I'm going to check out
this pyipopt.
--
Sebastian Walter
> Hi Sebastian,
>
> 22/06/09 @ 09:57 (+0200), thus spake Sebastian Walter:
>>
>> are you sure you can't reformulate the problem?
>
> Another approach would be to try to solve the system of
> equations resulting from equating the gradient to zero.
> Such equations are defined for all x. I have already tried
> that with fsolve(), but it only seems to find the obvious,
> useless solution of x=0. I was going to try with a
> Newton-Raphson alorithm, but since this would require the
> hessian matrix to be calculated, I'm leaving this option
> as a last resort :)
Ermmm, I don't quite get it. You have an NLP with linear equality
constraints and box constraints.
Of course you could write down the Lagrangian for that and define an
algorithm that satisifies the first and second order optimality
conditions.
But that is not going to be easy, even if you have the exact hessian:
you'll need some globalization strategy (linesearch, trust-region,...)
to guarantee global convergence
and implement something like projected gradients so you stay within
the box-constraints.
I guess it will be easier to use an existing algorithm...
And I just had a look at fmin_l_bfgs_b: how did you set the equality
constraints for this algorithm. It seems to me that this is an
unconstrained optimization algorithm which is worthless if you have a
constrained NLP.
remark:
To compute the Hessian you can always use an AD tool. There are
several available in Python.
My biased favourite one being pyadolc (
http://github.com/b45ch1/pyadolc ) which is slowly approaching version
1.0.
Ernest Adrogué
> 2009/6/22 Ernest Adrogué <eadr...@gmx.net>:
> > Hi Sebastian,
> >
> > 22/06/09 @ 09:57 (+0200), thus spake Sebastian Walter:
> >>
> >> are you sure you can't reformulate the problem?
> >
> > Another approach would be to try to solve the system of
> > equations resulting from equating the gradient to zero.
> > Such equations are defined for all x. I have already tried
> > that with fsolve(), but it only seems to find the obvious,
> > useless solution of x=0. I was going to try with a
> > Newton-Raphson alorithm, but since this would require the
> > hessian matrix to be calculated, I'm leaving this option
> > as a last resort :)
>
> Ermmm, I don't quite get it. You have an NLP with linear equality
> constraints and box constraints.
> Of course you could write down the Lagrangian for that and define an
> algorithm that satisifies the first and second order optimality
> conditions.
> But that is not going to be easy, even if you have the exact hessian:
> you'll need some globalization strategy (linesearch, trust-region,...)
> to guarantee global convergence
> and implement something like projected gradients so you stay within
> the box-constraints.
>
> I guess it will be easier to use an existing algorithm...
Mmmm, yes, but the box constraints are merely to prevent the
algorithm from evaluating f(x) with values of x for which f(x)
is not defined. It's not a "real" constraint, because I know
beforehand that all elements of x are > 0 at the maximum.
> And I just had a look at fmin_l_bfgs_b: how did you set the equality
> constraints for this algorithm. It seems to me that this is an
> unconstrained optimization algorithm which is worthless if you have a
> constrained NLP.
You're right. I included the equality constraint within the
function itself, so that the function I omptimised with fmin_l_bfgs_b
had one parameter less and computed the "missing" parameter
internally as a function of the others.
The problem is that this dependent parameter, was unaffected by
the box constraint and eventually would take values < 0.
Fortunately, Siddhardh Chandra has told me the solution, which
is to maximise f(|x|) instead of f(x), with the linear
constraint incorporated into the function, using a simple
unconstrained optimisation algorithm. His message hasn't made it
to the list though.
I have just done this and it seems to work! After 10.410
function evaluations and 8.904 iterations fmin has found the
solution and it looks sound at first sight.
Thanks for your help.
> remark:
> To compute the Hessian you can always use an AD tool. There are
> several available in Python.
> My biased favourite one being pyadolc (
> http://github.com/b45ch1/pyadolc ) which is slowly approaching version
> 1.0.
I will have a look.
Thanks again :)
Robert Kern
> Mmmm, yes, but the box constraints are merely to prevent the
> algorithm from evaluating f(x) with values of x for which f(x)
> is not defined. It's not a "real" constraint, because I know
> beforehand that all elements of x are > 0 at the maximum.
Unfortunately, that's not how constraints work in most optimizers.
Usually, the infeasible region is necessarily also explored.
--
Robert Kern
"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
-- Umberto Eco
Sebastian Walter
> 2009/6/22 Ernest Adrogué <eadr...@gmx.net>:
>
>> Mmmm, yes, but the box constraints are merely to prevent the
>> algorithm from evaluating f(x) with values of x for which f(x)
>> is not defined. It's not a "real" constraint, because I know
>> beforehand that all elements of x are > 0 at the maximum.
>
> Unfortunately, that's not how constraints work in most optimizers.
> Usually, the infeasible region is necessarily also explored.
Do you have a certain solver in mind?
>From what I read in the literature, I thought that for simple box
constraints usually some projection to the set of feasible search
directions is performed.
And only for the nonlinear constraints one uses a merit function and
an active set strategy for the inequality constraints which may yield
infeasible steps in the process.
Sebastian Walter
I'm curious: Could you elaborate how you incoroporated the linear
constraints into the objective function?
>
> I have just done this and it seems to work! After 10.410
> function evaluations and 8.904 iterations fmin has found the
> solution and it looks sound at first sight
>
Ernest Adrogué
There was only one constraint:
max: f(x)
s.t.: sum(a) - sum(b) = 0 (1)
where 'a' is the first half of vector x, and 'b' the second half.
What this constraint really says is that one function parameter
is dependent on the others, as you can rewrite the constraint as:
a0 = sum(b0, b1 ... b_n) - sum(a1, a2 ... a_n)
Therefore, if it is possible to incorporate the constraint into
the function, by writing a function g(x) that takes (2*n)-1
parameters, calculates the dependent parameter, and calls the
original f(x) with 2*n parameters. This g(x) function will have
the constraint "incorporated" and can be optimised with an
unconstrained algorithm.
Rob Falck
In my experience fmin_slsqp does not attempt to evaluate regions of the independent variables outside of the box bounds, but YMMV.
--
- Rob Falck