Design a heart with DJC

In connection with the approaching Christmas, we are receiving many support questions about the use of mixed-integer conic optimization in the design of the traditional Danish "braided heart" Christmas tree decorations:

Let us demonstrate how this can be done with MOSEK. First, we need to define the heart. A suitable model is given by

$$x^2+(y-p|x|)^2\leq r^2,$$

where $p,r$ are parameters which can be adjusted for different shapes and sizes:

This is not a convex shape, but it can still be modelled in conic form if we employ a mixed-integer model of $|x|$. This can be done using a classical big-M approach or, in MOSEK 10, using a disjunctive constraint (DJC):

$$\begin{array}{l}(x\geq 0 \ \mathrm{AND}\ z=x)\ \mathrm{OR}\ (x\leq 0 \ \mathrm{AND}\ z=-x) \\  x^2+(y-pz)^2\leq r^2\end{array},$$

combining a DJC with a quadratic cone constraint in affine conic (ACC) form.

Having defined the domain, it is time to draw some lines with slope $\pm 45^\circ$ to determine the braided pattern. To do this, we will choose, one at a time, line segments with integer endpoints, say $x=(x_1,x_2)$, $y=(y_1,y_2)$, inside the heart. Such a line segment has the correct slope if

$$x_1-x_2=y_1-y_2\ \mathrm{OR}\ x_1+x_2=y_1+y_2,$$

which, again, is a disjunctive constraint. We will only choose sufficiently long line segments, for example by requiring

$$|x_1-y_1|\geq c$$

for some constant $c$. Here the absolute value can again be modelled in the same mixed-integer fashion as before. We set no objective, so this becomes a feasibility problem, but one could also decide, for instance, to maximize the length of the segment or optimize some other measure of beauty for the design.

Having found one segment, we will construct a new one by adding additional constraints to guarantee that the new segments do not interfere too much with the previous ones. For the sake of simplicity, we will demand that new endpoints do not belong to any line of slope $\pm 45^\circ$ containing the previous endpoints, but one could consider other elimination conditions as well. This iterative procedure is similar to solving the travelling salesman problem by cycle elimination, and it terminates once there is no more place for new segments under all the constraints. The implementation requires us to be able to write the "not equals" condition $x\neq a$, which, assuming we work with integers, can again be cast as a DJC:

$$x\leq a-1 \ \mathrm{OR}\ x\geq a+1.$$

The iterative procedure will construct more and more segments:

To obtain different designs, we can vary the parameters $p,r$. Still, there is enough freedom in the choice of segments that even just varying the random seed ("mioSeed") for the mixed-integer optimizer produces different final solutions for the same choice of $p,r$:

And on this note, we leave you with the Fusion API code below and wish you all the best. (You can also read our last Christmas special from way back in 2018).

MOSEK Team

	from mosek.fusion import *
	import numpy as np
	import matplotlib.pyplot as plt
	import sys, itertools
	# For a scalar variable x returns a new variable
	# exactly representing abs(x)
	def abs(M, x):
	absx = M.variable()
	# (x>=0 AND absx==x) OR (x<=0 AND absx==-x)
	M.disjunction(DJC.AND(DJC.term(x, Domain.greaterThan(0)),
	DJC.term(Expr.sub(absx, x), Domain.equalsTo(0))),
	DJC.AND(DJC.term(x, Domain.lessThan(0)),
	DJC.term(Expr.add(absx, x), Domain.equalsTo(0))))
	return absx
	# Point belongs to a heart with parameters (p,r)
	def inHeart(M, x, p, r):
	x1, x2 = x.index(0), x.index(1)
	# r^2 >= x1^2 + (x2-p*absx1)^2
	M.constraint(Expr.vstack(Expr.constTerm(r),
	x1,
	Expr.sub(x2, Expr.mul(p, abs(M, x1)))),
	Domain.inQCone())
	# An integer variable or expression not equals a
	def different(M, x, a):
	# x>=a+1 OR x<=a-1
	M.disjunction(DJC.term(x, Domain.greaterThan(a+1)),
	DJC.term(x, Domain.lessThan(a-1)))
	##################################
	# Model construction starts here #
	##################################
	# Set some parameters
	p = 1.1
	r = 8.1
	seed = 77
	mindist = 7
	# Initialize the model
	M = Model("Christmas Heart")
	M.setSolverParam("mioSeed", seed)
	M.setLogHandler(sys.stdout)
	lines = []
	# Coordinates of the two points
	x = M.variable(2, Domain.integral(Domain.unbounded()))
	y = M.variable(2, Domain.integral(Domain.unbounded()))
	inHeart(M, x, p, r)
	inHeart(M, y, p, r)
	# Make each line sufficiently long
	dist = abs(M, Expr.sub(x.index(0), y.index(0)))
	M.constraint(dist, Domain.greaterThan(mindist))
	# Each x - y line has slope 45deg in either direction
	M.disjunction(DJC.term(Expr.dot(Var.vstack(x, y),[1,-1,-1,1]), Domain.equalsTo(0)),
	DJC.term(Expr.dot(Var.vstack(x, y),[1,1,-1,-1]), Domain.equalsTo(0)))
	# Add new lines as long as possible
	M.solve()
	while M.getProblemStatus() != ProblemStatus.PrimalInfeasible:
	lines.append((x.level(),y.level()))
	# Forbid endpoints on any 45deg sloped line passing through an existing point
	for (v1, v2, v3) in itertools.product([x,y],[[1,-1],[1,1]],[x.level(),y.level()]):
	different(M, Expr.dot(v1, v2), np.dot(v2, v3))
	# Resolve with the new constraints
	M.solve()
	# Final plot
	u = np.linspace(-1.5*r, 1.5*r, 101)
	v = np.linspace(-1.5*r, 1.5*r, 101)
	z = np.zeros(shape=(len(u), len(v)))
	for (i,j) in itertools.product(range(len(u)),range(len(v))):
	z[i,j] = u[i]**2 + (v[j]-p*np.abs(u[i]))**2 - r**2
	plt.contour(u, v, z.T, [0])
	for ((a,b),(c,d)) in lines:
	plt.scatter([a], [b])
	plt.scatter([c], [d])
	plt.plot((a,c),(b,d))
	plt.grid()
	plt.show()

view raw heart.py hosted with ❤ by GitHub

Christmas 2022

During the Christmas period our sales and support will be closed 23-27 December 2022, inclusive. 

The MOSEK Team

MIP 2023 - 20th Anniversary

MOSEK is proud to be one of the sponsors of MIP 2023 - Mixed Integer Programming workshop. MIP 2023 marks the 20th year of the conference and, as such will be particularly special. It will be held May 22-25, 2023, at the University of Southern California, Los Angeles.

More details at https://www.mixedinteger.org/2023/

A customer reported one order of magnitude of speed up for his problems.

A happy customer just reported back that Mosek can solve some of their mixed-integer problems one order of magnitude faster using the recently released version 10!

The origin of that speed improvement was that the customer a couple of years ago provided a set of typical problems that they solved i.e. problems they would like Mosek to perform better on. Now during the development of Mosek version 10 the mixed-integer team repeatedly looked at how the proposed changes to the mixed-integer optimizer worked for the provided benchmark problems. This has lead to numerous small improvements summing up to a large impromvement for the problems of this customer while not degrading the performance on average.

The take away message form this story is that if you provide a set of benchmark problems, then you are much more likely to observe an improvement in performance of Mosek over time for your problems.It may take some time before happens, but it is likely to happen.

A good benchmark set contain problems of diffrent difficulty and the most of the problems should represent the normal difficulty. Fell free to contact Mosek support to provide a set of benchmark problems.

2022 INFORMS Annual Meeting

MOSEK is proud to be one of the sponsors of the INFORMS Optimization Society at the 2022 INFORMS Annual Meeting, taking place on  October 16th–19th, 2022 at the JW Marriott Hotel Convention Center, Indianapolis, USA.

More details at https://meetings.informs.org/wordpress/indianapolis2022/

MOSEK 10 is out

You can now install and use MOSEK 10, which is the current stable release.

More information and links can be found in our September 2022 newsletter.

MOSEK 10 release approaching

MOSEK version 10 with all its new features has been available as a beta version for a while now. See https://www.mosek.com/products/version-10/ for the latest details and download and installation instructions.

Huge thanks to all the users who were so kind to test this beta release and share their suggestions, bug reports and comments. We are working on incorporating your input into the final release. Do not hesitate to write to us with further remarks.

MOSEK 10 is planned to become the current stable release latest on August 15th 2022.

In the meantime, we wish all MOSEK users a good summer, enjoyable vacation, and hope to meet you with version 10 after the summer.

The MOSEK team

On the performance improvement of the semi-definite optimizer when switching from version 9 to 10

The most reason Mosek version 10.0.16(beta) includes a very significant improvement of the optimizer for semi-definite problems. The improvement is most notable for problems having many matrix variables. In order to illustrate the potential performance gain we use a problem supplied by Mosek user. Some statistics about the problem can be seen in the Mosek log output:

Optimizer  - solved problem         : the primal
Optimizer  - Constraints            : 8808
Optimizer  - Cones                  : 0
Optimizer  - Scalar variables       : 0                 conic                  : 0
Optimizer  - Semi-definite variables: 49                scalarized             : 76493
Factor     - setup time             : 8.26              dense det. time        : 0.00
Factor     - ML order time          : 0.50              GP order time          : 0.06
Factor     - nonzeros before factor : 3.38e+07          after factor           : 3.38e+07
Factor     - dense dim.             : 0                 flops                  : 2.45e+11
Factor     - GP saved nzs           : 4.97e+06          GP saved flops         : 4.60e+10

Hence, the problem has 8808 constraints and 49 matrix variables. In addition it should be mentioned that the contraint matrix is not that sparse.

Below the computational results is shown.

CPU	9.3.20	9.3.20 (MT)	10.0.16	10.0.16 (MT)
Intel Xeon Gold 6126 (12 core)	1165	1171	833	207
AMD EPYC 7402P (24 core)	1189	1049	777	108

The table shows the computational time in wall clock seconds for running version 9.3.20 and 10.0.16 single and multi-threaded respectively. Moreover, the timing result is shown for an Intel and an AMD CPU.

The main conclusions from the table are

• version 10 is much faster than version 9 and the difference is more notable on the AMD CPU.
• version 10 offers much better speed-up when using multiple threads. In fact there are virtually no speed-up when using multiple threads in version 9.

It should emphasized that the performance improvement for solving semi-definite problems when switching from version 9 to version 10 is dependent on the structure of the problem. However, as illustrated above a big performance boost can be realized for some problems.

PS. Despite the AMD CPU being the fastest in the table above then this post do not claim AMD CPUs are better than Intel CPUs or the other way around.

Symmetry detection for MIPs in MOSEK 10

Remember counting ladybug spots as a child? If you realized ladybugs were symmetric, you could save some work by counting on one wing only and multiplying by two. In Mixed Integer Programming, something similar applies. If we detect symmetries in a problem formulation, we may be able to save some computational effort and thus time when solving it. But first of all, what are symmetries in a MIP formulation? Different from the geometric symmetry of a ladybug, the notion of symmetries in MIP may be thought of as the presence of equivalent solutions.

They can arise, for example, when identical physical objects are described distinctively in a model. As an example, take the MISOCP models for geometric facility location problems described here. Roughly speaking, the task there is to cover a set of given points in the plane with disks to be placed at will. Say we have 10 points and 3 fixed-radius disks, and want to cover as many points as possible. In a straightforward model, we would assign indexes $i = 1,2,3$ to the disks, and indexes $j = 1,...,10$ to the points, and then define all sorts of variables like $x_{ij} \in \{0,1\}$, stating whether disk $i$ covers point $j$. What we did here is assign labels to the disks, making them a blue, a red and a green one, say. But if all three disks have the same radius, this gives rise to equivalent solutions, see Figure 1. For the cost of a solution, it is only important whether a certain point is covered, but not by which disk.

Figure 1: Equivalent disk covering solutions

Symmetry breaking inequalities, orbital branching and fixing, or Isomorphism pruning are just some of the various techniques that have been studied over the years to better address symmetries in MIP. Variants of some of these various techniques have been implemented in the latest MOSEK version 10. Figure 2 sketches the effect on solution times of completely switching off symmetry detection and handling with the user parameter MSK_IPAR_MIO_SYMMETRY_LEVEL on the so-called Margot instances, known to be highly symmetric.

Figure 2: Solving Margot's instances w/ and w/o symmetry detection

Some of the techniques in question can be extended from Mixed Integer Linear to Mixed Integer Conic Programming, or more generally Mixed Integer Nonlinear Programming. We sometimes advocate that one nice feature of (Mixed Integer) Conic Programming is the separation of data and structure, and implementing symmetry detection is an aspect where this comes in handy. All data involved in a model is contained in linear constraints like $a^Tx\leq b$, and the symmetry detection routines needed here are already needed in MILP. Now extending symmetry detection to conic constraints, e.g., to a conic-quadratic constraint $x_0\geq \Vert x_{1:n}\Vert_2$, is relatively easy because one does not have to hassle with data, but only structure. Adding symmetry detection for possibly newly added cone types to a MICP solver is thus a small burden. And detecting and exploting symmetry in MICP can have a decisive effect. The following logs show MOSEK's solution paths to solve one of the aforementioned geometric facility location problems, with and without symmetry detection and handling, respectively.

Figure 3: Solving disk covering instances w/ and w/o symmetry detection

Easter 2022

Due to Easter support and sales are closed from Thursday, April 14th, until Monday, April 18th, inclusive.

The MOSEK team wishes everyone a peaceful Easter.

Improved performance of the interior-point optimizer when using an AMD CPU

One of the improvements in the recent released Mosek version 10 is better performance when run on a computer that has a recent AMD CPU. In order to illustrate that the table below is shown.

Problem	9.3	10.0
pds-100	93s	66s

The table shows the time taken in seconds to solve some well-known linear optimization problems on a AMD EPYC 7543 CPU from the Zen3 familie running the interior-point optimizer on 1 thread. On this AMD Zen3 family CPU Mosek version 10 is approximately 30% faster which is a significant improvement. It should be emphasized that on easier problems than pds-100 then the speed improvement is less than reported above. Nevertheless an upgrade to Mosek version 10 should be particularly worthwhile when using an AMD CPU.

MOSEK 10 Beta release

We are excited to announce the release of MOSEK 10 Beta. 

This version introduces many new features, in particular:

• Support for Apple M1.
• Improved control over multithreading.
• Efficiency improvements in the interior-point and mixed-integer optimizers, improved tuning for various platforms.
• New strategies in the mixed-integer optimizer.
• Disjunctive constraints.
• More convenient syntax for conic constraints in the Optimizer API.
• And many more improvements.

To read more and try the new version see https://www.mosek.com/products/version-10/

In the coming weeks we will publish a series of posts highlighting selected new features of MOSEK 10.

Your feedback will be most appreciated.

The MOSEK team

MIP 2022

Similarly to previous years MOSEK is proud to be one of the sponsors of MIP2022 - Mixed Integer Programming Workshop 2022, taking place May 23–26, 2022 at DIMACS, Rutgers University.

More details at https://www.mixedinteger.org/2022/

Planned end of support for version 8

On 31-may-2022 we will be stopping active support for version 8 which was released back in 2016.

For details of release and support policy see https://www.mosek.com/products/release-policy/