Search for Stationary Points of Constrained Nonlinear Function

I have a nonlinear function F(x1,x2). I have evaluated analytically the gradient vector of F and also the Hessian matrix. The function, its gradient and Hessian are stored in separate m-files: NUMF, NUMGRAD, and NUMHESS, respectively.

My task now is to perform a search for all stationary points of F subject to the following 3 constraints on x1 and x2 (i.e. within the region bounded by the following inequalities):

x1>=0

x2>=0

x1+x2<=1

Can someone please suggest a suitable optimization/other tool for solving the above problem? I had thought I may be able to use fsolve at various starting points to obtain stationary points one at a time, but I am not sure how and where to specify the constraints.

Thanks for your time,

Dean

[Walter Roberson]

I'm not sure what you mean by stationary point in this context? Do you mean local minimum?

I am accustomed to stationary point of a function being defined as the points that the function evaluates to the point (e.g., for x^2, the stationary points would be x=0 and x=1). In two dimensions, that would require the function to have two outputs or else output a vector with two elements.

[Bruno Luong]

What you have described Walter seems to me a fixed point. A stationary points are point where the gradient is 0 (in large sense), or local optima in a more restrict sense.

[Walter Roberson]

You are right, Bruno, I was thinking of "fixed point".

Answer 1

[Bruno Luong]

A generic tool is to solve stationary points is FMINCON. If we know more about objective function, we might orient you towards more appropriate method.

Bruno

[Dean]

The objective function is a very large, unwieldy nonlinear function. It is essentially a sum of terms in which one of the variables, eg x1 is multiplied by the logarithm of a basic function of x1 and x2 (by basic function I mean the only operations involved are +,-,* and /). It is called the tangent plane distance function (TPDF) based on the UNIFAC equation of state, a chemical engineering problem. Essentially the minima of the TPDF represent potentially stable compositions of a multicomponent mixture so it can be used to predict the compositions of phases in a partially miscible mixture of liquids. I would ideally like to find all stationary points within the constrained region, including local minima, local maxima and points of inflection, however all local minima within the region will suffice. Can FMINCON be used to find all stationary points, or just minima? If just minima, perhaps I could use FMINCON with the objective function being the norm of the gradient.

I meant to say in my question that I was considering using fsolve on the Gradient as Walter suggests. That way the main problem becomes (dF/dx1,dF/dx2)=(0,0) which is a system of two nonlinear equations in two unknowns, x1 and x2.

I don't know how to incorporate the constraints into the problem definition using fsolve, it doesn't seem to be allowed.

[Bruno Luong]

As I wrote earlier, use FMINCON. The stationary points are minimum or maximum under the constraints. The gradient are 0 if no constraint is active, and verifies KKT otherwise.

Answer 2

[Walter Roberson]

If you have the gradient analytically, then if you have the Symbolic Toolkit, you might be able to solve for gradient 0 in terms of one of the two variables (symbolically). You could then find the boundaries of the constraint region.

Not every function is suitable for analytic solutions for 0, of course, but you might find that what the toolbox is able to say about the form of the solutions might narrow the search considerably -- might turn it in to the solution of an expression that is comparatively easy to solve numerically. And if not that, then at least the symbolic toolbox might be able to isolate one of the variables, allowing you to fsolve() against an expression in one variable instead of minimizing in two variables.

Answer 3

[John D'Errico]

I would suggest a completely different approach. How do you define a stationary point? Clearly, this is where both elements of gradient are zero. So simply do TWO contour plots, overlaid on top of each other.

Look for contours of dfdx, where the contour level is zero. Color those lines red.

On top of that plot, overlaid by using hold on, plot a second set of contours perhaps in green, of the function df/dy, again where that operator is zero.

Where red and green lines intersect, this is a solution in that domain. A nice thing is, all solutions will magically appear. The accuracy of the solution is defined by the fineness of the grid used to generate the contour plots.

If you wish to actually find the list of solutions, just generate the contours temselves, using a tool like contourc. Then use a tool off the file exchange that will find intersections of piecewise linear polygonal curves. There are many such tools, but I like the one provided by Doug Schwarz:

http://www.mathworks.com/matlabcentral/fileexchange/11837-fast-and-robust-curve-intersections

Edit: I forgot that your problem is defined on a TRIANGULAR domain. This is easy enough to solve anyway, by defining the upper triangular part of the rectangle to be NaNs when you build the contour plots. Alternatively, I have tools that can work on triangulated domains.