Optimal parade route

Topic references

Iterated weighted :math:\ell_1 heuristic: slides <https://stanford.edu/class/ee364b/lectures/l1_ext_slides.pdf>__

Problem statement

In this notebook, we'll tackle the problem of trying to protect a parade route by placing a limited number of guards.

The parade route is discritized into :math:m points. There are :math:n possible guard locations with associated decision variable :math:x \in \lbrace 0,1\rbrace^n, where :math:x_i = 1 if and only if a guard is placed at location :math:i. Associated with guard location :math:i is a coverage vector :math:a_i \in \mathbf{R}^m, which describes how well a guard placed at location :math:i would 'cover' each point in the parade route. We will assume that guard coverage is additive, so that the vector describing the total coverage of every edge is given by :math:Ax, where :math:A \in \mathbf{R}^{m \times n} has :math:a_i as its :math:i\ th column.

The parade route is only as secure as its least well-covered point. Our goal is to place :math:k guards to maximize the minimum coverage over the points in the route.

Optimization formulation

We can formulate this as the optimization problem

.. math::

\begin{array}{ll} \mbox{maximize} & t \ \mbox{subject to} & t \leq Ax\ & 0 \leq x \leq 1 \ & \mathbf{1}^Tx = k, \end{array}

This problem is nonconvex and, in general, NP-hard due to the Boolean decision variable.

Relaxation

We can try to approach the problem with convex optimization by first forming the convex relaxation

.. math::

\begin{array}{ll} \mbox{maximize} & t - w^Tx\ \mbox{subject to} & t \leq Ax\ & 0 \leq x \leq 1 \ & \mathbf{1}^Tx = k, \end{array}

by constraining :math:x \in [0,1]^n.

In general, the solution to this problem, :math:x^\star, will have fractional values. As we want a Boolean allocation, we can use an iterated weighted :math:\ell_1 heuristic to try to recover a Boolean solution.

Iterated weighted :math:`\ell_1` heuristic

To try and recover a Boolean solution, we will solve a sequence of convex problems where we add a linear term :math:-w^Tx to the objective, picking the weight vector :math:w \in \mathbf{R}^n_+ at each iteration to try and induce a sparse solution vector :math:x^\star. The details can be found in the Stanford EE364B lecture notes <https://stanford.edu/class/ee364b/lectures/l1_ext_slides.pdf>__.

The algorithm consists of initializing :math:w = 0 and repeating the two steps

.. math::

\begin{array}{ll} \mbox{maximize} & t - w^Tx\ \mbox{subject to} & t \leq Ax\ & 0 \leq x \leq 1 \ & \mathbf{1}^Tx = k, \end{array}

Let :math:w_i = \alpha/(\tau + x_i) \forall i

until we reach a Boolean solution. Here, :math:\alpha and :math:\tau are adjusted to promote a sparse solution. Typical choices would be :math:\alpha = 1 and :math:\tau = 10^{-4}.

Intuitively, the weight vector :math:w is incentivizing elements of :math:x which were close to zero in the last iteration towards zero in the next iteration.

Example

We create a parade route in the unit square :math:[0,1] \times [0,1] by generating points along a connected sequence of line segments. Possible guard locations are a set of randomly placed points in the unit square. Guards' coverage of points in the parade route is a function of distance between the two points. We also add buildings in the unit square to obstruct the guards' view. The guard has no effect on a point if his line of sight is obstructed. We generate the :math:A matrix for this problem instance below.

.. code:: python

import numpy as np
import matplotlib.pyplot as plt
import cvxpy as cp

def form_path(points,n):
    x, y = [], []
    pold = points[0]
    for p in points[1:]:
        x += list(np.linspace(pold[0],p[0],n))
        y += list(np.linspace(pold[1],p[1],n))
        pold = p

    path = np.array([x,y]).T
    return path

def form_grid(k):
    xs = list(np.linspace(0,1,k))
    ys = list(np.linspace(0,1,k))

    locations = []
    for x in xs:
        for y in ys:
            locations.append(np.array((x,y)))
    return np.array(locations).T

def guard_sets(k,num,noise):
    guard_set = []
    grid = form_grid(k)
    for i in range(num):
        pert = noise*np.random.randn(*grid.shape)
        guard_set.append( grid+pert )
    return np.hstack(guard_set)

def inRect(p,rect):
    x,y,w,h = rect
    return x <= p[0] and p[0] <= x + w and y <= p[1] and p[1] <= y + h

def remove_guards(guards,buildings):
    '''Remove guards inside buildings and outside unit square.'''
    outside = []
    for i, guard in enumerate(guards.T):
        inside = False
        for build in buildings:
            if inRect(guard,build):
                inside = True
                break
            if not inRect(guard,(0,0,1,1)):
                inside = True
                break
        if not inside:
            outside.append(i)
    
    return guards[:,outside]

def intersect(p1,p2,xmin,xmax,ymin,ymax):
    '''determine if a rectangle given by xy limits blocks the line of sight between p1 and p2'''

    block = False
    
    # if either point inside block
    for p in [p1,p1]:
        if xmin <= p[0] and p[0] <= xmax and ymin <= p[1] and p[1] <= ymax:
            return True
    
    # if the two points are equal at this stage, then they are outside the block
    if p1[0] == p2[0] and p1[1] == p2[1]:
        return False
    
    
    if p2[0] != p1[0]:
        for x in [xmin,xmax]:
            alpha = (x-p1[0])/(p2[0] - p1[0])
            y = p1[1] + alpha*(p2[1] - p1[1])

            if 0 <= alpha and alpha <= 1 and ymin <= y and y <= ymax:
                return True
            
    if p2[1] != p1[1]:
        for y in [ymin,ymax]:
            alpha = (y-p1[1])/(p2[1] - p1[1])
            x = p1[0] + alpha*(p2[0] - p1[0])

            if 0 <= alpha and alpha <= 1 and xmin <= x and x <= xmax:
                return True
        
    return False

def p_evade(x,y,r=.5,minval=.1):
    d = np.linalg.norm(x-y)
    if d > r:
        return 1
    return (1-minval)*d/r + minval

def get_guard_effects(path, guards, buildings, evade_func):
    guard_effects = []
    for guard in guards.T:
        guard_effect = []
        for p in path:
            prob = 1
            if not np.any([intersect(p,guard,x,x+w,y,y+h) for x,y,w,h in buildings]):
                prob = evade_func(p,guard)
            guard_effect.append(prob)
        guard_effects.append(guard_effect)
    return np.array(guard_effects).T

    locations = []
    for x in xs:
        for y in ys:
            point = np.array((x,y))
            detect_p = []
            for r in path:
                detect_p.append(p_evade(point,r,r=.5,m=0))
            locations.append((point,np.array(detect_p)))

.. code:: python

np.random.seed(0)

buildings = [(.1,.1,.4,.1),
             (.6,.1,.1,.4),
             (.1,.3,.4,.1),
             (.1,.5,.4,.1),
             (.4,.7,.4,.1),
             (.8,.1,.1,.3),
             (.8,.5,.2,.1),
             (.2,.7,.1,.3),
             (.0,.7,.1,.1),
             (.6,.9,.1,.1),
             (.9,.7,.1,.2)]

n = 10

points = [(.05,0),(.05,.25),(.55,.25),(.55,.6),(.75,.6),(.75,.05),(.95,.05), (.95,.45),(.75,.45), (.75,.65),(.85,.65),
          (.85,.85),(.35,.85),(.35,.65),(.15,.65),(.15,1)]

path = form_path(points,n)

g = guard_sets(12,4,.02)
g = remove_guards(g,buildings)

guard_effects = get_guard_effects(path, g, buildings, p_evade)

A = 1 - np.log(guard_effects)

fig = plt.figure(figsize=(10,10))
ax = plt.subplot(111,aspect='equal')
for x,y,w,h in buildings:
    rect = plt.Rectangle((x,y),w,h,fc='y',alpha=.3)
    ax.add_patch(rect)

ax.plot(path[:,0],path[:,1],'o')

ax.plot(g[0,:],g[1,:],'ro',alpha=.3)

.. parsed-literal::

[<matplotlib.lines.Line2D at 0x105d17cc0>]

.. image:: parade_route_files/parade_route_2_1.png

We perform the iterative algorithm below. At each step, we plot the vector :math:x, demonstrating that it becomes increasingly sparse at each iteration.

.. code:: python

num_guards = 12
tau = 1e-2

m,n = A.shape

w = np.zeros(n)

for i in range(3):
    x = cp.Variable(shape=n)
    t = cp.Variable(shape=1)

    objective = cp.Maximize(t - x.T*w)
    constr = [0 <=x, x <= 1, t <= A*x, cp.sum(x) == num_guards]
    cp.Problem(objective, constr).solve(verbose=False)
    x = np.array(x.value).flatten()
    w = 2/(tau+np.abs(x))
    fig = plt.figure(figsize=(5,5))
    ax = fig.add_subplot(111)
    ax.plot(x,'o')

xsol = x
print("final objective value: {}".format(objective.value))

.. parsed-literal::

final objective value: -10.27091799207174

.. image:: parade_route_files/parade_route_4_1.png

.. image:: parade_route_files/parade_route_4_2.png

.. image:: parade_route_files/parade_route_4_3.png

Below, we plot the final Boolean allocation. The blue dots represent the parade route. The red dots represent the possible guard placement locations. The green dots show the actual guard placements. Yellow rectangles are buildings which obstruct the guards' view.

.. code:: python

fig = plt.figure(figsize=(10,10))
ax = plt.subplot(111,aspect='equal')
for x,y,w,h in buildings:
    rect = plt.Rectangle((x,y), w, h, fc='y', alpha=.3)
    ax.add_patch(rect)

ax.plot(path[:,0], path[:,1], 'o')

ax.plot(g[0,:], g[1,:], 'ro', alpha=.3)
ax.plot(g[0,xsol > .5], g[1,xsol > .5], 'go', markersize=20, alpha=.5)

.. parsed-literal::

[<matplotlib.lines.Line2D at 0xb1853f1d0>]

.. image:: parade_route_files/parade_route_6_1.png

Optimal parade route

Optimal parade route

Topic references

Problem statement

Optimization formulation

Relaxation

Iterated weighted :math:\ell_1 heuristic

Example

Iterated weighted :math:`\ell_1` heuristic