to appear in Journal of Global Optimization, 2011

**A continuation approach for the capacitated multi-facility Weber** **problem based on nonlinear SOCP reformulation**

Jein-Shan Chen ^{1}
Department of Mathematics
National Taiwan Normal University

Taipei 11677, Taiwan

Shaohua Pan^{2}

School of Mathematical Sciences South China University of Technology

Guangzhou 510640, China

Chun-Hsu Ko ^{3}

Department of Electrical Engineering I-Shou University

Kaohsiung 840, Taiwan

September 4, 2009

(revised on July 21, 2010, December 2, 2010)

**Abstract. We propose a primal-dual continuation approach for the capacitated multi-**
facility Weber problem (CMFWP) based on its nonlinear second-order cone program (SOCP)
reformulation. The main idea of the approach is to reformulate the CMFWP as a nonlinear
SOCP with a nonconvex objective function, and then introduce a logarithmic barrier term
and a quadratic proximal term into the objective to construct a sequence of convexiﬁed
subproblems. By this, this class of nondiﬀerentiable and nonconvex optimization problems
is converted into the solution of a sequence of nonlinear convex SOCPs. In this paper,
we employ the semismooth Newton method proposed in [17] to solve the KKT system of
the resulting convex SOCPs. Preliminary numerical results are reported for eighteen test

1Member of Mathematics Division, National Center for Theoretical Sciences, Taipei Oﬃce. The author’s work is partially supported by National Science Council of Taiwan. E-mail: jschen@math.ntnu.edu.tw, FAX: 886-2-29332342.

2The author’s work is supported by the Fundamental Research Funds for the Central Universities (SCUT), Guangdong Natural Science Foundation (No. 9251802902000001) and National Young Natural Science Foundation (No. 10901058). E-mail: shhpan@scut.edu.cn

3E-mail: chko@isu.edu.tw

instances, which indicate that the continuation approach is promising to ﬁnd a satisfying suboptimal solution, even a global optimal solution for some test problems.

**Key words: Capacitated multi-facility Weber problem, nondiﬀerentiable, nonconvex,**
second-order cone program, semismooth Newton method.

**1** **Introduction**

There are many diﬀerent kinds of facility location problems for which various methods have been proposed; see [10, 15, 21, 20, 22, 23, 26, 32, 35, 36] and references therein. The web-site maintained by EWGLA (Euro Working Group on Location Analysis) also serves as a good resource to look for related literature including survey papers, books, and journals. In this paper, we consider the newer but more diﬃcult capacitated multi-facility Weber problem (CMFWP) which plays an important role in the operation and management science.

The CMFWP, also called the capacitated Euclidean distance location-allocation prob- lem in other contexts, is concerned with locating a set of facilities and allocating their capacity to satisfy the demand of a set of customers with known locations so that the total transportation cost is minimized. Supply centers such as plants and warehouses may constitute the facilities, while retailers and dealers may be considered as customers. The mathematical model of CMFWP can be stated as follows:

min

∑*m*
*i=1*

∑*n*
*j=1*

*c**ij**w**ij**∥x**i**− a**j**∥*

*s.t.*

∑*n*
*j=1*

*w*_{ij}*= s*_{i}*,* *i = 1, 2, . . . , m* (1)

∑*m*
*i=1*

*w*_{ij}*= d*_{j}*,* *j = 1, 2, . . . , n*

*w*_{ij}*≥ 0, i = 1, 2, . . . , m; j = 1, 2, . . . , n*
*x*_{i}*∈ IR*^{2}*, i = 1, 2, . . . , m.*

*In the formulation of CMFWP, m is the number of facilities to be located, n is the number*
*of customers, s**i* *is the capacity of facility i, and d**j**is the demand of customer j. Throughout*
this paper, we without loss of generality assume that the total capacity of all facilities equals
the total demand of all customers, i.e.,

∑*m*
*i=1*

*s** _{i}* =

∑*n*
*j=1*

*d*_{j}*.* (2)

*In addition, the allocations w** _{ij}* are unknown variables denoting the amount to be shipped

*from facility i to customer j with the unit shipment cost per unit distance being c*

*ij*. If all

*w** _{ij}* are known, the CMFWP reduces to the traditional convex multi-facility location prob-
lem for which many eﬃcient algorithms (see [10, 13, 22, 33, 34, 26]) have been proposed,

*whereas if all x*

*are ﬁxed, it reduces to the ordinary transportation problem. For the sake*

_{i}*of notation, in the sequel, we denote by x*

_{i}*= (x*

_{i1}*, x*

*) the unknown coordinates of facility*

_{i2}*i, and by a*

_{j}*= (a*

_{j1}*, a*

_{j2}*) the given coordinates of customer j.*

*We see that the objective function of (1) is nondiﬀerentiable at any points where x*_{i}*= a*_{j}*for some i* *∈ {1, 2, . . . , m} and j ∈ {1, 2, . . . , n}, which precludes a direct application of*
*eﬀective gradient-based methods for ﬁnding its solution. Also, it is nonconvex since w** _{ij}* are
unknown decision variables. Therefore, the CMFWP belongs to a class of nondiﬀerentiable

*nonconvex optimization problems subject to (m + n) linear constraints and mn nonnega-*tivity constraints. In fact, Sherali and Nordai [30] have shown that this class of problems

*is NP-hard, even if all demand points a*

*are located on a straight line.*

_{j}For this class of problems, Cooper in his seminal work [7] ﬁrst proposed an exact solu- tion method by explicit enumeration of all extreme points of the transportation polytope, deﬁned by the ﬁrst three groups of constraints of (1). Selim [29] later presented a biconvex cutting plane procedure in his unpublished dissertation. The exact solution method, simi- lar to Cooper’s complete enumeration, can eﬀectively deal with very small instances only.

Recently, Sherali and Nordai [28] developed a branch-and-bound algorithm which is based on a partitioning of all the allocation space and ﬁnitely converges to a global optimum within a speciﬁed percentage tolerance. Apart from these exact methods, some heuristic methods have also been proposed by the reformulation linearization technique [27] or an approximating mixed integer linear programming formulation [1]. We observe that most of these methods are combinatorial primal ones which are designed by exploiting the structure of the problem itself and do not provide any information about the dual solution.

In this paper, we propose a continuous primal-dual approach by converting the CMFWP into the solution of a sequence of nonlinear convex second-order cone programs (SOCPs).

Speciﬁcally, we ﬁrst reformulate the CMFWP as a nonlinear SOCP with a nonconvex cost function, and then introduce a logarithmic barrier term and a quadratic proximal term into the objective to circumvent the nonconvex diﬃculty. Among others, the strict convexity of logarithmic function and the strong convexity of quadratic proximal term are fully used to convexify the objective function of the resulting SOCP. Such a technique is not new which is often used in the literature; see [4] for example. The SOCP reformulation has recently attracted much attention for engineering and operations research problems. However, to our best of knowledge, they are all formed into linear SOCPs for which some softwares using interior point methods [18, 31] can be applied. In contrast, the nonlinear SOCP reformula- tion has little been used since all the aforementioned softwares are only able to solve linear SOCPs. This paper is concerned with the application of the nonlinear SOCP reformulation in the CMFWP, and its main purpose is to propose an alternative continuous approximate

method to handle the class of diﬃcult problems, instead of introducing a highly specialized method in a competition for the best solution and the fasted computation time.

This paper is organized as follows. In Section 2, we review the general convex SOCP and the semismooth Newton method [17] for solving it. Section 3 presents a detailed process of reformulating (1) as a nonlinear SOCP. Section 4 proposes a primal-dual continuation approach for the CMFWP by solving approximately a sequence of convexiﬁed SOCPs. In Section 5, we report the preliminary numerical results for some test problems from [3, 28]

with the continuation approach, and compare the ﬁnal objective values with those yielded by the global optimization methods [28]. Finally, we conclude this paper.

Throughout this paper, *∥ · ∥ denotes the Euclidean norm, IR*^{n}*denotes the space of n-*
dimensional real column vectors, and IR^{n}^{1}*× · · · × IR*^{n}* ^{m}* is identiﬁed with IR

^{n}^{1}

^{+}

^{···+n}*. Thus,*

^{m}*(x*

_{1}

*, . . . , x*

*)*

_{m}*∈ IR*

^{n}^{1}

*×· · ·×IR*

^{n}*is viewed as a column vector in IR*

^{m}

^{n}^{1}

^{+}

^{···+n}

^{m}*. The notations I*

**and 0 denote an identity matrix and a zero matrix of suitable dimension, respectively, and**

*diag(x*

_{1}

*, . . . , x*

_{n}*) means a diagonal matrix with x*

_{1}

*, . . . , x*

*as the diagonal elements. Given*

_{n}*a ﬁnite number of square matrices Q*

_{1}

*, . . . , Q*

*, we denote the block diagonal matrix with*

_{n}*these matrices as block diagonal by diag(Q*

_{1}

*, . . . , Q*

_{n}*). For a diﬀerentiable function f , we*denote by

*∇f(x) and ∇*

^{2}

*xx*

*f (x) the gradient and the Hessian matrix of f at x, respectively.*

*For a diﬀerentiable mapping G : IR*^{n}*→ IR*^{m}*, we denote by G*^{′}*(x)∈ IR*^{m}* ^{×n}* the Jacobian of

*G at x. LetO be an open set in IR*

^{n}*. If G :O → IR*

*is a locally Lipschitz continuous, then*

^{n}*∂*_{B}*G(x) :=*^{{}*H* *∈ IR*^{m}^{×n}*| ∃{x*^{k}*} ⊆ D**G* *: x*^{k}*→ x, G*^{′}*(x** ^{k}*)

*→ H*

^{}}

*is nonempty and called the B-subdiﬀerential of G at x* *∈ O, where D**G* denotes the set of
*points at which G is diﬀerentiable. We assume that the reader is familiar with the concepts*
of (strongly) semismooth functions, and refer to [24, 25] for details.

**2** **Preliminaries**

The convex SOCP is to minimize a convex function over the intersection of an aﬃne linear manifold with the Cartesian product of second-order cones, which can be described as

minimize *g(x)*

*subject to Ax = b,* *x∈ K,* (3)

*where g : IR*^{n}*→ IR is a twice continuously diﬀerentiable convex function, A is an m × n*
*matrix with full row rank, b* *∈ IR** ^{m}* and

*K is the Cartesian product of second-order cones*(SOCs), also called Lorentz cones. In other words,

*K = K*^{n}^{1} *× K*^{n}^{2} *× · · · × K*^{n}^{q}*,* (4)

*where q, n*_{1}*, . . . , n*_{q}*≥ 1, n*1+*· · · + n**q* *= n, and* *K*^{n}* ^{i}* denotes the SOC in IR

^{n}*deﬁned by*

^{i}*K*

^{n}*:=*

^{i}{

*x*_{i}*= (x*_{i1}*, x*_{i2}*, . . . , x*_{in}* _{i}*)

*∈ IR*

^{n}

^{i}*|*

^{√}

*x*

^{2}

_{i2}*+ . . . + x*

^{2}

_{in}*i* *≤ x**i1*

}

(5)
with *K*^{1} denoting the nonnegative real number set IR_{+}. A special case of (4) corresponds
to the nonnegative orthant cone IR^{n}_{+}*, i.e., q = n and n*1 =*· · · = n**q* *= 1. When g is linear,*
clearly, (3) becomes a linear SOCP which has been investigated in many previous works
and the interested reader is referred to the survey papers [2, 19] and the books [5, 6] for
many important applications and theoretical properties.

The treatment of nonlinear convex SOCPs is much more recent and mainly focuses on
the research of eﬀective solution methods. Notice that *K is a closed convex cone which*
is self-dual in the sense that *K equals its dual cone K** ^{∗}* :=

*{y ∈ IR*

^{n}*| ⟨y, x⟩ ≥ 0 ∀x ∈ K}.*

Thus, it is easy to write the Karush-Kuhn-Tucker (KKT) optimality conditions of (3) as

*∇g(x) − A*^{T}*y− λ = 0*
*Ax− b = 0*

*⟨x, λ⟩ = 0, x ∈ K, λ ∈ K.*

(6)

*These conditions are also suﬃcient for optimality since g is convex. Based on the KKT*
system (6), there have been several methods proposed for solving (3), which include the
smoothing Newton methods [8, 12, 16], the merit function method [9], and the semismooth
Newton method [17]. As mentioned in the introduction, this paper is concerned with the
application of the nonlinear SOCP in the CMFWP.

Since the resulting convex SOCPs in Sec. 4 will be solved with the semismooth Newton
*method in [17], we next review it. Let P** _{K}* : IR

^{n}*→ IR*

*denote the Euclidean projection operator onto the cone*

^{n}*K, i.e., P*

*K*

*(z) := argmin*

_{y}

_{∈K}*{∥z − y∥} for any z ∈ IR*

*. Then, from [12, Prop. 4.1], we have that*

^{n}*x− P*_{K}*(x− λ) = 0 ⇐⇒ x ∈ K, λ ∈ K, ⟨x, λ⟩ = 0.* (7)
Consequently, the solution of the convex SOCP (3) is equivalent to ﬁnding the zeros of

*Φ(ω) := Φ(x, y, λ) :=*

*∇g(x) − A*^{T}*y− λ*
*Ax− b*
*x− P**K**(x− λ)*

*.* (8)

*Since P** _{K}* is strongly semismooth by [16, Prop. 4.5], using [11, Theorem 19] then yields
that the operator Φ is at least semismooth, and furthermore, it is strongly semismooth if

*∇*^{2}*xx**g(x) is locally Lipschitz continuous at any x∈ IR** ^{n}*. The semismooth Newton method in
[17] ﬁnds a zero of Φ by applying the nonsmooth Newton method [24, 25] to the semismooth

*system Φ(ω) = 0. In other words, it generates the iterate sequence*

*{ω*

^{k}*= (x*

^{k}*, y*

^{k}*, λ*

*)*

^{k}*} by*

*ω*^{k+1}*:= ω*^{k}*− W**k*^{−1}*Φ(ω*^{k}*),* (9)

*where W*_{k}*is an arbitrary element from the B-subdiﬀerential ∂*_{B}*Φ(ω** ^{k}*) and has the form of

*W** _{k}*=

*∇*^{2}_{xx}*g(x** ^{k}*)

*−A*

^{T}*−I*

*A* **0** **0**

*I− V*^{k}**0** *V*^{k}

*for a suitable block diagonal matrix V*^{k}*= diag(V*_{1}^{k}*, . . . , V*_{q}^{k}*) with V*_{i}^{k}*∈ ∂**B**P*_{K}*ni**(x*^{k}_{i}*− λ*^{k}*i*).

The following two technical lemmas respectively provide the formula to compute the
*value of P*_{K}*at any point and the representation of each element in V* *∈ ∂**B**P*_{K}^{n}*(z).*

**Lemma 2.1 [17, Lemma 2.2] For any given z = (z**_{1}*, z*_{2})*∈ IR × IR*^{n}^{−1}*, it holds that*
*P*_{K}^{n}*(z) = max{0, µ*1*(z)}u*^{(1)}*z* + max*{0, µ*2*(z)}u*^{(2)}*z* *,*

*where µ*_{1}*(z), µ*_{2}*(z) and u*^{(1)}_{z}*, u*^{(2)}_{z}*are the spectral values and the spectral vectors of z, respec-*
*tively, given by*

*µ*_{1}*(z) = z*_{1}*− ∥z*2*∥,* *µ*_{2}*(z) = z*_{1}+*∥z*2*∥;*

*u*^{(1)}* _{z}* = 1
2

(

*1,−¯z*2

)

*,* *u*^{(2)}* _{z}* = 1
2

(

*1, ¯z*_{2}^{)}
*with ¯z*_{2} = _{∥z}^{z}^{2}

2*∥* *if z*_{2} *̸= 0 and otherwise ¯z*2 *being any vector in IR*^{n−1}*satisfying* *∥¯z*2*∥ = 1.*

* Lemma 2.2 [17, Lemma 2.6] Given a general point z = (z*1

*, z*2)

*∈ IR×IR*

^{n}

^{−1}*, each element*

*V*

*∈ ∂*

*B*

*P*

_{K}

^{n}*(z) has the following representation:*

**(a) If z**_{1} *̸= ±∥z*2*∥, then P*_{K}^{n}*(z) is continuously diﬀerentiable with*

*V = P*_{K}^{′}*n**(z) =*

**0** *if z*_{1} *<−∥z*2*∥*
*I* *if z*1 *>∥z*2*∥*
1

2

( 1 *z*¯_{2}^{T}

¯
*z*_{2} *H*

)

if *− ∥z*2*∥ < z*1 *<∥z*2*∥,*
*where*

¯

*z*_{2} := *z*_{2}

*∥z*2*∥, H :=*

(

1 + *z*_{1}

*∥z*2*∥*

)

*I−* *z*_{1}

*∥z*2*∥z*¯_{2}*z*¯_{2}^{T}*.*
**(b) If z**_{2} *̸= 0 and z*1 =*∥z*2*∥, then*

*V* *∈*

{

*I,* 1
2

( 1 *z*¯_{2}^{T}

¯
*z*_{2} *H*

)}

*, where ¯z*_{2} := *z*_{2}

*∥z*2*∥* *and H := 2I− ¯z*2*z*¯_{2}^{T}*.*

**(c) If z**_{2} *̸= 0 and z*1 =*−∥z*2*∥, then*
*V* *∈*

{

* 0,* 1
2

( 1 *z*¯^{T}_{2}

¯
*z*2 *H*

)}

*, where ¯z*_{2} := *z*_{2}

*∥z*2*∥* *and H := ¯z*_{2}*z*¯_{2}^{T}*.*
**(d) If z = 0, then either V = 0 or V = I or V belongs to the set**

{1 2

( 1 *z*¯_{2}^{T}

¯
*z*_{2} *H*

) *H = (w*_{0}*+ 1)I− w*0*z*¯_{2}*z*¯_{2}* ^{T}* for some

*|w*0

*| ≤ 1 and ∥¯z*2

*∥ = 1*

}

*.*

**3** **Nonlinear SOCP reformulation**

In this section, we will present the detailed process of reformulating the capacitated multi-
facility Weber problem (1) as a nonlinear SOCP with the form of (3). First, by introducing
*mn new variables t**ij* *for i = 1, 2, . . . , m, j = 1, 2, . . . , n, we can transform (1) into*

min

∑*m*
*i=1*

∑*n*
*j=1*

*c*_{ij}*w*_{ij}*t*_{ij}

*s.t.* *∥x**i**− a**j**∥ ≤ t**ij**, i = 1, 2, . . . , m; j = 1, 2, . . . , n*

∑*n*
*j=1*

*w*_{ij}*= s*_{i}*,* *i = 1, 2, . . . , m* (10)

∑*m*
*i=1*

*w*_{ij}*= d*_{j}*,* *j = 1, 2, . . . , n*

*w*_{ij}*≥ 0, i = 1, 2, . . . , m; j = 1, 2, . . . , n*
*x*_{i}*∈ IR*^{2}*, i = 1, 2, . . . , m.*

We write out all the constraints of *∥x**i**− a**j**∥ ≤ t**ij**, i = 1, 2, . . . , m; j = 1, 2, . . . , n as below:*

√

*(x*_{i1}*− a*11)^{2}*+ (x*_{i2}*− a*12)^{2} *≤ t**i1**, i = 1, 2, . . . , m,*

√

*(x**i1**− a*21)^{2}*+ (x**i2**− a*22)^{2} *≤ t**i2**, i = 1, 2, . . . , m,* (11)
... ... ... ...

√

*(x*_{i1}*− a**n1*)^{2}*+ (x*_{i2}*− a**n2*)^{2} *≤ t**in**, i = 1, 2, . . . , m.*

Let _{{}

*u*_{ij}*:= x*_{i1}*− a**j1**, i = 1, 2, . . . , m,* *j = 1, 2, . . . , n;*

*v**ij* *:= x**i2**− a**j2**, i = 1, 2, . . . , m,* *j = 1, 2, . . . , n.* (12)
Then, the constraints in (11) turn into

*(u** _{i1}*)

^{2}

*+ (v*

*)*

_{i1}^{2}

*≤ (t*

*i1*)

^{2}

*, t*

_{i1}*≥ 0, i = 1, 2, . . . , m,*

*(u** _{i2}*)

^{2}

*+ (v*

*)*

_{i2}^{2}

*≤ (t*

*i2*)

^{2}

*, t*

_{i2}*≥ 0, i = 1, 2, . . . , m,*(13)

... ... ... ... ...

*(u**in*)^{2}*+ (v**in*)^{2} *≤ (t**in*)^{2}*, t**in* *≥ 0, i = 1, 2, . . . , m.*

Let

ˆ
*x** _{ij}* :=

*t**ij*

*u*_{ij}*v*_{ij}

*∈ IR*^{3}*,* *i = 1, 2, . . . , m, j = 1, 2, . . . , n.* (14)

From (13) and the deﬁnition of*K*^{3}, we have ˆ*x*_{ij}*∈ K*^{3} *for i = 1, 2, . . . , m and j = 1, 2, . . . , n.*

It should be pointed out that we have created additional constraints through the above reformulation procedure. In particular, from (12), it follows that

*u*_{11}*− u**1k* *= a*_{k1}*− a*11*,* *k = 2, 3, . . . , n,*
... ...

*u*_{m1}*− u**mk* *= a*_{k1}*− a*11*,* *k = 2, 3, . . . , n* (15)
and

*v*_{11}*− v**1k* *= a*_{k2}*− a*12*,* *k = 2, 3, . . . , n,*
... ...

*v*_{m1}*− v**mk* *= a*_{k2}*− a*12*,* *k = 2, 3, . . . , n.* (16)
We will see that (15)–(16) can be recast as a linear system. Note that (15) and (16) are
equivalent to

[0 1 0]

*t*_{11}
*u*_{11}
*v*_{11}

+ [0 *−1 0]*

*t*_{1k}*u*_{1k}*v*_{1k}

*= a*_{k1}*− a*11*,* *k = 2, 3, . . . , n,*

... ... ...

[0 1 0]

*t*_{m1}*u*_{m1}*v**m1*

+ [0 *−1 0]*

*t*_{mk}*u*_{mk}*v**mk*

*= a*_{k1}*− a*11*, k = 2, 3, . . . , n*

and

[0 0 1]

*t*_{11}
*u*11

*v*_{11}

+ [0 0 *−1]*

*t*_{1k}*u**1k*

*v*_{1k}

*= a*_{k2}*− a*12*,* *k = 2, 3, . . . , n*

... ... ...

[0 0 1]

*t*_{m1}*u*_{m1}*v*_{m1}

+ [0 0 *−1]*

*t*_{mk}*u*_{mk}*v*_{mk}

*= a*_{k2}*− a*12*, k = 2, 3, . . . , n.*

We can simplify (15)–(16) by introducing the following notations. More speciﬁcally, let

*A** _{u}* :=

0 1 0 0 *−1 0*

0 1 0 0 *−1 0*

... . ..

0 1 0 0 *−1 0*

*∈ IR*^{(n−1)×3n}*,* (17)

*A** _{l}* :=

0 0 1 0 0 *−1*

0 0 1 0 0 *−1*

... . ..

0 0 1 0 0 *−1*

*∈ IR*^{(n}^{−1)×3n}*,* (18)

and denote

*b** _{u}* :=

*a*21*− a*11

*a*_{31}*− a*11

...
*a**n1**− a*11

*∈ IR*^{n}^{−1}*,* *b** _{l}* :=

*a*22*− a*12

*a*_{32}*− a*12

...
*a**n2**− a*12

*∈ IR*^{n}^{−1}*.*

Then, equations (15)–(16) can be recast as the following system of linear constraints

*A**u*

*A*_{u}

. ..

*A**u*

*−− − − − − − − − − −*
*A*_{l}

*A*_{l}

. ..

*A*_{l}

[ˆ*x*_{11}]
...
[ˆ*x**1n*]

*−−*

[ˆ*x*_{21}]
...
[ˆ*x** _{2n}*]

*−−*...

*−−*

[ˆ*x** _{m1}*]
...
[ˆ

*x*

*]*

_{mn}

=

*b**u*

*b** _{u}*
...

*b*

*u*

*−−*

*b*_{l}*b** _{l}*
...

*b*

_{l}

*,* (19)

*where the dimensions in the linear system are 2m(n− 1) × 3mn, 3mn × 1, 2m(n − 1) × 1*
for the matrix, column of variables, and column of constants, respectively.

We next look into the constraints on demand and capacity, namely the two groups of

equality constraints in (10). We notice that they can be recast as a linear system as below:

[1 1 1 *· · · 1]*

[1 1 1 *· · · 1]*

. ..

[1 1 1 *· · · 1]*

*− − − − − − − − − − − − − − − − − − − −*
[1 0 0*· · · 0] [1 0 0 · · · 0]* *· · ·* [1 0 0*· · · 0]*

[0 1 0*· · · 0] [0 1 0 · · · 0]* *· · ·* [0 1 0*· · · 0]*

... ... ... ...

[0 0 0*· · · 1] [0 0 0 · · · 1]* *· · ·* [0 0 0*· · · 1]*

*w*_{11}
...
*w*_{1n}

*−−*

*w*_{21}
...
*w**2n*

*−−*...

*−−*

*w** _{m1}*
...

*w*

*mn*

=

*s*_{1}
*s*_{2}
...
*s*_{m}

*−−*

*d*_{1}
*d*_{2}
...
*d*_{n}

*.* (20)

*Again, we point it out that the dimensions of the system (20) are (m + n)×(mn), (mn)×1,*
*(m + n)× 1 for the matrix, column of variables and column of constants, respectively.*

Moreover, the coeﬃcient matrix has not full row rank due to the assumption (2), but its
*any m + n− 1 rows are all linear independent. For convenience, in the rest of this paper,*
*A*_{w}*denotes the matrix composed of the ﬁrst m + n−1 rows in the coeﬃcient matrix of (20).*

In summary, we reformulate the CMFWP (1) as the following nonlinear SOCP:

minimize

∑*m*
*i=1*

∑*n*
*j=1*

*c*_{ij}*w*_{ij}*t*_{ij}*subject to Ax = b*

*x∈ (K*^{3})^{mn}*× (K*^{1})^{mn}

(21)

where (*K*^{3})^{mn}*× (K*^{1})^{mn}*denotes the Cartesian product of mn* *K*^{3} *and mnK*^{1}, and

*x := (ˆx*_{11}*, . . . , ˆx*_{1n}*, . . . , ˆx*_{m1}*, . . . ˆx*_{mn}*, w*_{11}*, . . . , w*_{1n}*, . . . , w*_{m1}*, . . . , w*_{mn}*) ,* (22)

*A :=*

*A*_{u}*|*

*A**u* *|*

. .. *|*

*A*_{u}*|*

*− − − − − − − − − − − − |* **0**

*A*_{l}*|*

*A*_{l}*|*

. .. *|*

*A*_{l}*|*

*− − − − − − − − − − − − | − − −*

**0** *|* *A*_{w}

*,* (23)

*b := (b*_{u}*, b*_{u}*, . . . , b*_{u}*, b*_{l}*, b*_{l}*, . . . , b*_{l}*, s*_{1}*, . . . , s*_{m}*, d*_{1}*, . . . , d*_{n}_{−1}*) .* (24)
*Notice that, by the expression (22) of x, the objective function of (21) can be rewritten as*
*a quadratic function x*^{T}*Qx, where Q = [Q** _{kl}*]

_{4mn}*is a nonsymmetric matrix with*

_{×4mn}*Q**kl*=

{ *c*_{ij}*if k = 3(i− 1)n + 3(j − 1) + 1, l = 3mn + (i − 1)n + j;*

0 otherwise (25)

*for k, l = 1, 2, . . . , 4mn and i = 1, 2, . . . , m; j = 1, 2, . . . , n. Hence, (21) is equivalent to*
minimize *x*^{T}*Qx*

*subject to Ax = b*

*x∈ (K*^{3})^{mn}*× (K*^{1})^{mn}

(26)

*where Q is a 4mn× 4mn matrix given by (25), A is a (2mn − m + n − 1) × 4mn matrix*
*with full row rank (2mn− m + n − 1), and the dimension of b is (2mn − m + n − 1) × 1.*

*Since x*^{T}*Qx = x*^{T}*Q*^{T}*x for any x∈ IR** ^{4mn}*, the SOCP (21) is further equivalent to
minimize

*f (x) :=*1

2*x*^{T}*(Q + Q*^{T}*)x*
*subject to Ax = b*

*x∈ (K*^{3})^{mn}*× (K*^{1})^{mn}

(27)

The SOCP (27) has an advantage over (26) that its Hessian matrix is symmetric, although its objective function is still nonconvex. In the next section, we develop a primal-dual algorithm for solving the CMFWP based on the nonlinear SOCP reformulation (27).

*We should point out that the Hessian matrix of the function f , i.e., ¯Q = (Q + Q** ^{T}*) and

*the constraint coeﬃcient matrix A of (27) are both sparse. For example, when m = n = 2,*

*Q = [ ¯*¯

*Q*

*]*

_{kl}_{16}

*has only eight nonzero entries:*

_{×16}*Q*¯

*= ¯*

_{1,13}*Q*

_{13,1}*= c*

_{11}

*, ¯Q*

*= ¯*

_{4,14}*Q*

*=*

_{14,4}*c*

_{12}

*, ¯Q*

*= ¯*

_{7,15}*Q*

_{15,7}*= c*

_{21}

*, ¯Q*

*= ¯*

_{10,16}*Q*

_{16,10}*= c*

_{22}

*, whereas A = [A*

*]*

_{ij}_{7}

*has sixteen nonzero*

_{×16}*entries: A*12

*= 1, A*15 =

*−1, A*28

*= 1, A*

*2,11*=

*−1, A*33

*= 1, A*36 =

*−1, A*49

*= 1, A*

*4,12*=

*−1, A**5,13* *= A*_{5,14}*= 1, A*_{5,15}*= A*_{5,16}*= 1, A*_{6,13}*= A*_{6,15}*= 1, A*_{7,14}*= A** _{7,16}* = 1. In
addition, we observe that each row of ¯

*Q has at most a nonzero entry c*

_{i}_{1}

_{,j}_{1}

*with c*

_{i}_{1}

_{,j}_{1}

*∈*

*{c*11

*, . . . , c*

_{1n}*, . . . , c*

_{m1}*, . . . , c*

_{mn}*}, and all diagonal entries are zero.*

**4** **Continuation approach for the CMFWP**

This section designs a primal-dual approximate algorithm for the nonconvex SOCP (27).

The main idea is transforming (27) into the solution of a sequence of convexiﬁed SOCP subproblems.

*First, note that the variables w** _{ij}* are nonnegative and restricted by the linear constraints

∑_{n}

*j=1**w*_{ij}*= s** _{i}* and

^{∑}

^{m}

_{i=1}*w*

_{ij}*= d*

*, and hence the nonconvex SOCP (27) is equivalent to*

_{j}*min f (x)*

s.t. *Ax = b*

0*≤ w**ij* *≤ min{s**i**, d**j**}, i = 1, . . . , m; j = 1, . . . , n*
*x∈ (K*^{3})^{mn}*× (K*^{1})^{mn}*.*

(28)

Since the logarithmic barrier function *−*^{∑}^{m}_{i=1}^{∑}^{n}_{j=1}*[ln(w** _{ij}*) + ln(min

*{s*

*i*

*, d*

_{j}*} − w*

*ij*)] is well

*deﬁned when 0 < w*

_{ij}*< min{s*

*i*

*, d*

_{j}*} for all i = 1, 2, . . . , m and j = 1, 2, . . . , n, and moreover,*its limit is +

*∞ if some w*

*ij*tends to 0 or min

*{s*

*i*

*, d*

*j*

*}, we can dispense with the bound*

*constraints on w*

*in (28) and obtain the following transformed problem:*

_{ij}*min f (x)− τ*^{∑}^{m}

*i=1*

∑*n*
*j=1*

[

*ln(w** _{ij}*) + ln(min

*{s*

*i*

*, d*

_{j}*} − w*

*ij*)

^{]}s.t.

*Ax = b*

*x∈ (K*^{3})^{mn}*× (K*^{1})^{mn}

(29)

*where τ > 0 is a barrier parameter that ﬁnally tends to 0. The above operation seems to*
purposely make the original problem (27) more complex. However, we will see that the
introduction of the logarithmic barrier term gives a help in locating the feasible solution of
*(27), as well as plays a certain role in convexifying f (x).*

Now, we introduce a quadratic proximal term 1

2*∥x − z∥*^{2} *with z* *∈ IR** ^{4mn}* being a given

*vector into the objective of (29) to further convexify f (x). Deﬁne*

*g(x, z, τ, ε) := f (x)− τ*^{∑}^{m}

*i=1*

∑*n*
*j=1*

[

*ln(w** _{ij}*) + ln(min

*{s*

*i*

*, d*

_{j}*} − w*

*ij*)

^{]}+

*ε*

2*∥x − z∥*^{2} (30)
*where ε > 0 is a proximal parameter that will become larger until over some threshold.*

**Proposition 4.1 Given a vector z***∈ IR*^{4mn}*, let g be deﬁned as in (30) for any τ, ε > 0.*

*Then, there exists ε*_{0} *> 0 such that g is strictly convex for any ε > ε*_{0} *and τ > 0.*

**Proof. For given z**∈ IR^{4mn}*and any τ, ε > 0, we compute the Hessian matrix of g as*

*∇*^{2}*xx**g(x, z, τ, ε) = Q + Q*^{T}*+ diag(ε, . . . , ε*

| {z }

*3mn*

*, P*11*+ε, . . . , P**1n**+ε, . . . , P**m1**+ε, . . . , P**mn**+ε)*

*with P*_{ij}*for i = 1, 2, . . . , m and j = 1, 2, . . . , n given by*
*P** _{ij}* =

*τ*

*(w** _{ij}*)

^{2}+

*τ*

(min*{s**i**, d*_{j}*} − w**ij*)^{2} *> 0.*

*Since all diagonal entries of Q + Q*^{T}*are zero and each row of Q + Q** ^{T}* has at most a nonzero

*entry c*

_{i}_{1}

_{,j}_{1}

*with c*

_{i}_{1}

_{,j}_{1}

*∈ {c*11

*, . . . , c*

_{1n}*, c*

_{21}

*, . . . , c*

_{2n}*, . . . , c*

_{m1}*, . . . , c*

_{mn}*}, we have that*

[*∇*^{2}*xx**g(x, z, τ, ε)*^{]}

*kk* *>*

*4mn*∑

*l=1,l̸=k*

[*∇*^{2}*xx**g(x, z, τ, ε)*^{]}

*kl*

*for all k = 1, 2, . . . , 4mn whenever ε > ε*_{0} with
*ε*_{0} := max

1*≤i≤m,1≤j≤n**{c**ij**}.* (31)

In other words, the matrix *∇*^{2}*xx**g(x, z, τ, ε) is strictly diagonally dominant for any ε > ε*_{0}
*and τ > 0. From Corollary 7.2.3 of [14], it then follows that* *∇*^{2}_{xx}*g(x, z, τ, ε) is positive*
*deﬁnite, and consequently g(x, z, τ, ε) is strictly convex, for any ε > ε*_{0} *and τ > 0.* *2*

*Proposition 4.1 states that for a given z* *∈ IR*^{4mn}*, the function g(x, z, τ, ε) is strictly*
*convex for any ε > ε*_{0} *and τ > 0. In fact, it is also strongly convex when ε > ε*_{0}. In view of
this, our continuation approach seeks for an approximate optimal solution of the CMFWP
by solving a sequence of the following subproblems

*min g(x, ˆx*^{k}*, τ*_{k}*, ε** _{k}*)
s.t.

*Ax = b*

*x∈ (K*^{3})^{mn}*× (K*^{1})^{mn}

(32)

with a decreasing sequence*{τ**k**} and a increasing sequence {ε**k**}, where ˆx** ^{k}* is a vector given

*by the solution of the last subproblem. For a ﬁxed k, since the subproblem (32) is a convex*

*SOCP whenever ε*

_{k}*> ε*

_{0}, its solution is easy, which from Section 2 is equivalent to solving

*∇**x**g(x, ˆx*^{k}*, τ*_{k}*, ε** _{k}*)

*− A*

^{T}*y− λ = 0*

*Ax− b = 0*

*x− P**K**(x− λ) = 0*

(33)

with *K = (K*^{3})^{mn}*× (K*^{1})* ^{mn}*. Deﬁne the mapping Φ

*: IR*

_{k}

^{10mn}

^{−m+n−1}*→ IR*

^{10mn}*by*

^{−m+n−1}Φ_{k}*(ω) = Φ*_{k}*(x, y, λ) :=*

*∇**x**g(x, ˆx*^{k}*, τ*_{k}*, ε** _{k}*)

*− A*

^{T}*y− λ*

*Ax− b*

*x− P**K**(x− λ)*

*.* (34)

Then, solving the nonsmooth system (33) is equivalent to ﬁnding the zero of the operator
Φ_{k}*(ω). We will attain this goal by using the semismooth Newton method in [17].*

Next we describe the iteration steps of the continuation approach in which two ﬁxed
*constants c*_{1} *∈ (0, 1) and c*2 *> 1 are used to reduce and increase the dynamic parameters τ*
*and ε. Let Ψ*_{k}*(w) :=* 1

2*∥Φ**k**(ω)∥*^{2} denote the natural merit function of system Φ_{k}*(ω) = 0.*

**Algorithm 4.1 (Continuation Approach )**

**(S.0) Given the constants c**_{1} *∈ (0, 1), c*2 *> 1 and ˆτ , ˆε > 0. Select τ*_{0}*, ε*_{0} *> 0 and a starting*
*point ¯ω*^{0} = (¯*x*^{0}*, ¯y*^{0}*, ¯λ*^{0}*) with the last mn elements ¯w*^{0}_{ij}*of ¯x*^{0} *satisfying 0 < ¯w*^{0}_{ij}*<*

*min{s**i**, d*_{j}*}. Let ˆx*^{0} := ¯*x*^{0}*, and set k := 0.*

**(S.1) Compute (by Algorithm 4.2 below) an approximate optimal solution ω**^{k}*= (x*^{k}*, y*^{k}*, λ** ^{k}*)

*of (33) with the starting point ¯ω*

*= (¯*

^{k}*x*

^{k}*, ¯y*

^{k}*, ¯λ*

^{k}*).*

**(S.2) If τ**_{k}*< ˆτ and ε*_{k}*> ˆε, then stop, and otherwise go to (S.3).*

**(S.3) Modify the parameters τ**_{k}*and ε*_{k}*by τ*_{k+1}*:= c*_{1}*τ*_{k}*and ε*_{k+1}*:= c*_{2}*ε*_{k}*, and let*

¯

*x*^{k+1}*:= x*^{k}*,* *y*¯^{k+1}*:= y*^{k}*,* *λ*¯^{k+1}*:= z*^{k}*,* *x*ˆ^{k+1}*:= x*^{k}*.*
**(S.4) Set k := k + 1, and go to (S.1).**

*Note that the continuation approach is a primal-dual one and the last mn components*
*of the ﬁnal iterate y** ^{k}* provide an approximate dual solution of the CMFWP, which usually
has a certain economic meaning associated with this class of transportation problems. In
addition, the main computation work of Algorithm 4.1 is to seek an approximate solution

*of (33). Such a solution can be easily obtained when ε*

*is large enough since, on the one hand, the mapping*

_{k}*∇*

*x*

*g(x, ˆx*

^{k}*, τ*

_{k}*, ε*

*) is strongly monotone in this case, which together with Lemma 3.1 of [16] implies that the function Ψ*

_{k}

_{k}*(ω) has bounded level sets; and on the other*hand, from Section 2, the operator Φ

*is at least semismooth which guarantees in theory a fast algorithm with a superlinear (or quadratic) convergence rate, to seek the solution of (33). Of course, we should emphasize that, when computing the approximate solution of (33), we must restrict the maximum steplength of the iterates (see Algorithm 4.2 below).*

_{k}We next describe the speciﬁc iteration steps of the semismooth Newton method [17]

*when applying it for solving the semismooth system (33). For a given k, from Section 2 it*
follows that the main iteration step of the method is as follows:

*ω*^{l+1}*:= ω*^{l}*− W**l** ^{−1}*Φ

_{k}*(ω*

^{l}*), W*

_{l}*∈ ∂*

*B*Φ

_{k}*(ω*

^{l}*) for l = 0, 1, 2, . . . ,*

*where ∂**B*Φ*k**(ω** ^{l}*) denotes the B-subdiﬀerential of the semismooth mapping Φ

*k*at the point

*ω*

^{l}*and any element W*

_{l}*of ∂*

*Φ*

_{B}

_{k}*(ω*

*) has the expression*

^{l}*W** _{l}* :=

*∇*^{2}_{xx}*g(x*^{l}*, ˆx*^{k}*, ε*_{k}*, τ** _{k}*)

*−A*

^{T}*I*

*A* **0** **0**

*I− V*^{l}**0** *V*^{l}

(35)

*for a suitable block diagonal matrix V*^{l}*= diag(V*_{1}^{l}*, . . . , V*_{4mn}^{l}*) with V*_{i}^{l}*∈ ∂**B**P*_{K}^{3}*(x*^{l}*− λ** ^{l}*) for

*i = 1, 2, . . . , 3mn and V*

_{i}

^{l}*∈ ∂*

*B*

*P*

_{K}^{1}

*(x*

^{l}*− λ*

^{l}*) for i = 3mn + 1, . . . , 4mn.*

**Algorithm 4.2 (Solving the subproblem (33) for Step (S.1) of Algorithm 4.1)**