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

Equity in optimization models

In optimization models, we often use an aggregate measure in the objective function, such as total profit, the sum of tardiness of jobs, and countrywide GDP. This can lead to particularly bad results for some individuals or groups.

Here is an example I have used on several occasions. 

Problem Statement

We have \(P\) persons. They must be assigned to \(M\) groups or teams. For simplicity, we can assume \(n\) is a multiple of \(m\), and the group size is \[\frac{N}{M}\] Each person \(p_1\) specifies some preferences to be placed in the same group as a person \(p_2\). A negative preference can be used to indicate that I prefer to be in a different group. Find an optimal assignment taking into account these preferences. 

String Art

 

In [1], a greyscale picture is approximated by strings (lines) between points around the image. Here, I will try something similar with a formal optimization model.

Portfolio optimization with risk as constraint?

The standard mean-variance portfolio optimization models have the form:

Model M1
\[\begin{align}\min\>&\color{darkred}{\mathit{Risk}}=\color{darkred}x^T\color{darkblue}Q\color{darkred}x\\ & \color{darkblue}r^T\color{darkred}x \ge \color{darkblue}{\mathit{MinimumReturn}} \\  &\sum_i \color{darkred}x_i = 1\\ & \color{darkred}x \ge 0\end{align}\]

or may be:

Model M2
\[\begin{align}\min\>&\color{darkred}z=\color{darkred}x^T\color{darkblue}Q\color{darkred}x - \color{darkblue}\lambda \cdot\color{darkblue}r^T\color{darkred}x\\ &\sum_i \color{darkred}x_i = 1\\ & \color{darkred}x \ge 0\end{align}\]

where

• \(\color{darkblue}Q\) is a variance-covariance matrix (we assume here it is positive semi-definite [1]),
• \(\color{darkblue}\lambda\) is an exogenous constant that sometimes is varied to draw an efficient frontier, and
• \(\color{darkblue}r\) are the returns.

A strange regression problem: multiple-line fitting

Consider the following data:

----     14 PARAMETER x  data

i1  17.175,    i2  84.327,    i3  55.038,    i4  30.114,    i5  29.221,    i6  22.405,    i7  34.983,    i8  85.627
i9   6.711,    i10 50.021,    i11 99.812,    i12 57.873,    i13 99.113,    i14 76.225,    i15 13.069,    i16 63.972
i17 15.952,    i18 25.008,    i19 66.893,    i20 43.536,    i21 35.970,    i22 35.144,    i23 13.149,    i24 15.010
i25 58.911,    i26 83.089,    i27 23.082,    i28 66.573,    i29 77.586,    i30 30.366,    i31 11.049,    i32 50.238
i33 16.017,    i34 87.246,    i35 26.511,    i36 28.581,    i37 59.396,    i38 72.272,    i39 62.825,    i40 46.380
i41 41.331,    i42 11.770,    i43 31.421,    i44  4.655,    i45 33.855,    i46 18.210,    i47 64.573,    i48 56.075
i49 76.996,    i50 29.781


----     14 PARAMETER y  data

i1    32.320,    i2   238.299,    i3    87.081,    i4    17.825,    i5   -10.154,    i6   -35.813,    i7    14.506
i8   204.903,    i9  -103.014,    i10   79.521,    i11 -153.679,    i12  -47.971,    i13  277.483,    i14  160.084
i15   28.687,    i16  133.650,    i17  -63.820,    i18    4.028,    i19  145.434,    i20   48.326,    i21   49.295
i22   35.478,    i23  -58.602,    i24   34.463,    i25  -50.666,    i26 -115.522,    i27  -19.540,    i28  134.829
i29  198.074,    i30   -7.351,    i31   40.786,    i32   82.107,    i33   27.661,    i34  -82.669,    i35    5.081
i36   -1.278,    i37  118.345,    i38  172.905,    i39  -57.943,    i40   56.981,    i41  -45.248,    i42   34.517
i43   11.780,    i44  -98.506,    i45   -1.064,    i46   33.323,    i47  139.050,    i48  -35.595,    i49 -102.580
i50    3.912

When we plot this, we see:

QuadProg: usage issues

The QuadProg QP (Quadratic Programming) solver, based on an algorithm by Goldfarb and Idnani [1,2], is quite popular. It is available for R [3] and Python [4]. However, it has some real drawbacks. Some of it is just "user interface" and documentation: how to call the solver.

The mathematical model the solver works with is [3]: \[\begin{align} \min\>& -\color{darkblue}d^T\color{darkred}b + \frac{1}{2}\color{darkred}b^T\color{darkblue}D\color{darkred}b \\ & \color{darkblue}A^T\color{darkred}b \begin{bmatrix}= \\ \ge \end{bmatrix} \color{darkblue}b_0 \end{align} \]

This immediately poses some questions:

• Variables are typically called \(\color{darkred}x\), and \(\color{darkblue}b\) is usually a constant. Why the different notation here?
  I don't know. Maybe to make sure you pay attention.
• How to handle bounds?
  The variable type the solver supports is called "free" (i.e. unrestricted). This is different than for many LP solvers, where the default variable type is positive. For non-negative variables, you need to explicitly add the constraints \(\color{darkred}b \ge 0\). Similar for other bounds.
• Why the different signs for the linear part and the quadratic part in the objective?
  I don't know.
• Why the constant \(1/2\) in the objective?
  We see this more often in QP solvers. This is because the gradient of the quadratic term \(\color{darkred}x^T\color{darkblue}Q\color{darkred}x\) is \(2\color{darkblue}Q\color{darkred}x\). The algorithm prefers to work with \(\color{darkblue}Q\color{darkred}x\) however. So in order to stay really close to the algorithm, the user is asked to multiply all quadratic coefficients \(\color{darkblue}D_{i,j}\) by 2.
• Most QP solvers call the matrix with quadratic coefficients \(\color{darkblue}Q\). Why is it called \(\color{darkblue}D\) here?
  Don't know.
• Does \(\color{darkblue}D\) have to be symmetric?  
  The documentation does not tell, but the answer is: only the diagonal and upper-triangular part of \(\color{darkblue}D\) are used. This means the matrix is assumed to be symmetric. E.g. when passing \[\color{darkblue}D = \begin{pmatrix} 1 & 0 \\ 100 & 1 \end{pmatrix}\] the solver thinks it received \[\color{darkblue}D = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}\]
• Can \(\color{darkblue}D\) be positive semi-definite?
  In QuadProg, \(\color{darkblue}D\) has to be strictly positive definite.
  Defnition: A matrix \(Q\) is called positive semi-definite when \(x^TQx\ge 0\) for all \(x\). It is called positive definite when \(x^TQx\gt 0\) for all \(x \ne 0\). A positive definite matrix has strictly positive eigenvalues. I.e., the smallest eigenvalue of \(\color{darkblue}D\) must be \(\lambda_{min}\gt 0\). For a positive semi-definite matrix all eigenvalues are non-negative (i.e., zero is allowed).
  The fact that QuadProg requires a strict positive definite matrix, means that all variables need to be quadratic and also that this algorithm can not solve LPs. Most QP solvers (in fact all the ones I know beyond QuadProg) only require positive semi-definite quadratic coefficient matrices. I.e., for linear variables we can use coefficients \(\color{darkblue}q_{i,j}=0\) (where \(i\) or \(j\) corresponds to a linear variable). Quadprog is stricter. The whole matrix \(\color{darkblue}D\) most be strictly positive definite. You should be aware of these restrictions. And the documentation [3] should have mentioned this.
• Can I trick QuadProg to accept linear variables?
  Yes, suppose we have 2 quadratic variables and 2 linear variables. Then form the \(\color{darkblue}D\) matrix as \[\color{darkblue}D = \begin{pmatrix} \color{darkblue}d_{1,1} &\color{darkblue}d_{1,2} & 0 & 0  \\ \color{darkblue}d_{2,1} &\color{darkblue}d_{2,2} & 0 & 0  \\ 0 & 0 &\color{darkblue} \varepsilon & 0 \\ 0 & 0 & 0 &\color{darkblue} \varepsilon \end{pmatrix}\] where \(\color{darkblue} \varepsilon\gt 0\) is a small constant. This trick would also allow us to feed the algorithm linear problems.
• Usually, solvers express linear constraints as \(\color{darkblue}A\color{darkred}x=\color{darkblue}b\). Why do we have a transposed matrix \(\color{darkblue}A^T\) here?
  I suspect this is closer to what the dual algorithm prefers. Again, making the input closer to what the solver wants.
• How to input \(\le\) constraints?
  Multiply these by \(-1\) and make them \(\ge\) constraints.

These are some of the questions I have seen in the past about QuadProg.  Hopefully, this at least demystifies some of them.

As an alternative to QuadProg, you may want to have a look at CVXR or CVXPY for modeling QPs. These tools can talk to different solvers and will apply the necessary transformations to make the model suited for the selected solver. The modeling is also a bit more natural than the matrix input for QuadProg. For non-trivial quadratic programming models, QuadProg may not be the best choice.  

References

1. D. Goldfarb and A. Idnani (1982). Dual and Primal-Dual Methods for Solving Strictly Convex Quadratic Programs. In J. P. Hennart (ed.), Numerical Analysis, Springer-Verlag, Berlin, pages 226–239.
2. D. Goldfarb and A. Idnani (1983). A numerically stable dual method for solving strictly convex quadratic programs. Mathematical Programming, 27, 1–33.
3. R package quadprog, https://cran.r-project.org/web/packages/quadprog/quadprog.pdf
4. Python package quadprog, https://pypi.org/project/quadprog/