n-queens in Fusion

In a nice post on LinkedIn by Alireza Soroudi the author presents a short GAMS integer model for the classical n-queens puzzle: place the maximal number of non-attacking queens on an $n\times n$ chessboard.

We couldn't resist writing it up in Python Fusion as well.

And here is the solution we get for $n=8$:

[[_ ♕ _ _ _ _ _ _]

 [_ _ _ _ _ ♕ _ _]

 [♕ _ _ _ _ _ _ _]

 [_ _ _ _ _ _ ♕ _]

 [_ _ _ ♕ _ _ _ _]

 [_ _ _ _ _ _ _ ♕]

 [_ _ ♕ _ _ _ _ _]

  [_ _ _ _ ♕ _ _ _]]

Grouping in Fusion

We sometimes get asked how to efficiently perform a "group by" operation on a variable. That is, we have a variable $x=(x_0,\ldots,x_{n-1})$ naturally divided into groups and we want to add constraints for each group or constraints on aggregate values within groups. This arises for instance in portfolio optimization where the groups are assets, each consisting of a (varying) number of tax lots.

To keep the discussion more concrete, suppose we have a variable $x=(x_0,x_1,x_2,x_3,x_4,x_5)$ which consists of 3 groups $(x_0,x_1,x_2)$, $(x_3)$ and $(x_4,x_5)$. Let's say we want to express:

• an upper bound $b_i$ on the total value within each group,
• a joint volatility constraint such as $$\gamma\geq\left\|G\cdot \left[\begin{array}{c} x_0+x_1+x_2-i_0 \\ x_3-i_1 \\ x_4+x_5-i_2\end{array}\right]\right\|_2$$ for some matrix $G$ and constant index weights $i_0,i_1,i_2$
• the constraint that all values within one group have the same sign.

Pick, slice

A straightforward solution is to pick (Expr.pick) the content of each group into a separate view and add relevant constraints in a loop over the groups. One could also take slices (Expr.slice) when the groups form contiguous subsequences. This approach is implemented below in Python Fusion.

Loop-free

The previous solution requires picking and stacking expressions in a loop over all groups. This can sometimes be slow. A nicer and more efficient solution uses a bit of linear algebra to perform the grouping. We first encode the groups via a sparse membership matrix $\mathrm{Memb}$, where the rows are groups, columns are entries of $x$, and there is a $1$ whenever an entry belongs to a group. In our example $$\mathrm{Memb}=\left[\begin{array}{cccccc}1&1&1&0&0&0\\0&0&0&1&0&0\\0&0&0&0&1&1\end{array}\right] .$$

This matrix is easy to construct in Fusion. 

Note that $\mathrm{Memb}\cdot x$ is the vector of all group sums, in our case $$\mathrm{Memb}\cdot x = \left[\begin{array}{c} x_0+x_1+x_2 \\ x_3 \\ x_4+x_5\end{array}\right].$$

That means we can express both of the previous models in a single call without looping. Since $\mathrm{Memb}$ is very sparse, this becomes a very efficient representation. Of course it is important to keep $\mathrm{Memb}$ as a sparse matrix.

Loop-free same sign

We can now use the same matrix to model the last problem from the introduction: all entries within each group must have the same sign. We introduce a sequence of binary variables $z=(z_0,\ldots,z_{g-1})$, one for each of the $g$ groups. The $j$-th variable will determine the sign of elements in the $j$-th group. That is imposed by constraints $$-M(1-z_j)\leq x_i\leq Mz_j,$$ whenever $x_i$ belongs to $j$-th group. We can use the pick/slice strategy per group, as before, or observe that $\mathrm{Memb}^T\cdot z$ produces the vector with the correct binary variable for each entry in $x$. I our case, if $z=(z_0,z_1,z_2)$ then $$\mathrm{Memb}^T\cdot z = (z_0,z_0,z_0,z_1,z_2,z_2)^T.$$ Now each inequality in (4) can be written as a single constraint:

Harmonic mean, p-norms and extremely disciplined conic modeling

Contrary to what we inadverently claim in our modeling cookbook, there IS a conic quadratic model of the harmonic mean. More precisely, we can model the constraint
$$0\leq t\leq \frac{n}{\frac{1}{x_1}+\cdots+\frac{1}{x_n}},\quad x_i>0$$
Here is the model, which at the same time demonstrates the non-obvious fact that the right-hand side of (1) is actually concave. We call this approach extremely disciplined modelling, since by definition it prevents users from creating a non-convex model. We first notice that (1) is equivalent to
$$\sum_{i=1}^n\frac{t^2}{x_i}\leq nt$$
which can be written using rotated quadratic cones $2uv\geq w^2$, $u,v\geq 0$, as follows:
$$
\begin{eqnarray}
t^2\leq 2x_iy_i, &\quad i=1,\ldots,n \\
2\sum_i y_i = nt. &
\end{eqnarray}
$$
For the special case $n=2$ we can model $t\leq \frac{2xy}{x+y}$ with just one rotated cone: this inequality is equivalent to $2(x-\frac12 t)(y-\frac12 t)\geq (\frac12t)^2$.

The same way we can model (and prove concativity of) a general "negative $p$-norm":
$$ 0\leq t\leq (\sum_i x_i^{-p})^{-1/p} $$
where $p>0$, $x_i>0$ and the harmonic mean corresponds to $p=1$ (up to a factor $n$). The equivalent version is
$$ \sum_{i=1}^n\frac{t^{p+1}}{x_i^p}\leq t$$
which leads to
$$
\begin{eqnarray}
|t|\leq x_i^{\frac{p}{p+1}} y_i^{\frac{1}{p+1}}, &\quad i=1,\ldots,n \\
\sum_i y_i = t. &
\end{eqnarray}
$$
Constraint (7) defines the three-dimensional power cone with exponents $(\frac{p}{p+1},\frac{1}{p+1})$. Conic models with the power cone will be expressible in Mosek 9, see https://themosekblog.blogspot.dk/2018/01/version-9-roadmap_31.html.

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!

What if the solver stalls?

In this post we will discuss why the solver get sometimes stalls and how does this relates with the solution quality. What does we mean for stall?

The solver get stalled whenever it cannot make any sensible progress towards the optimal solution.

An iteration of a continuous optimizer in MOSEK consists of
\[ x^+ = x + \alpha d_x\]
where

• $x$:  Is the current guess for the optimal solution.
• $d_x$: Is a direction computed in some way by the optimizer.
• $\alpha$: Is nonnegative step size.
• $x^+$: Is the updated guess for the optimal solution.

The interpretation is that the optimizer computes a search direction $d_x$ and then takes a step in that direction. The search direction is computed such that when moving in that direction the current solution will be moved closer towards an optimal solution. Normally the step size $\alpha$ is chosen such that the new trial solution $x^+$ is as close as possible to the optimal solution with respect to some merit function. 

If the search direction is badly chosen then the optimal step size might be very close to $0$, in which case the algorithm is stalled. So why does this happens in practice? Typically the search direction $d_x$ is given as the solution of a system of linear equations, say

\[H d_x = h,\]

where it can be proven that if the linear system is solved accurately then good progress is obtained. Unfortunately in practice due to fact that computations are done in finite precision then the search direction may be of poor quality leading to a stall.

To summarize due to numerical issues MOSEK may not be able to compute a good search direction and if that leads to an extremely short step then MOSEK stops and says STALLED. Thus it does not make sense to continue because the subsequent iterations will also leads to small steps and no progress at all.

Observe that when a stall occurs then usually the current solution is very close to an optimal solution and many cases the reported solution is good enough for most practical purposes. Recall that usually the optimization problem solved is approximation of the true problem and even in the best case MOSEK will only report an approximated solution to approximated problem. 

Let's visualize the algorithm execution using a simplified scheme: the feasible region (primal and dual) is the yellowish one and the solution is the green dot on its corner. Algorithm execution generates a sequence of iterates, the green dots, that hopefully will end up in 10.

Interior-point algorithm execution sketch.

For each iterates, a normalized primal and dual feasibility violation is computed, as well as a measure of optimality violation. Let's call them $r_p,r_d$ and $r_g$ respectively. Ideally, $||(r_p,r_d,r_g)||$ should be driven to zero, but in practice the solver stops when
\[ ||(r_p,r_d,r_g)|| \leq \epsilon_f.\]

Thus $\epsilon_f$ defines a region in which a solution is optimal, i.e. the purple region in the figure. When the solver terminates it provides a solution with a certain solution status. Hopefully, it will be optimal. Otherwise, the solver may still be able to provide a useful solution: for this reason, and to improve flexibility, MOSEK denotes some solutions as near-optimal: the solution would be optimal if we scale up $\epsilon_f$ by a factor $\alpha>1$, i.e.

\[ ||(r_p,r_d,r_g)|| \leq \epsilon_f\alpha.\]

This corresponds to the reddish region. As the algorithm proceeds, we reach the near-optimal region at step (9) and then the optimal one at step (10). If all goes well, we stop in the latter. If we get stalled at (9) we will get a near-optimal solution. But what if the solver get stalled at step (8)? The solution status is unknown in this case.

In general, if the solver is stalled, the solution status is still meaningful since the solver is able to check the solution feasibility and the residual values. In particular a solution could be near-optimal. In our example, step (8) is not even near optimal and could be not very useful, but it is still worth to check it.
In case of stall, user should

• look at solution summary: provides the primal and dual violation under the column PFEAS and DFEAS, and the MU column measure the optimality violation;
• retrieve the solution status from the solver.

What can be done to avoid stalls? It is hard to say something about in general. A few good rules to follow are:

• Avoid nearly feasible or nearly infeasible problems.
• Make sure the problem is scaled reasonably.
• Make sure the primal and dual solutions are reasonably bounded i.e. not too large in absolute size.

Users should consider revising their model is stalling happens regularly. If modeling does not help, you are encouraged to contact our support.

MOSEK modeling manual

A new manual about how to model optimization problems is now available at docs.mosek.com. It provides many tips and tricks about how to formulate optimization problems so they can be solved by MOSEK.