Chapter 5
Dynamic programming
methods
5.1
Optimal solution for multi-target
track-ing
In the previous chapter the linear programming problem for multi-target tracking was introduced. To solve the problem it is necessary to verify whether a solution exists. The following theorems provide the conditions for the existence of a solution. Vanderbei (1997) has the following theorem: Theorem 1. For an arbitrary linear program, the following statements are true:
1. If there is no optimal solution, then the problem is either infeasible or unbounded.
2. If a feasible solution exists, then a basic feasible solution exists. 3. If an optimal solution exists, then a basic optimal solution exists.
Equation 4.7 subject to equation 4.8 cannot be solved with linear program-ming solvers such as the Simplex method due to the size of the problem. In addition, it is unlikely that the solution will be optimal. Rewrite the linear
programming problem (equation 4.7 subjected to equation 4.8) as maximize X t,i log pt i 1− pt i X jǫN(i) fi,jt subject to f ≥ 0 Cf ≤ [1, ..., 1, 0, ..., 0]T
with f a flow vector and C a constraint matrix. The constraint matrix is totally unimodular as will be shown, and the right hand side consists of integers which will ensure an optimal solution. A matrix is unimodular when the determinant is 0, 1 or −1 (Van Mieghem, 2011, 17; Berclaz et al., 2010). The set{x : Cx ≤ b, x ≥ 0} is the polytope of feasible solutions. The optimal solution is on the boundaries. The polytope, also known as convex closure or convex hull, can be defined as follows: For a subset X ⊆ Rnarbitrary chosen
the intersection of all the subsets in Rnwhich contain X and are convex, is a
convex hull (Edelberg, 1970). The convex hull, H(X), is the smallest subset containing X in Rn which is convex. Further, H(X) = X when X is also
convex. One of the vertices is the optimal solution unless there is more than one optimal solution. In the case of multi optimal solutions at least two of the solutions are vertices of the polytope.
The constraint matrix C is unimodular. This can be proven with the use of the following two theorems. The two theories are adapted from Berclaz et al. (2010).
Theorem 2. A matrix A = {aij} ∈ Zm×n is totally unimodular if and
only if for every subset R ⊆ {1, ..., m} of rows, there exists a partition R = R1∪ R2, R1∩ R2 =∅ such that ∀j = 1, ..., n X i∈R1 aij − X i∈R2 aij ∈ {0, −1, 1}.
Theorem 3. If C is a totally unimodular matrix, then the vertices of the polytope {x : Cx ≤ b, x ≥ 0} are integral, that is, have integer coordinates, for any integral vector b.
Divide the rows of matrix C into two subsets U1 and U2, where U1 consists
X
jǫN(i)
fi,jt 61
and U2 consists of the last two constraints
X jǫN(i) ft i,j− X k:iǫN(k) fk,it−1 ≤ 0 X j∈N(vsource) fvtsource,j− X k:vsink∈N(k) fk,vt sink ≤ 0 (Jongen et al., 2004, 261).
Without loss of generality the columns can be ordered according to time. At i the node has the value−1 if jǫ N (i) or else 0. Define f = [f1
1,1,· · · , fK,KT ] T,
with T the last time frame, which contains all the flow values. The constraint can be rewritten as
C· f ≤ [1, · · · , 1, 0, · · · , 0]T. (5.1) Construct with R arbitrarily chosen rows a matrix CR, which is a sub-matrix
of C and let R1 = U1∩ R and R2 = U2∩ R. Then for each column there will
be 8 possibilities; thus for column j the possibilities are {aij|i ∈ R1} {aij|i ∈ R2} Pi∈R1aij − P i∈R2aij {0, ..., 0, 1} {0, ..., 0} 1 {0, ..., 0, 1} {0, ..., 0, 1} 0 {0, ..., 0, 1} {0, ..., 0, −1} 2 {0, ..., 0, 1} {0, ..., 0, 1, −1} 1 {0, ..., 0} {0, ..., 0} 0 {0, ..., 0} {0, ..., 0, 1} -1 {0, ..., 0} {0, ..., 0, −1} 1 {0, ..., 0} {0, ..., 0, 1, −1} 0 (Berclaz et al., 2010).
Except for row three which contains the value 2 in the table above, theorem 2 is satisfied. Berclaz et al., (2010) suggested the idea of moving the rows in R1 which contain a non-zero element to R2; then theorem 2 will be satisfied.
5.2
Solution using k -shortest paths
Having concluded that a solution for the formulated multi-target tracking problem exists, different approaches to determine the solution are now dis-cussed.
The shortest path problem can be defined as follows: For a graph G = (V, E), a digraph, where s is a distinguished vertex and s∈ V and l : E → ℜ+ is the
length function (Jongen et al., 2004, 233), determine a directed (s− v)-path with minimum length with v the source node (Suurballe & Tarjan 1984). Therefore the k -shortest path problem is based on this principle of finding k such minimum length paths where k can be any integer value.
The reformulation of the linear programming problem, equation 4.7 subjected to equation 4.8, as a k -shortest path problem on a directed acyclic graph will be as follows (Berclaz et al., 2010): Using figure 4.5, which represents a directed graph G with the pair of nodes vsource and vsink the aim is to determine the k node disjoint paths {p1, p2,· · · pk} with minimum total cost.
In figure 4.5 each path represents a single object’s flow. The cost of each directed edge (et i,j) is c(et i,j) =− log pt i 1− pt i (5.2) with the flow for one time interval (t to t+1) being from location i to location j (Berclaz et al., 2010). The k -shortest paths optimal solution will then be
f∗ = argmin f ∈η X t,i c(et i,j) X j∈N(i) fi,jt (5.3)
with η a set of feasible solutions of equation 4.7, which also satisfy equations 4.8 (Berclaz et al., 2010). For k arbitrarily chosen, η will contain all the k paths between vsource and vsink.
5.2.1
Remark
The multi-target tracking can be formulated as a linear programming prob-lem as demonstrated in equations 4.7 and 4.8. The cost of each directed edge can be determined by means of equation 4.8 (see equation 5.2). For equation 5.2 if
0 < p t i 1− pt i <1, (5.4)
then the cost of the directed edge is positive and if pt
i
1− pt i
>1, (5.5)
then the cost of the directed edge is negative. Thus for Dijkstra’s algorithm inequality 5.4 must hold. For the Floyd-Warshall algorithm, both inequalities 5.4 and 5.5 are applicable. Thus for inequality 5.4 to hold
0 < pti <
1
2, (5.6)
and for inequality 5.5
pti > 1 2.
5.3
Algorithm to find k -shortest paths
For the k -shortest path problem, dynamic programming methods such as Dijkstra’s and Floyd-Warshall algorithms can be used to solve equation 5.3. The Dijkstra algorithm only allows non-negative edges, whereas Floyd-Warshall allows non-negative and negative edges. These two algorithms will now be discussed.
5.3.1
Dijkstra’s algorithm
In 1959 Dijkstra solved this problem with the following algorithm (Lipschutz & Lipson 2007, 219). For a cycle-free graph G:
1. For a graph similar to figure 4.5, the aim of the algorithm is to deter-mine the shortest path from the source node, u, to the sink node, w, w using the notation:
• Let l(v) be the length from any node to v and p(v) be the path from the same node to v.
• Initially let p(u) = u and l(u) = 0. 2. For each step examine an edge e = (v′′
, v′
); if the length from v′′
to v′
is a, then calculate l(v′′
) + a. (a) If l(v′′
) + a < l(v′
) then a shortest path is determined and then l(v′) will be replaced by l(v′′) + a and p(v′) with p(v′′)v′.
(b) Else l(v′
) and p(v′
) remain unchanged.
3. The algorithm ends when p(w) with l(w), the smallest possible value, is determined.
Example
As adapted from Lipschutz & Lipson (2007, 218-220). Refer to figure 5.1 In matrix form it will be
W = 0 4 6 2 ∞ ∞ ∞ ∞ ∞ 0 3 ∞ 4 ∞ ∞ ∞ ∞ ∞ 0 ∞ ∞ 2 1 ∞ ∞ ∞ 3 0 ∞ ∞ 5 ∞ ∞ ∞ ∞ ∞ 0 2 ∞ 3 ∞ ∞ ∞ ∞ ∞ 0 ∞ 3 ∞ ∞ ∞ ∞ ∞ ∞ 0 3 ∞ ∞ ∞ ∞ ∞ ∞ ∞ 0 .
From the source (node u) there are three successive vertices, nodes x, y and
z with
l(x) = 4 p(x) = ux
l(y) = 6 p(y) = uy
Figure 5.1: Tree graph for specific example (Adapted from Lipschutz & Lip-son 2007, 218)
• First only use the branches from x :
– Node x has two successors, r and y, from which l(r) = 4 + 4 = 8 with p(r) = p(x)r = uxr l1(y) = 4 + 3 = 7 with p1(y) = p(x)y = uxy
Since l1(y) > l(y), l(y) and p(y) remain unchanged.
– Node r has the two possible successors, s and w, with l(s) = 4 + 4 + 2 = 10, p(s) = p(r)s = uxrs and
l(w) = 4 + 4 + 3 = 11, p(w) = p(r)w = uxrw. – Node s has only one successor, w, with length
l1(w) = 4 + 4 + 2 + 3 = 13
and path
p1(w) = p(s)w = uxrsw.
Since 13 = l1(w) > l(w) = 11, l(w) and p(w) stay unchanged.
Since 7 = l1(y) > l(y) = 6 there is no need to determine the length
of the rest of the branches. • Branches from y:
– Two successors, s and t, are possible from y, with l1(s) = 6 + 2 = 8, p1(s) = p(y)s = uys
l(t) = 6 + 1 = 7, p(t) = p(y) = uyt.
Since 8 = l1(s) < l(s) = 10, the shorter path p1(s) with length
l1(s) will be used.
– Form s there is only one successor, w
l2(w) = 6 + 2 + 3 = 11, p2(w) = uysw
– The only successor of t is w, with
l3(w) = 6 + 1 + 3 = 10, p3(w) = uytw.
Since 10 = l3(w) < l(w) = 11, further on l3(w) and p3(w) will be
used.
• Branches from z :
– From z there are two successors, y and t, the length and paths being
l2(y) = 2 + 3 = 5, p2(y) = uzy
and
l1(t) = 2 + 5 = 7, p1(t) = uzt
Since 5 = l2(y) < l(y) = 6, the path and associated length l2(y)
and p2(y) will be used. Since l1(t) = l(t) = 7, both are kept for
now.
– From y there are two successors, s and t, resulting in l2(s) = 2 + 3 + 2 = 7, p2(s) = p(y)s = uzys
and
l2(t) = 2 + 3 + 1 = 6, p2(t) = p(y) = uzyt
As previously l2(s) and p2(s) is used, since 7 = l2(s) < l1(s) = 8.
Also since 6 = l2(t) < l1(t) = 7, further on l2(t) and p2(t) will be
used.
– From s there is only one successor, w, for which
l4(w) = 2 + 3 + 2 + 3 = 10, p4(w) = uzysw.
Since l4(w) = l3(w) = 10, l3(w) and p3(w) will be unchanged.
– From t there is only one successor, w
l5(w) = 2 + 3 + 1 + 3 = 9, p5(w) = uzytw
and since 9 = l5(w) < l3(w) = 10, l5(w) and p5(w) will be used.
Since 6 = l2(t) < l1(t) = 7, there is no need to calculate the length
of the rest of the branches.
Therefore, from the source u to the sink w the shortest path is p5(w) = uzytw with length l5(w) = 9.
Results with Matlab
See chapter 7. The input for the Matlab code 7.1.1 is the weight matrix
W =
0 4 6 2 inf inf inf inf
inf 0 3 inf 4 inf inf inf
inf inf 0 inf inf 2 1 inf
inf inf 3 0 inf inf 5 inf
inf inf inf inf 0 2 inf 3
inf inf inf inf inf 0 inf 3
inf inf inf inf inf inf 0 3
inf inf inf inf inf inf inf 0. (5.7)
With the command
[spcost] = dijkstra9(W,1,8)
Matlab performs the calculations as shown above. Solution:
spcost = 9
path = 9 1 4 1 2 3 3 7
With
Node number Node symbol
1 u 2 x 3 y 4 z 5 r 6 s 7 t 8 w
the results are the same, the shortest path is 1→ 4 → 3 → 7 → 8 with cost 9.
5.3.2
Floyd-Warshall Algorithm
For each iteration j the matrices Wj and Pj are calculated with W the
weight matrix and initially P0 = P and W = W0 (Balakrishnan, 1997, 120).
Further, let wj(u, v) be the element of the matrix Wj, which represents the
edge between the node u and v and pj(u, v) the element of Pj and the path
between node u and v with u the adjacent of v. If there are no edges between the nodes then the weight element is ∞.
The triangle inequality for vectors states that for two vectors a and b |a + b| ≤ |a| + |b|
with |a| and |b| the length of a and b respectively (Clapham & Nicholson, 2009). Thus for the weight matrix the elements for each iteration can be calculated as follows:
wj−1(u, v)≤ wj−1(u, j) + wj−1(j, v) (5.8)
(Balakrishnan, 1997, 120-121).
If inequality 5.8 is satisfied, then wj(u, v) = wj−1(u, v) and pj(u, v) = pj−1(u, v);
else wj(u, v) = wj−1(u, j) + wj−1(j, v) and pj(u, v) = pj−1(u, j). The input
is an n× n weight matrix with the output being the new weight and path matrices. The new weight matrix represents the shortest distance and the new path matrix the shortest path.
Example
As adapted from Balakrishnan (1997, 120). For figure 5.2 the matrices will be
W = 0 4 −3 ∞ −3 0 −7 ∞ ∞ 10 0 3 5 6 6 0 (5.9)
Figure 5.2: Network (Balakrishnan, 1997, 120) and P = 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 . (5.10)
Then for iteration j = 0:
• For w1(1, 1) test if the inequality 5.8 holds. Now
w0(1, 1) + w0(1, 1) = 0 + 0,
thus the inequality holds since w0(1, 1) = w0(1, 1) + w0(1, 1) and
there-fore w1(1, 1) = w0(1, 1) and p1(1, 1) = p0(1, 1).
• Similarly for w1(1, 2),
w0(1, 1) + w0(1, 2) = 0 + 4 = 4,
thus the inequality holds since w0(1, 2) = w0(1, 1) + w0(1, 2) and
• For w1(1, 3),
w0(1, 1) + w0(1, 3) = 0− 3 = −3,
thus the inequality holds since w0(1, 3) = w0(1, 1) + w0(1, 3) and
there-fore w1(1, 3) = w0(1, 3) and p1(1, 3) = p0(1, 3).
• For the rest the calculations will be the same except for w1(4, 3)
w0(4, 1) + w0(1, 3) = 5− 3 = 2
the inequality does not hold since w0(4, 3) > w0(4, 1) + w0(1, 3) and
therefore w1(4, 3) = w0(4, 1) + w0(1, 3) = 2 and p1(4, 3) = p0(4, 1). Thus W1 = 0 4 −3 ∞ −3 0 −7 ∞ ∞ 10 0 3 5 6 2 0 and P1 = 1 2 3 4 1 2 3 4 1 2 3 4 1 2 1 4 (Balakrishnan, 1997, 120). For iteration j = 2: W2 = 0 4 −3 ∞ −3 0 −7 ∞ 7 10 0 3 3 6 −1 0 and
P2 = 1 2 3 4 1 2 3 4 2 2 3 4 2 2 2 4 (Balakrishnan, 1997, 120-121). For j = 3 W3 = 0 4 −3 0 −3 0 −7 −4 7 10 0 3 3 6 −1 0 and P3 = 1 2 3 3 1 2 3 3 2 2 3 4 2 2 2 4 (Balakrishnan, 1997, 120-121). For j = 4 W4 = 0 4 −3 0 −3 0 −7 −4 6 9 0 3 3 6 −1 0 and P4 = 1 2 3 3 1 2 3 3 4 4 3 4 2 2 2 4
(Balakrishnan, 1997, 120-121).
Matrix W4 is the shortest distance matrix and matrix P4 is the shortest path
matrix and therefore the shortest paths are 1→ 1 1→ 2 1→ 3 1→ 3 → 4 2→ 1 2→ 2 2→ 3 2→ 3 → 4 3→ 4 → 2 → 1 3 → 4 → 2 3→ 3 3→ 4 4→ 2 → 1 4→ 2 4→ 2 → 3 4→ 4
The (1,1) entries in the two matrices are respectively 0 and 1; therefore the shortest distance from node one to one is zero and the path is simply 1 → 1. The (1, 4) entries are respectively 0 and 3; therefore the shortest path is not directly from 1→ 4, but via node 3. To determine the rest of the path from node 3; use (3, 4), which indicates in P4 that from node 3 we go to node
4. Thus the total shortest path from node 1 to 4 is 1 → 3 → 4 with total distance 0.
Results with Matlab
See chapter 7. Using Matlab code 7.1.2 the input is w =
0 4 -3 Inf
-3 0 -7 Inf
Inf 10 0 3
5 6 6 0
which is the weight matrix. With the command
[new, path] = floyd(w)
Matlab, with the use of for loops, performs the same calculations as indicated above. Matlab only displays the final results.
weight = 0 4 −3 0 −3 0 −7 −4 6 9 0 3 3 6 −1 0 path = 1 2 3 3 1 2 3 3 4 4 3 4 2 2 2 4
