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