Minimum enclosing ellipse

Minimum encompassing circle and ellipse

This is again about finding the smallest geometric shape containing all our data points. Here I focus on easy, convex cases: a circle and an ellipse. After the last post, where I was struggling mightily with a non-convex version of this model, this should be a breeze. First let's do circles.

Revisiting a crazy global NLP problem

Minimum encompassing triangle

This looks like a simple problem. Given \(n\) 2d points, find the smallest encompassing triangle. I follow the formulation from [1]. This post is self-contained, so you don't have to go back to [1]. The main contribution here is how one could work around the performance issues of the original formulation. 

Summary

The problem can be formulated as a nonconvex quadratic problem. However, this leads to models that are very difficult for global NLP solvers: global optimality can not be proven in reasonable time, even for very small data sets. One possible way to attack the problem is to use a different measure for size of the triangle: circumference instead of the more obvious area of a triangle. Although there are explicit algorithms for this problem [3,4], here we are trying to use an optimization model.

Experience with NLP solvers on a simple economic growth model

Growth Models

In this post, we formulate and solve a relatively simple economic growth model where we want to find an optimal savings rate. Consumers can either use their income for direct consumption or they can save. Savings lead to investment and to more income down the road. The resulting NLP (nonlinear programming) model is not too complicated, and it is interesting to see how we can implement this model using different tools.

I often get the question: why not use Python or R for my (economic) models? This note tries to give an answer using a rather smallish model. Although the model is small and simple, it actually illustrates some of the problems we can face when using say Python or R code. Some of the code I present here is not very obvious, or self-explanatory. One reason is that some algorithm developers look at NLP problems as something like: \[\begin{align}\min\>&f(x) \\ &g(x)=0 \\& h(x)\ge 0\end{align}\] Unfortunately, such a representation is close to the solver, but far away from practical models. It is just the wrong abstraction level. The model implementation impedance becomes even larger if we need to provide analytic gradients. Time we could have spent on developing better models is instead wasted on error-prone low-level detailed programming. This problem is even more costly, when we think about maintenance such as bug-fixing and adding new features to the model. We'll see how this looks like for this model.  

The underlying economic model is often called a Ramsey growth model, after Frank Plumpton Ramsey [1].

Crack the passcode

In this puzzle [1,2], we need to determine what the 3-digit passcode is, using a few hints. Each digit is an integer between 0 and 9. The hints are:

Let's see if we can shoehorn this into a MIP model.

Sorting: minimize number of swaps

In this post, I want to delve further into sorting. In a question on or.stackexchange.com [1], the subject was minimizing the number of swaps in sorting algorithms. A swap, i.e., an interchange of two items, is a basic operation in sorting. We usually don't pay much attention to this. First, we assume swaps are cheap. If they are not, we can sort not the real data (which can be complex, and thus more complicated and expensive to swap), but rather pointers to the data or keys. But what if we really have to reorder bulky physical things? Then the number of swaps is a more meaningful concept. Can we minimize the work to reorder say \(n\) boxes (or containers to make it more dramatic) by minimizing the number of swaps (assuming interchanging the position of two boxes is the way reordering works).

Clock problem

From [1]:

The hour, minute and second hands of this clock are all the same length and move smoothly in a circle. The dial contains hour and minute markers, but the numbers are missing. Therefore, it’s impossible to tell which one of the 12 hour markers belongs to the 12. The two hands on the left are positioned exactly on hour markers, and the hand on the right is positioned between a minute and an hour marker. What time does the clock show?

It is possible to solve this without really any math, but, of course, here I try to model this as a mathematical programming model.

diag(x)

When using \({\bf diag}(x)\) in a text, there is always the nagging feeling that there must be a nice way to express this in standard matrix algebra (i.e., some combination of identity matrices, all-ones vectors, and standard matrix multiplications). To remind ourselves, the \({\bf diag}(x)\) function creates a diagonal matrix with \(x_i\) as diagonal elements: \[{\bf diag}(x) = \begin{bmatrix} x_1 & & & \\ & x_2 & & \\ & & \ddots & \\ & & & x_n \end{bmatrix} \] Obviously, my intuition is wrong here, as I have never seen such formulation, and \({\bf diag}(x)\) is used all over the place. I have seen some attempts, but they require exotic notation, and, as a result, don't really improve upon straight use of \({\bf diag}(x)\).

Graph connectivity as constraints

I was generating, for an example model, a random, sparse, directed graph. Unfortunately, when sparse enough, this is likely to yield a graph that is not connected. Here is my first attempt.

Nodes

I generate simply \(n=25\) nodes with \((x,y)\) coordinates from a uniform distribution \(U(0,100)\). This looks like:

\(n=25\) nodes randomly placed

Revisiting a continuous facility location problem

I am revisiting here a problem from [1]:

We have \(n\) demand points and their locations. How many facilities do we need to service these customers? And where do we place them? We have a restriction: there is a maximum distance between customer and facility.   

Data

We randomly generated 75 demand points inside the \([0,1]\times[0,1]\) square. The maximum allowed distance between a demand point and a facility is 0.25. 

Random Demand Points

Nonconvex problem: local vs multistart vs global

In [1] a somewhat abstract non-convex problem is given:

\[\begin{align}\min_x & - x_1^2 - x_2^2 - x_3^2 - x_4^2\\ & Ax \le b \end{align}\]

This is a nonconvex objective with some linear constraints. Of course, the easiest way to gain some insight into how to solve this is to perform some quick and dirty tests with a global QP solver designed for these types of models. As the constraints are linear, a solver like Cplex or Gurobi may be a good starting point (pun intended).

If you have access to a local NLP solver, it may not always be easy to find a good starting point. One possible approach is to use a multistart algorithm (some NLP solvers, like Baron and Knitro, have this built-in). This will not guarantee a global solution (and most likely, it doesn't give you one), but at least we can prevent really bad solutions.