A Proximal Point Algorithm for Location Problems

H.Benker  A.Hamel  J.Spitzner  Chr.Tammer

1 Introduction

Many authors have studied generalizations of the well-known Fermat-Weber problem from the theoretical and computational point of view (see [1], [2], [3], [4], [5]). The aim of this paper is to give a MAPLE-program (MAPLE-worksheet) for the numerical solution of the following location problem in the plane:

f(x) = ∑_{i=1}^{n} λ_{i} ∥x − a^{i}∥_{i} → min_{x∈D},  (1)

where a^{i} and x = (x_{1}, x_{2}) ∈ R^{2}, λ_{i} are weighting factors ≥ 0, (i = 1, ..., n), R^{2} is the twodimensional Euclidean space, D_{j} ⊆ R^{2} (j = 1, ..., m) are closed and convex sets, D_{0} ⊆ R^{2} is a linear subspace, D = ⋂_{j=0}^{m} D_{j}, ∥...∥_{i} denotes a

norm in the space R^{2} and B_{i}^{0} the unit ball with respect to the dual norm ∥...∥_{i}^{∗}. In the following we assume that anyone of constraint qualifications (generalized Slater condition or stability) is fulfilled. Furthermore, P_{D}(e) denotes the projection of e onto a set D. N_{D_{j}}(x^{0}) denotes the normal cone to D_{j} at x^{0} ∈ R^{2}. In chapter 2 we consider optimality conditions for the problem (1) which we apply in chapter 3 in order to derive a proximal point algorithm. A MAPLEprogram for this algorithm is presented in chapter 4. Finally, in chapter 5 the given MAPLE-program is applied in order to solve a location problem in town planning.

2 Optimality Conditions

The given approach follows the recent papers of Michelot and Lefebvre [5] and Idrissi, Lefebvre and Michelot [3] (see also [1]).

Using the indicator functions χ_{D_{j}} of D_{j} defined by χ_{D_{j}}(x) = 0 if x ∈ D_{j} and χ_{D_{j}}(x) = +∞ otherwise the problem (1) is equivalent to the following unconstrained minimization problem

F (x) = ∑_{i=1}^{n} λ_{i} ∥x − a^{i}∥_{i} + ∑_{j=0}^{m} χ_{D_{j}}(x) → min_{x∈R^{2}} .  (2)

The functional F (x) is convex under the given assumptions. Therefore, the subdifferential-condition

0 ∈ ∂F (x^{0})  (3)

is necessary and sufficient for the optimality of x^{0}. Applying the rule for sums of subdifferentials and the known formula for the subdifferential of a indicator function (∂χ_{D_{j}}(·) = N_{D_{j}}(·)) on the functional F in (2) we can write the optimality conditions (3) equivalently in the form (4) - (6):

q_{i} ∈ λ_{i} · ∂ ∥x^{0} − a^{i}∥ ,  i = 1, 2, ..., n,  (4)
r_{j} ∈ N_{D_{j}}(x^{0}) ,  j = 1, 2, ..., m,  (5)
∑_{i=1}^{n} q_{i} + ∑_{j=1}^{m} r_{j} ∈ D_{0}^{⊥}.  (6)

In order to reformulate the optimality conditions (4) - (6) in a more practical way we introduce:

• The space E := (R^{2})^{n+m} with the inner product ∑_{i=1}^{n+m} u_{i} · v_{i}

for u = (u_{1}, ..., u_{n+m}) ∈ E and v = (v_{1}, ..., v_{n+m}) ∈ E.

• The subspaces

A := {y ∈ E | y = (y_{1}, ..., y_{n+m}) ; y_{j} = x , j = 1, ..., n + m; x ∈ D_{0}} ,

B := { p ∈ E  |  p = (p_{1}, ..., p_{n+m}) ,  ∑_{i=1}^{n+m} p_{i} ∈ D_{0}^{⊥}} .

• The operator T defined on E by

T (y) := ∏_{i=1}^{n+m} T^{i}(x) = (T^{1}(x), ..., T^{n+m}(x)),  (7)

where

T^{i}(x) := λ_{i} · ∂ ∥x − a^{i}∥ (i = 1, ..., n) and T^{n+j}(x) := N_{D_{j}}(x)

(j = 1, ..., m).

With this notations the optimality conditions (4)-(6) can be rewritten as:

Find (y^{0}, p^{0}) ∈ A × B such that p^{0} ∈ T (y^{0}).  (8)

Remark : It is easy to see that the operator T, defined in (7), realizes a maximal monotone multifunction since T is composed by subdifferentials. For a numerical solution we need the condition (8) in the form 0 ∈ T (y^{0}) − p^{0} but the new multifunction T̃(y, p) := T (y)−p is not maximal monotone. Therefore, using an idea of Spingarn [6], we obtain the algorithm given in chapter 3.

3 Application of the Proximal Point Algorithm

Now we can apply the proximal point algorithm f or solving (8) (compare Michelot/Lefebvre [5] and Idrissi/Lefebvre/Michelot [3]). Problem (8) is equivalent to:

Find z^{0} ∈ E such that 0 ∈ T_{A}(z^{0}),

where T_{A} is the partial inverse of T with respect to A (see [1],[2]).

For this problem the proximal point algorithm has the following form

z^{k+1} = (I + c_{k}T_{A})^{−1}z^{k},  k = 1, 2, ...,  (9)

where the starting point z^{1} is arbitrary and {c_{k}} is a sequence of positive real numbers . If {c_{k}} is bounded away from zero (we set c_{k} = 1) then either {z^{k}} converges weakly to a zero of T_{A} or ∥z^{k}∥ → ∞ and T_{A} has no zeros. The algorithm can be transformed in the following final form:

Algorithm :

• Choose the starting points

x^{1} ∈ D_{0} and p^{1} ∈ E  with ∑_{i=1}^{n+m} p_{i}^{1} ∈ D_{0}^{⊥}.  (10)

• Compute x^{k+1} and p^{k+1} from

x^{k} − x^{k+1} + p_{i}^{k+1} = λ_{i} · P_{B^{0}}(\frac{x^{k} + p_{i}^{k} − a^{i}}{λ_{i}}),  i = 1, ..., n, (11)

x^{k+1} + p_{n+j}^{k} − p_{n+j}^{k+1} = P_{D_{j}}(x^{k} + p_{n+j}^{k}),  j = 1, ..., m,  (12)

with

∑_{i=1}^{n+m} p_{i}^{k+1} ∈ D_{0}^{⊥}  and x^{k+1} ∈ D_{0}.  (13)

That means that we must first calculate the projections on B^{0} and D_{j} in the right-hand side of the equations (11) and (12). Then we have to solve the linear equations (11),(12) such that the relations (13) are satisfied. The solution of the linear equations (11),(12) is very simple.

• The algorithm is finished if the difference of two iterations is less than a given error bound eps.

Remark : We remark that the general convergence results for Spingarn’s algorithm (compare [6]) also hold for our application of the proximal point algorithm.

4 MAPLE-program

In the following we present a MAPLE-program (using MAPLE V Release 3) for the algorithm given in chapter 3:

• In this MAPLE-program we can choose between the Euclidean norm, the maximum norm and the Lebesgue norm in the twodimensional space R^{2} .

• The set of constraints D can be formed by intersection of halfplanes given by linear inequalities of the form α_{j} x_{1} + β_{j} x_{2} + γ_{j} ≤ 0.

MAPLE-Subroutine for the given Proximal-Point-Algorithm in order to compute the solution of the location problem in the plane (1):

• The subroutine needs the following parameters:

1. The location vectors a^{i} and the factors lambda_{i} (= λ_{i}) in the form of the list parameter := [[a^{1}, lambda_{1}], ...,[a^{n}, lambda_{n}]]; for example parameter := [[1, 0, 1], [1, 2, 1]]; if a^{1} := [1, 0], a^{2} := [1, 2] and lambda_{1} = lambda_{2} := 1.

2. The norms used in the cost functional (1) in the form of the list norm := [norm_{1}, ..., norm_{n}]; where we choose norm_{i} = 1 if we have the Lebesgue norm, norm_{i} = 2 if we have the Euclidean norm, norm_{i} = 3 if we have the maximum norm.

3. The constraints D given by linear inequalities α_{j}x_{1}+β_{j}x_{2}+γ_{j} ≤ 0 in the form of the list of the coefficients α_{j},β_{j},γ_{j} : meck := [[α_{1}, β_{1}, γ_{1}], ..., [α_{j}, β_{j}, γ_{j}]]; for instance meck := [[1, 0, −1], [−1, 0, −1], [0, 1, −1], [0, −1, −1]].

4. Error bound eps for the result, for instance eps := 10^{−10} .

• The numerical result erg obtained by the algorithm is given in a vector with three components

erg = [i, x_{1}, x_{2}],

where i is the number of iterations and x_{1}, x_{2} are the components of the computed result.

• The starting point for the iteration is zero.

• Call of the subroutine: s(parameter, norm, meck, eps);

• subroutine: s := proc(param, norm, dger, eps)

local i, par, n, dg, m, parama, paraml, ppv1, ppv2, l, differ, xyks, erg; with(linalg); n := vectdim(param); parama := [[param[1][1], param[1][2]]]; for i from 2 to n do par:=[param[i][1],param[i][2]]; parama := [op(parama), par] od; paraml := [param[1][3]]; for i from 2 to n do paraml:=[op(paraml), param[i][3]] od; m := vectdim(dger); xyks := [0,0]; ppv1:=[xyks]; for i from 2 to n+m+1 do ppv1:=[op(ppv1), xyks] od; ppv2 := ppa(parama, paraml, norm, ppv1, dger, n, m); differ := max(abs(ppv1[1][1]-ppv2[1][1]), abs(ppv1[1][2]-ppv2[1][2])); for i from 2 to n+m+1 do differ := max(differ, abs(ppv1[i][1]-ppv2[i][1]), abs(ppv1[i][2]-ppv2[i][2])) od; l := 1; while differ < eps do ppv1:=ppv2; ppv2 := ppa(parama, paraml, norm, ppv1, dger, n, m); l := l+1; differ := max(abs(ppv1[1][1]-ppv2[1][1]), abs(ppv1[1][2]-ppv2[1][2])); for i from 2 to n+m+1 do differ := max(differ, abs(ppv1[i][1]-ppv2[i][1]),

abs(ppv1[i][2]-ppv2[i][2])) od od; erg := [l, evalf(ppv2[1][1]), evalf(ppv2[1][2])]; erg end

• Proximal point algorithm:

ppa := proc(parama, paraml, norm, xyk, dger, n, m) local x, y, he, hep, hev, lv, grd, s, lot, t, lsg, erg, yh, xh, xmat, ymat, xykp, xykph, i, j; xmat := [xyk[1][1], xyk[1][2]]; ymat := [[xyk[2][1], xyk[2][2]]]; for i from 3 to n+m+1 do ymat:=[op(ymat), [xyk[i][1], xyk[i][2]]] od; x:=ymat; y:=ymat; xykp:=xyk; for i from 1 to n do yh := add(xmat, parama[i], 1, -1); yh := evalm(add(yh, ymat[i], 1, 1) / paraml[i]); yh := evalm(paraml[i] * projki(yh,norm[i])); y := subsop(evalm(i=yh), y) od; for i from 1 to n do xh := add(xmat, ymat[i], 1, 1); xh := add(xh, y[i], 1, -1); x := subsop(evalm(i=xh), x) od; for j from 1 to m do xh := add(xmat, ymat[n+j], 1, 1); x := subsop(evalm(n+j=xh), x); he := parfo(dger[j]); hep := [he[1], he[2]]; hev := [he[3], he[4]]; lv := [hev[2], -hev[1]]; grd := add(hep, hev, 1, s); lot := add(x[n+1], lv, 1, t); lsg := solve(grd[1]=lot[1], grd[2]=lot[2], s,t); erg := subs(lsg, evalm(grd)); if 0 < (dger[j][1]*x[n+j][1] + dger[j][2]*x[n+j][2] + dger[j][3]) then xh:=erg else xh:=evalm(xh) fi; x := subsop(evalm(n+j=xh),

x) od; for j from 1 to m do yh := add(xmat, ymat[n+j], 1, 1); yh := add(yh, x[n+j], 1, -1); y := subsop(evalm(n+j=yh), y) od; xykph := x[1]; for i from 2 to n+m do xykph:=add(xykph, x[i], 1, 1) od; xykph := evalm(xykph / (n+m)); xykp := subsop(evalm(1=xykph), xykp); yh := y[1]; for i from 2 to n+m do yh:=add(yh, y[i], 1, 1) od; yh := evalm(yh / (n+m)); for i from 1 to n+m do xykph:=add(y[i], yh, 1, -1); xykp := subsop(evalm(i+1=xykph), xykp) od; xykp end

• Parametrization of the feasible set:

parfo := proc(ger) local x0, x1, y0, y1, gerade,gp, gv; if ger[2] <> 0 then x0:=0 else x0:=evalf(-ger[3]/ger[1]) fi; if ger[2] <> 0 then x1:=1 else x1:=x0 fi; if ger[2] <> 0 then y0:=evalf(-ger[3]/ger[2]) else y0:=0 fi; if ger[2] <> 0 then y1:=evalf(-ger[3]/ger[2]-ger[1]/ger[2]) else y1:=1 fi; gp := [x0,y0]; gv := [x1-x0, y1-y0]; gerade := [op(gp), op(gv)]; gerade end

• Projection onto the unit sphere for the corresponding norms: projki := proc(ny, nrn) local betr, ney, ney1, ney2, neyh; ney := ny;if nrn=1 then ney1 := ny[1] + ny[2]; ney2 := - ny[1] + ny[2]; ney := [ney1, ney2]; neyh := projki(ney,3); ney1 := (neyh[1]

- neyh[2]) /2; ney2 := (neyh[1] + neyh[2]) /2; ney := [ney1, ney2]; elif nrn=2 then betr := evalf(sqrt(ny[1]ˆ 2+ny[2]ˆ 2));if 1 < betr then ney:=evalm(ny/betr) else ney:=evalm(ny) fi; else if abs(ny[1]) <= 1 then ney1:=ny[1] else ney1:=signum(ny[1]) fi; if abs(ny[2])<=1 then ney2:=ny[2] else ney2:=signum(ny[2]) fi; ney := [ney1, ney2] fi; ney end;

5 Example

We will study the location problem in town planning to determine a location for a children’s playground in a new part of a town (i.e., the town Halle (Saale) in Germany) and apply the algorithm from Section 4 in order to solve this problem. In a certain district we have given facilities of blocks of flats, schools, daycare centers and railway stations. The blocks of flats are located in a^{1} = (3, 5; 3, 5), a^{2} = (4, 5; 4), a^{3} = (4; 2, 5), a^{4} = (5; 1, 5), a^{5} = (6; 1), a^{6} = (6, 5; 1, 5), a^{7} = (6; 4), a^{8} = (7; 3, 5), a^{9} = (8; 3), a^{10} = (9, 5; 3), a^{11} = (10, 5; 3, 5), a^{12} = (10, 5; 5), a^{13} = (6; 5, 5), a^{14} = (6, 5; 6, 5), a^{15} = (11; 6, 5), a^{16} = (10; 7, 5), a^{17} = (9; 8), the schools are located in a^{18} = (2; 4, 5), a^{19} = (2; 6), the day nurseries are located in a^{20} = (3; 4, 5), a^{21} = (6, 3) and railway stations are located in a^{22} = (8; 6, 5), a^{23} = (3; 13, 5). We note that in the blocks located in a^{3} and a^{17} a great number of children at the age of 5 up to 12 live. Now the problem is to determine a location for a children’s playground in such a way that the distances to the given facilities are minimal in a certain sense. It is possible to formulate the problem as a real-valued location problem (Fermat-Weber-Problem) which is the problem to determine a location x of a new facility such that the weighted sum of the distances between the given facilities a^{i} (i = 1, ..., 23) and x is minimal (see problem (1)).

If we choose the weights λ_{i} = 1 for all i=1,...,23, then in order to apply the MAPLE-program given in chapter 4 we must define:

• parameter:=[[3.5,3.5,1],[4.5,4,1],[4,2.5,1],[5,1.5,1],[6,1,1],[6.5,1.5,1], [6,4,1],[7,3.5,1],[8,3,1],[9.5,3,1],[10.5,3.5,1],[10.5,5,1],[6,5.5,1], [6.5,6.5,1],[11,6.5,1],[10,7.5,1],[9,8,1],[2,4.5,1],[2,6,1],[3,4.5,1], [6,3,1],[8,6.5,1],[3,13.5,1]]; for the location vectors,

• meck:=[[-1,0,2],[1,0,-9],[0,-1,5],[0,1,-8]]; for instance for the constraints,

• eps := 10^{−3}; for instance for the error bound.

(1) If we take for all norms the maximum norm we must define NORM:=[3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3]; The application of the subroutine s(parameter,NORM,meck,eps); yields the result [137, 6.000004245, 4.999026097].

(2) If we take for all norms the Euclidean norm we must define NORM:=[2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2]; and we get the result [80, 6.322029126, 4.999000119].

(3) If we take for all norms the Lebesgue norm we must define NORM:=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]; and we get the result [179, 6.499083675, 4.999680615].

(4) If we take different norms, i.e., NORM:=[2,1,2,2,2,2,3,2,2,2,1,2,2,3,2,2,2,1,2,2,3,2,2]; we get the result [96, 6.130685912, 5.000654600].

If we choose the different weights λ_{3} = 5 , λ_{17} = 5 and λ_{i} = 1 for the other i, we must define:

• parameter:=[[3.5,3.5,1],[4.5,4,1],[4,2.5,5],[5,1.5,1],[6,1,1],[6.5,1.5,1], [6,4,1],[7,3.5,1],[8,3,1],[9.5,3,1],[10.5,3.5,1],[10.5,5,1],[6,5.5,1], [6.5,6.5,1],[11,6.5,1],[10,7.5,1],[9,8,5],[2,4.5,1],[2,6,1],[3,4.5,1], [6,3,1],[8,6.5,1],[3,13.5,1]]; for the location vectors, • meck:=[[-1,0,2],[1,0,-9],[0,-1,5],[0,1,-8]];

for instance for the constraints,

• eps := 10^{−3}; for instance for the error bound.

(5) If we take for all norms the Euclidean norm we must define NORM:=[2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2]; and we get the result [66, 6.314127419, 5.000998134].

Remark: The choice of the concrete norm in the location problem depends on the courses of the roads in the city. Typically for this problem is the weighted sum of the so called ”Manhattan norm” (maximum norm) in the objective function of the location problem.

6 Comments

In the examples (1) - (4) given above we have choosen the weights λ_{i} = 1 for all i = 1, ..., n. This means that all distances between the location for the children’s playground and the given facilities a^{i} (i = 1, ..., n) are equally important. In practice, it is often of interest to use different weights in the location problem. In our example this is the case for the distances between the location for the children’s playground and the blocks of flats. Here the weights for blocks of flats with a great number of children at the age of 5 up to 12 (in our case the blocks located in a^{3} and a^{17}) should be choosen greater than for other blocks. In the example (5) we have choosen λ_{3} = 5 , λ_{17} = 5 and λ_{i} = 1 for the other i. Moreover, it is also practically important to use different norms in the items of the objective function of the location problem. This is the case when the roads in the city have different courses.

References

[1] Benker, H.; Hamel, A.; Tammer, Chr.(1996): A Proximal Point Algorithm for Control Approximation Problems. ZOR-Mathematical Methods of Operations Research, Vol.43, Issue 3 (to appear).

[2] Benker, H.; Hamel, A.; Tammer, Chr.(1995): An Algorithm for Vectorial Control Approximation Problems. Proceedings of the Twelfth International Conference on Multiple Criteria Decision Making Hagen, Springer Verlag (to appear).

[3] Idrissi, H.; Lefebvre, O.; Michelot, C.(1988): A Primal-Dual Algorithm for a constrained Fermat-Weber Problem Involving Mixed Norms. Recherche operationnelle / Operations Research 22 (4), 313-330.

[4] Kuhpfahl, I.; Patz, R.; Tammer, Ch.(1995): Location problems in town planning. Wiss. Zeitschrift Universitaet Kaiserslautern (to appear).

[5] Michelot, C.; Lefebvre, O.(1987): A Primal-dual Algorithm for the Fermat-Weber Problem Involving Mixed Gauges. Mathematical Programming 39, 319-335.

[6] Spingarn, J.E.(1983): Partial inverse of a monontone operator. Applied Mathematics and Optimization (10), 247-265.

[7] Spingarn, J.E.(1985): Applications of the Method of partial inverses to convex programming: decomposition. Mathematical Programming (32), 199-223.