Modeling non-convex absolute value constraints with MILP

We are being asked if one can use MOSEK to express a constraint of the form $$\label{eq:sumabs}\sum_i|x_i|=c$$ for some $x\in\mathbb{R}^n$ and $c\in\mathbb{R}$. This condition is non-convex; for instance for $n=1$ it is equivalent to $x=\pm c$. Such constraints appear in long-short portfolio optimization problems (see stackoverflow, paper).

The idea is to introduce a binary vector indicating the positions of positive/negative entries in $x$. Concretely, we want to split $x$ into the positive and negative part $x=s-t$ with $s,t\geq 0$.

If we know an upper bound $M$ on $s_i$ and $t_i$, then the following system:
$$\begin{align*}
0\leq s_i&\leq My_i,\\
0\leq t_i&\leq M(1-y_i),\\
& y_i\in \{0,1\}
\end{align*}$$ has the property that at most one of $s_i,t_i$ is non-zero. Indeed:
$$
\begin{align*}
s_i>0 \implies y_i=1 \implies t_i=0,\\
t_i>0 \implies y_i=0 \implies s_i=0.
\end{align*}
$$ Given that we are looking at $\sum_i|x_i|=c$ we can choose $M=c$ as an upper bound for all $|x_i|$. Then an equivalent version of the condition $\eqref{eq:sumabs}$ as a mixed-integer linear program is:
$$
\begin{align*}
x&=s-t,\\
0\leq s&\leq cy,\\
0\leq t&\leq c(1-y),\\
0\leq y&\leq 1,\\
 c&=\sum_is_i+\sum_it_i,\\
& x,s,t\in\mathbb{R}^n,\ y\in\mathbb{Z}^n.\\
\end{align*}
$$

New phone number

We have gotten a new phone number. It is

• +45 7174 9373

Also the fax number has been closed. 

Improving Python API performance

Recently a MOSEK user pointed out out that sometimes the MOSEK Python API can be quite slow when it comes to input large amount of data. In particular the problem is related to the way MOSEK handles 32 and 64 bit integers.


If you use Numpy with a 64-bit Python, integer arrays will usually end up being arrays of numpy.int64 since the native Python integer is an int64.


Most MOSEK functions require arrays of numpy.int32, e.g. for variable or constraint indexing. In this case it can often significantly speed up the function call to explicitly convert the types.

To this extent is worth using numpy.array class to convert arrays of indexes to the right size: the numpy implementation is highly efficient and saves a lot of time!

Tests

To show the performance gain in Python 2.7 we implement a simple tests involving

• n=1000000 variables,
• linear objective function terms and
• quadratic objective function terms.

using Python timeit module. Each tests is repeated num=100 times and we removed from the timing unecessart contributions such as module import and task/environment declaration.

Note we show source code without timeit instrumentation to keep it simple. 

The complete testing code is available at the end of the post.

Objective function linear coefficents.

Say we want to input the coefficients for n variables in the objective function. We can do it using Python arrays or  numpy arrays instead. The code to time the running time is

We obtain the following results

• direct input: 1.113329s
• using numpy : 0.017299s

Objective function quadratic coefficients.

Say we want to input the coefficients of the linear part of the objective function. The benchmark code is

It provides the following results

• direct input: 2.239344s
• using numpy : 0.028208s

Conclusion

The MOSEK Python API converts by default 64bit integers to 32 bit ones, and this operation could introduce significant overhead. Using numpy array to convert integer arrays to 32 bit is a simple simple workaround that allows to speed up the input process significantly. 

The improvement is quite remarkable. However, it must be noticed that not all users will experience such an improvement. Nevertheless it is a good idea to check your code and see whether you can gain some performance.

Remember that this issue is solved in the new MOSEK 8 release.

Complete code

MOSEK 8 beta now available!

We are happy to announce that MOSEK version 8 beta is now available for download from our web site:

https://mosek.com/products/beta-versions

Main highlights:

• The conic optimizer ha been improved in many respect: pre-solving, accuracy, speed.
• The semidefinite optimizer in particular has been greatly improved in terms of accuracy and stability.
• The Fusion API is now available for Python 3 and C++.
• An automatic dualizer for conic quadratic problem is now available.
• A simple optimization server to solve problems remotely is now included.
• Linux 32bit is no more supported.

We have also updated may other part of the solver and polished its interface. The documentation has been also revised.

Want to try this new release? Download the package and get an updated license:

• Trial and academic licenses can be used with version 8 beta as well. 
• Commercial users that are entitled for upgrade can contact us to obtain a new license file.

We are very much confident this new release will be a further improvement for our product. 

Help us improving even more: let us know your feedback! 

Do you want to know more about the new release? Contact us!

The MOSEK Team.

Sponsorship: 2016 MIP Workshop in Coral Gables, Florida (US)

We are happy to announce that MOSEK is among the sponsors of the 2016 Mixed Integer Programming Workshop that will be held May 23rd - 26th in Coral Gables, Florida (US).

The program includes talks from many distinguished specialists in the field of mixed-integer programming both from the academic and industrial world. No doubt the level of the presentations will be outstanding.

This year MOSEK not only support the workshop as sponsor, but we are part of the Program Committee.

We look forward for such an exciting event!

The MOSEK team

Sponsorship: AOO 2016 - Copenhagen (DK) - May the 3rd 2016

MOSEK is once more one of the sponsors of the Applications of Optimization 2016 day (AOO 2016), the annual meeting of the Danish Operational Research Society (DORS).

The AOO day is meant to be an occasion for all operation research practitioners to meet, exchange experiences and find new inspiration. The program includes six application oriented talks delivered by experienced OR specialists. Maritime logistic and oil industry are among the topics.

The quality of the talks and speakers, the informal atmosphere and the clear focus on practical applications make AOO a nice and fruitful  day.

The meeting will be held at Industriens Hus, in the very heart of Copenhagen (DK) on May 3rd, 2016. Please check the official web site for more information.

We look forward to meet all the AOO 2016 participants and speakers!

The MOSEK team.

Reformulating a non convex MINLP as MISOCP

Few weeks ago we read a nice and interesting post here.

The author took inspiration from a previous post on an open forum to show how the intuitive formulation of a simple problem can lead to a non convex MINLP when a MIQP can be used instead.

In this post we will

• show how to reformulate the proposed MIQP problem as an MISOCP,
• implement the MISOCP model using MOSEK Fusion API,
• derive a more compact MISOCP formulation and
• implement the second model using MOSEK Fusion API as well.

The problem is the following: we are given $N$ power unit (in the original post ovens) that must be used to provide a required amount $T$ of energy. Each plant contributes an amount $x_i$ of energy.

The inefficiency of each plant is given by

\[
g_i(x_i) = \left\lbrace \begin{array}{ll} 100 \sum_ i \frac{(a_i - x_i)^2}{a_i} & x_i>0,\\
0& x_i=0.  \end{array}\right.
\]

where $a_i$ is the optimal operational level of the plant. Moving away from $a_i$ leads to inefficient use of the plant. Here there is an example for $a_i = 10 $.

The overall inefficiency is then

\[
f(x) = \sum_i g_i(x_i).
\]

We want to minimize the inefficiency still providing the amount $T$ of energy.

A simple approach is to use an indicator variable $y_i$ to represent the $i-$th plant activation. Assuming we know an upper bound $M_i$ to the power each plant can produce, the model takes the form

\begin{equation}
\begin{array}{lll}
\min & 100\sum_i y_i \frac{(a_i - x_i)^2}{a_i} \\
&\\
s.t. & \sum_i x_i = T\\
& x_i \leq y_i M_i & i=1, \ldots, N\\
& x_i \geq 0 & i=1, \ldots, N,
& y_i \in \{0,1 \} & i=1, \ldots, N.\\ \end{array}
\end{equation}


As noted in the original post, this is a non-convex NLP. A simple reformulation allows to remove the bilinear term $x_i y_i$:



\begin{equation}
\begin{array}{lll}
\min & \sum_i \frac{d_i^2}{a_i} \\
&\\
s.t. & \sum_i x_i = T\\
& x_i \leq y_i M_i & i=1, \ldots, N\\
& -M_i(1- y_i) \leq s_i \leq M_i(1-y_i)  & i=1, \ldots, N\\
& d_i = a_i - x_i +s_i  & i=1, \ldots, N\\
& x_i \geq 0 & i=1, \ldots, N,\\
& y_i \in \{0,1 \} & i=1, \ldots, N.\\ \end{array}
\end{equation}

The trick is simple: an active plant $y_i=1$ implies $s_i=0$ and therefore $d_i = a_i - x_i$. But for $y_i =0$ the slack variable $s_i$ is free, while $x_i=0$. This implies that $d_i$ is free to assume any value, and since we are minimizing its second power, it will be driven to zero as well.

In conic form the problem reads



\begin{equation}
\begin{array}{lll}
\min & t\\
&\\
s.t. & \sum_i x_i = T\\
& x_i \leq y_i M_i & i=1, \ldots, N\\
& -M_i(1- y_i) \leq s_i \leq M_i(1-y_i)  & i=1, \ldots, N\\
& d_i = a_i - x_i +s_i  & i=1, \ldots, N\\
& (1/2, t, d_1/\sqrt{a_1}, \ldots,  d_n/\sqrt{a_N}) \in \mathcal{Q_r}\\
& x_i \geq 0 & i=1, \ldots, N,\\
& y_i \in \{0,1 \} & i=1, \ldots, N.\\ \end{array}
\end{equation}

Those not used to conic formulation may refer to our modeling cookbook.

The Fusion implementation is the following 

The code is very compact and readable, no need for many comments. The only operators one may wonder about are

• vstack - it returns an array stacking its argument vertically,  
• mulElm - it performs element wise multiplications.

Now we ask: is there a more compact formulation? 

The answer is yes: the purpose of the slack variables is to remove the contribution to the inefficiency that a deactivated oven will provide anyway. That contribution is simply equal to $a_i$. Therefore we only need to subtract $a_i$ if $y_i=1$.

\begin{equation}
\begin{array}{lll}
\min & t - \sum_i a_i ( 1- y_i)\\
&\\
s.t. & \sum_i x_i = T\\
& x_i \leq y_i M_i & i=1, \ldots, N\\
& (1/2, t, (a_1 -x_1)/\sqrt{a_1}, \ldots, (a_N - x_N)/\sqrt{a_N}) \in \mathcal{Q_r}\\
& x_i \geq 0 & i=1, \ldots, N,\\
& y_i \in \{0,1 \} & i=1, \ldots, N.\\ \end{array}
\end{equation}



which translates to the following code:





As one can see, the code is more compact, clean and readable. Fusion helps implementing conic optimization models in short time with a remarkable simplicity.


Last observation: the proposed problem is indeed simple to formulate, but the intuition leads to a non-convex MINLP. That is typical, but luckily in this case with a bit of work we can actually move to a MIQP formulation first, and a MISOCP then. And the final MISOCP comes with no additional variables at all.

Investing time on the model pays!

Personal Academic Licenses

We have recently experienced technical issues that have prevented us to provide personal academic licenses (see our recent blog post). Academic users trying to apply or renew personal academic licenses have got an error message instead.

The service is now up and running and academic users can obtain fully-featured personal academic licenses at

https://mosek.com/resources/academic-license

Please let us know if any other issues arises.

We do apologize for the inconvenience.


The MOSEK Team.

Trial license server not working

Our academic trial license server is currently malfunctioning. 

Due to a technical issue, academic license requester will get an error from our server.

The error has been introduced by a recent update.

We are working on the issue and hope to resolve it as soon as possible. In the mean time, please grab 30 days commercial trial licenses instead at

https://mosek.com/resources/trial-license

Contact us at info@mosek.com if any other issue arises.

We apologize for the inconvenience.

The MOSEK team