#!/usr/bin/env python
# coding: utf-8
# ![MOSEK ApS](https://www.mosek.com/static/images/branding/webgraphmoseklogocolor.png )
# Smallest sphere enclosing a set of points.
# ===========================
#
# The aim of this tutorial is two-fold
#
# 1. Demostrate how to write a conic quadratic model in `Fusion` in a very simple and compact way.
# 2. Show how and way the dual formulation may solved more efficiently.
#
#
# Our problem is the defined as:
#
# **Find the smallest sphere that encloses a set of** $k$ **points** $p_i \in \mathbb{R}^n$.
# In[1]:
get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import random
def plot_points(p, p0=[], r0=0.):
n,k= len(p0), len(p)
plt.rc('savefig',dpi=120)
fig, ax = plt.subplots()
ax.set_aspect('equal')
ax.plot([ p[i][0] for i in range(k)], [ p[i][1] for i in range(k)], 'b*')
if len(p0)>0:
ax.plot( p0[0],p0[1], 'r.')
ax.add_patch( mpatches.Circle( p0, r0 , fc="w", ec="r", lw=1.5) )
plt.grid()
plt.show()
n = 2
k = 500
p= [ [random.gauss(0.,10.) for nn in range(n)] for kk in range(k)]
plot_points(p)
# The problem boils down to determine the sphere center $p_0\in \mathbb{R}^n$ and its radius $r_0\geq 0$, i.e.
#
#
# \begin{equation}
# \begin{aligned}
# \min \max_{i=1,\dots,k} \| p_0 - p_i\|_2 \\
# \end{aligned}
# \end{equation}
#
# The maximum in the objective function can be easily, i.e.
#
# \begin{equation}
# \begin{aligned}
# \min r_0 & & &\\
# s.t.& & &\\
# & r_0 \geq \| p_0 - p_i\|_2 ,& \quad & i=1,\ldots,k\\
# \end{aligned}
# \end{equation}
#
# The SOCP formulation reads
#
# \begin{equation}
# \begin{aligned}
# \min r_0 & & &\\
# s.t.& & &\\
# & \left[r_0,p_0 - p_i\right] \in Q^{(n+1)}, & \quad & i=1,\ldots,k.
# \end{aligned}
# \end{equation}
# Before defining the constraints, we note that we can write
#
#
# \begin{equation}
# R_0 = \left( \begin{array}{c} r_0 \\ \vdots \\ r_0 \end{array} \right) \in \mathbb{R}^k , \quad
# P_0 = \left( \begin{array}{c} p_0^T \\ \vdots \\ p_0^T \end{array} \right) \in \mathbb{R}^{k\times n}, \quad
# P = \left( \begin{array}{c} p_1^T \\ \vdots \\ p_k^T \end{array} \right) \in \mathbb{R}^{k\times n}.
# \end{equation}
#
# so that
#
# \begin{equation}
# \left[r_0,p_i - p_0\right] \in Q^{(n+1)}, \quad i=1,\ldots,k.
# \end{equation}
#
# can be compactly expressed as
#
# \begin{equation}
# \left[ R_0,P_0-P\right] \in \Pi Q^{(n+1)},
# \end{equation}
#
# that means, with a little abuse of notation, that each rows belongs to a quadratic cone of dimension $n+1$.
#
#
# Now we are ready for a compact implementation in `Fusion`:
# In[2]:
from mosek.fusion import *
import mosek as msk
def primal_problem(P):
print(msk.Env.getversion())
k= len(P)
if k==0: return -1,[]
n= len(P[0])
with Model("minimal sphere enclosing a set of points - primal") as M:
r0 = M.variable(1 , Domain.greaterThan(0.))
p0 = M.variable([1,n], Domain.unbounded())
R0 = Var.repeat(r0,k)
P0 = Var.repeat(p0,k)
M.constraint( Expr.hstack( R0, Expr.sub(P0 , P) ), Domain.inQCone())
M.objective(ObjectiveSense.Minimize, r0)
M.setLogHandler(open('logp','wt'))
M.solve()
return r0.level()[0], p0.level()
# We will also store the solver output in a file to use it later on. And then just solve the problem.
# In[3]:
r0,p0 = primal_problem(p)
print ("r0^* = ", r0)
print ("p0^* = ", p0)
# In[4]:
plot_points(p,p0,r0)
# Dual Formulation
# -----------------------
#
# The dual problem can be determined in few steps following basic rules. Introducing dual variables
#
# \begin{equation}
# Y = \left( \begin{array}{c} y_1^T\\ \vdots \\y_k \end{array}\right), \quad z = \left( \begin{array}{c} z_1\\ \vdots \\z_k \end{array}\right),
# \end{equation}
#
# the dual is:
#
# \begin{aligned}
# \max & \left\langle P, Y \right\rangle \\
# & e_k^T z = 1\\
# & Y^T e_k = \mathbf{0}_n \\
# & \left[z_i , y_i\right] \in \mathcal{Q}^{n+1}\\
# & z_i\in \mathbb{R}, y_i\in \mathbb{R}^n,
# \end{aligned}
#
# where $e_k\in \mathbb{R}^k$ is a vector of all ones.
#
# The ``Fusion`` code is the following:
#
# In[5]:
def dual_problem(P):
k= len(P)
if k==0: return -1,[]
n= len(P[0])
with Model("minimal sphere enclosing a set of points - dual") as M:
Y= M.variable([k,n], Domain.unbounded())
z= M.variable(k , Domain.unbounded())
M.constraint(Expr.sum(z), Domain.equalsTo(1.) )
e= [1.0 for i in range(k)]
M.constraint(Expr.mul(Y.transpose(), Matrix.ones(k,1)), Domain.equalsTo(0.) )
M.constraint( Var.hstack(z,Y), Domain.inQCone())
M.objective( ObjectiveSense.Maximize, Expr.dot( P, Y ))
M.setLogHandler(open('logd','wt'))
M.solve()
return
dual_problem(p)
# Looking at the solver output for the primal and dual implementation we get
# In[6]:
get_ipython().system('tail logp')
# In[7]:
get_ipython().system('tail logd')
# **The solver has got the same result doing exacly the same steps!**
#
# In fact, **MOSEK** has solve the same problem, thanks to the automatic dualizer that decides whether is more convenient to solve the primal or the dual of a given problem. To see that, just look at the output again:
# In[8]:
get_ipython().system('grep Optimizer logp')
# In[9]:
get_ipython().system('grep Optimizer logd')
# In primal problem output the solver reports it is going to solve the `dual`, i.e. our second formulation. In the dual problem output it says it is going to solve the `primal`. Therefore in both cases we actually solve the very same formulation.
#
This work is licensed under a Creative Commons Attribution 4.0 International License. The **MOSEK** logo and name are trademarks of Mosek ApS. The code is provided as-is. Compatibility with future release of **MOSEK** or the `Fusion API` are not guaranteed. For more information contact our [support](mailto:support@mosek.com).
#
# In[ ]: