---
title: "Group Sparse Optimal Transport for Sparse Process Flexibility Design"
authors:
  - "Dixin Luo"
  - "Tingting Yu"
  - "Hongteng Xu"
doi: "10.24963/ijcai.2023/679"
url: "https://ijcai.org/proceedings/2023/679"
pdf: "https://ijcai.org/proceedings/2023/0679.pdf"
venue: "Proceedings of the Thirty-Second International Joint Conference on Artificial Intelligence"
year: 2023
pages: "6121-6129"
track: "AI for Good"
keywords:
  - "AI for Good: Multidisciplinary Topics and Applications"
---

# Group Sparse Optimal Transport for Sparse Process Flexibility Design

Dixin Luo, Tingting Yu, Hongteng Xu

*Proceedings of the Thirty-Second International Joint Conference on Artificial Intelligence*, pages 6121-6129. DOI: [10.24963/ijcai.2023/679](https://doi.org/10.24963/ijcai.2023/679)

## Abstract

As a fundamental problem in Operations Research, sparse process flexibility design (SPFD) aims to design a manufacturing network across industries that achieves a trade-off between the efficiency and robustness of supply chains. 
In this study, we propose a novel solution to this problem with the help of computational optimal transport techniques.
Given a set of supply-demand pairs, we formulate the SPFD task approximately as a group sparse optimal transport (GSOT) problem, in which a group of couplings between the supplies and demands is optimized with a group sparse regularizer. 
We solve this optimization problem via an algorithmic framework of alternating direction method of multipliers (ADMM), in which the target network topology is updated by soft-thresholding shrinkage, and the couplings of the OT problems are updated via a smooth OT algorithm in parallel. 
This optimization algorithm has guaranteed convergence and provides a generalized framework for the SPFD task, which is applicable regardless of whether the supplies and demands are balanced. 
Experiments show that our GSOT-based method can outperform representative heuristic methods in various SPFD tasks.
Additionally, when implementing the GSOT method, the proposed ADMM-based optimization algorithm is comparable or superior to the commercial software Gurobi. 
The code is available at https://github.com/Dixin-s-Lab/GSOT.

## Full text

## 1 Introduction

Sparse process flexibility design (SPFD) is a fundamental Operations Research (OR) problem that aims to design a sparse bipartite manufacturing network to make plants shift production portfolios according to the dynamics of demands with little penalty in terms of time and cost. This problem originates in automotive manufacturing [Jordan and Graves, 1995] and commonly appears in various industrial scenarios, e.g., supply chain design [Graves and Tomlin, 2003; SimchiLevi et al., 2018], medical resource management [Chan et al., 2022], express delivery system optimization [Chou et al., 2008], and so on. Because of its significance in industry, the study of the SPFD problem has both theoretical values and social impacts, which is highly correlated with the 8th and 9th United Nations Sustainable Development Goals, i.e., promoting sustained, inclusive, and sustainable economic growth, building resilient infrastructure, promoting inclusive and sustainable industrialization, and fostering innovation. Especially for those developing countries, the flexibility of their supply chains is critical for making their economy robust to the uncertainty in the global market. Solving the SPFD problem helps optimize their existing supply chains and provides valuable guidance for establishing their future industrial infrastructures.

∗Correspondence author.

In general, the SPFD problem corresponds to a multistage optimization problem. In particular, we need to estimate the significance of edges in the current manufacturing network and update the network topology via adding or deleting edges based on the edge significance [Chou et al., 2008]. Repeating the two steps above leads to a heuristic search process. Most existing methods implement the above heuristic strategy in different ways. When estimating the edge significance, they often solve the maximum flow (max-flow) problem [Harris and Ross, 1956] or learn a neural network-based surrogate of the max-flow solver in the reinforcement learning framework [Wei et al., 2021; Chan et al., 2022]. Given the edge significance, the network topology is updated in a deterministic [Chou et al., 2010] or stochastic [Chen et al., 2015] way. These methods have been applied in many scenarios, whose performance even has theoretical guarantees in some settings. However, the methods and their theory are often developed under questionable assumptions, e.g., the same cost/profit per edge, same-sized supply and demand sets, specified topology type (like long chain or k-chain), and so on, which limits their practical applications. Moreover, none of the methods consider jointly optimizing the edge significance and the network topology in a scalable framework. As a result, they often lead to local optimum for large-scale applications [Feng et al., 2017].

To overcome the drawbacks mentioned above, we propose a novel approximate solution, called group sparse optimal transport (GSOT), to the SPFD problem based on computational optimal transport techniques [Peyr´e et al., 2019]. As illustrated in Figure 1, given predefined supplies and a set of demands sampled from some distributions (as existing work [Feng et al., 2017; Wei et al., 2021; Chan et al., 2022] did), we jointly approximate the SPFD problem by solving a set of optimal transport (OT) problems. Each OT problem corresponds to maximizing the expected profit of a manufacturing network given a supply-demand pair. All OT problems are optimized jointly with a group sparse regularizer [Bengio et al., 2009; Deng et al., 2013], leading to a set of couplings (or called optimal transport plans) with a consistent sparse structure. The aggregation of the group sparse couplings indicates the flexible network topology that is robust to the demand uncertainty.

Figure 1: An illustration of the proposed method.

![Image](https://ijcai.org/proceedings/2023/fulltext/679/image_000000_cd2401f11e9d1ca424399bdd628dd336fc706c15da8f9d4e5d72c4dc15027bd7.png)

Instead of applying heuristic multi-stage optimization strategies, our GSOT approximates the SPFD problem by solving a convex optimization problem, whose convergence to the global optimum is theoretically-guaranteed. In particular, we design an algorithmic framework based on the alternating direction method of multipliers (ADMM) [Boyd et al., 2011], solving the optimization problem efficiently. The ADMM framework updates the target network topology iteratively, in which a group-lasso operation [Simon et al., 2013] ensures the group sparsity of the couplings. Each coupling is updated by solving a quadratically-regularized optimal transport problem, leading to the smooth OT algorithm [Blondel et al., 2018]. Our GSOT method is applicable for both balanced and unbalanced supply-demand settings, and the ADMM framework is feasible for large-scale SPFD problems. The comparisons with existing heuristic methods on synthetic and real-world datasets demonstrate the effectiveness of our method.

## 2 RelatedWork

## 2.1 Sparse Process Flexibility Design

Focusing on the SPFD problem, many solutions have been proposed. The early work in [Jordan and Graves, 1995] proposed a concept called "long chain", which designed a flexible process empirically by constructing a long circular chain visiting as many nodes as possible. Taking the long chain as the design principle, many heuristic methods have been proposed [Graves and Tomlin, 2003]. These methods design a sparse manufacturing network by adding edges to a null graph [Chou et al., 2011; Feng et al., 2017] or deleting edges from a complete bipartite graph [Yan et al., 2018; Simchi-Levi et al., 2018] step by step.

In each step, different methods select one or multiple candidate edges to add or delete by various criteria, which can be coarsely categorized into three classes. The optimizationbased methods select the edges by solving a max-flow problem [Harris and Ross, 1956] or its variants (e.g., its dual problem [Yan et al., 2018] and its stochastic version [Feng et al., 2017]) given the current network topology. The node expansion methods adapt the concept of graph expanders in graph theory, adding edges iteratively to improve upon the node expansion ratio in a greedy manner [Chou et al., 2008]. Theoretically, the "probabilistic graph expander" introduced in [Chen et al., 2015] provides an upper bound for the number of edges that guarantees (1-ϵ)-optimality with a high probability in balanced and symmetric systems. The work in [Chen et al., 2019] further generalizes the theoretical result to unbalanced and asymmetric systems with O(n ln(1/ϵ)) edges. Recently, some learning-based SPFD methods have been developed to solve combinatorial problems [Khalil et al., 2017; Sultana et al., 2022; Barrett et al., 2020]. These methods formulate the steps of adding edges and their influences on the cost of the whole supply chain as an action-state sequence in the reinforcement learning (RL) framework. Instead of solving the max-flow problem per step, they train a policy model to imitate the optimizer [Wei et al., 2021; Chan et al., 2022]. The RL-based methods are more efficient than the optimization-based methods in the deploying phase because of avoiding iterative optimization. However, they often have poor scalability and limited generalization power.

Although the above methods achieve encouraging performance in some cases, their nature of heuristic search often leads to undesired local optimum in practical applications. In theory, the optimality of the long chain is only held in some ideal situations, as discussed in [D´esir et al., 2016]. As aforementioned, most existing methods ignore the grounding cost/profit associated with each edge. Moreover, none consider reformulating the two-stage strategy (i.e., the network update and evaluation) in a joint optimization framework.

## 2.2 Computational Optimal Transport

Optimal transport (OT) theory [Villani, 2008] has proven to be useful in machine learning tasks, e.g., distribution matching [Frogner et al., 2015], graph matching [Maretic et al., 2022], data clustering [Agueh and Carlier, 2011; Cuturi and Doucet, 2014], information fusion [Xu and Cheng, 2022], and generative modeling [Arjovsky et al., 2017]. Given the samples of two distributions, the discrete OT problem corresponds to a linear programming problem [Kusner et al., 2015], which optimizes a joint distribution of the samples (a.k.a., the coupling or the optimal transport plan) to minimize the expected pairwise distance between the two distributions' samples. Typically, this problem is solved by the conditional gradient algorithm [Titouan et al., 2019] or the Bregman ADMM algorithm [Wang and Banerjee, 2014; Ye et al., 2017]. To accelerate the optimization step, the work in [Cuturi, 2013] proposed a surrogate called Sinkhorn distance, which makes the OT problem strictly convex by adding an entropic regularizer to the coupling. The entropic OT problem can be solved efficiently by the Sinkhorn-scaling algorithm [Sinkhorn and Knopp, 1967; Benamou et al.,

2015]. To further improve the numerical stability and scalability of the Sinkhorn-scaling algorithm, some variants have been proposed, including the proximal point method [Xie et al., 2020] and the stochastic Sinkhorn algorithm (called Greenkhorn) [Altschuler et al., 2017]. Besides the entropic OT problem, the quadratically-regularized OT problem is also considered in many applications, which can be solved by the smoothed semi-dual algorithm [Blondel et al., 2018] or the Bregman ADMM algorithm [Xu and Cheng, 2022].

Note that although the above methods are designed for the balanced OT problem (i.e., the marginal distributions of the coupling are known and with equal measures), they can be extended to the cases where the marginal distributions of the coupling are unreliable or unavailable, e.g., the unbalanced OT problem [Chizat et al., 2018] and the partial OT problem [Benamou et al., 2015]. The optimal coupling indicates the correspondence between the two distributions' samples, which has the potential to various resource assignment problems [Zhang and Zhu, 2019; Hughes and Chen, 2021]. To our surprise, however, applying the optimal transport methods is seldom considered for the SPFD problem. In the following content, we fill this blank and propose a group sparse optimal transport method to solve the SPFD problem approximately.

## 3 ProposedMethod

## 3.1 An Optimal Transport Perspective on SPFD

Given M plants and N products, we would like to design a flexible manufacturing network to match the plants' supplies with the products' random demands in a robust and efficient way. Following existing SPFD work [Chou et al., 2008], we assume that the supplies are predefined and formulated as a vector µ ∈ [0,∞)M , and the demands are unknown but can be sampled from a known distribution, i.e., ν ∈ [0,∞)N ∼ νD, where νD is the distribution of the demands defined in the sample space D. As illustrated in Figure 1, the proposed SPFD problem corresponds to constructing a bipartite graph connecting the plants with the products. Mathematically, this problem can be formulated as the following bi-level optimization problem:

<!-- formula-not-decoded -->

Here, A = [amn] ∈ {0, 1}M×N is a binary matrix indicating the network topology. amn = 1 means that there exists an edge (m, n) between the plant m and the product n. pν(A) represents the maximum profit obtained when the network topology is A and the demand vector is ν. The upper-level optimization problem corresponds to optimizing the network topology to maximize the expected profit. In the lower level, a set of max-flow problems are solved to determine the profits given the current network topology. To avoid a trivial solution (i.e., A = 1M×N), we further consider the sparsity of A, restricting its number of nonzero elements to be equal to or less than a threshold L, i.e., ∥A∥0 ≤ L, where ∥ · ∥0 is the ℓ0-norm of the matrix.

Following the work in [Chan et al., 2022], we decouple the variables of each max-flow problem as the network topology A and the transport matrix X = [xmn] ∈ RM×N . Accordingly, A ⊙ X represents the Hadamard product of the two matrices. Its (nonzero) element amnxmn indicates the number of products passing through the edge (m, n). Obviously, A ⊙X should meet the constraints provided by the supplies and the demands, i.e., A⊙X ∈ Π(µ, ν) = {T ≥ 0|T1N ≤ µ, T T 1M ≤ ν}. The objective ⟨P,A⊙X⟩ = ∑m,n pmnamnxmn is the whole profit created by the network, and P = [pmn] is the grounding profit matrix recording the profit per product created by a plant (i.e., pmn represents the profit of making the product n from the plant m).

In this study, we consider simplifying the SPFD problem in (1) based on computational optimal transport techniques. In particular, we replace A⊙X with a single functional variable T (·). Accordingly, we can reformulate (1) as follows.

<!-- formula-not-decoded -->

Here, the functional variable T(·) : D ↦→ RM×N takes a demand vector as its input and outputs a transport matrix. Accordingly, T(ν) := argmaxT∈Π(µ,ν)⟨P, T⟩ corresponds to the output determined by the input ν, which is a typical optimal transport problem. 1 The expectation of the output, i.e., Eν∈D[T(ν)], indicates the network topology, whose element represents the expected amount of products per edge. It works as a surrogate of the A in (1), and we impose the sparse constraint on it, as shown in (2).

Introducing the functional variable, we reformulate the bilevel optimization problem (1) to a nested functional optimization problem (2). The inner problem optimizes a slice T(ν) given a specific ν, and the outer problem optimizes the whole function T(·) with a sparse constraint. Suppose that we have a set of demands sampled from νD, denoted as {νk}Kk=1 , as the work in [Wei et al., 2021; Chan et al., 2022] did. We can approximate the functional variable with the help of the kernel trick [Sch¨olkopf et al., 2001], i.e., T(ν) = ∑Kk=1 Tkκ(ν, νk), where Tk := T(νk) and κ can be a predefined kernel function. Then, we can rewrite (2) as

<!-- formula-not-decoded -->

which corresponds to solving a group of K optimal transport problems with a consistency regularizer on the sparsity of the couplings. Denote the coupling group as a tensor T = [Tk] ∈ RM×N×K . Following the work intensor [Bengio et al., 2009;

1 When Π(µ, ν) = {T ≥ 0|T 1N = µ, T T 1M = ν}, the problem is a balanced OT problem [Villani, 2008]. When Π(µ, ν) = {T ≥ 0|T1N ≤ µ, T T 1M ≤ ν}, the problem is an unbalanced/partial OT problem [Benamou et al., 2015].

Deng et al., 2013], we relax the strict sparse constraint in (3) to a group sparse regularizer and derive the final group sparse optimal transport (GSOT) problem as follows.

<!-- formula-not-decoded -->

where ∥T ∥1,2 = ∑Mm=1 ∑Nn=1 √∑Kk=1 t2mnk represents the group sparse regularizer. This regularizer encourages the couplings to share the same sparse structure, whose significance is controlled by the hyperparameter α ≥ 0.

We achieve sparse process flexibility design approximately by solving the GSOT problem in (4). In particular, after obtaining the optimal couplings {T ∗ k }Kk=1 , we can determine the network topology as I(∑Kk=1 T ∗ k ), where I(·) is an indicator returning 1 when the input is nonzero and returning 0 otherwise. Compared to the original SPFD problem in (1), the GSOT problem is convex, which can be solved efficiently and have a theoretical guarantee of convergence. In the next subsection, we develop an efficient ADMM algorithm to solve this problem iteratively.

## 3.2 An ADMM-Based Optimization Algorithm

For the GSOT problem in (4), we introduce an auxiliary variable Z = [Zk] ∈ RM×N×K and a dual variable U = [Uk] ∈ RM×N×K , respectively, and reformulate it in the following augmented Lagrangian form:

<!-- formula-not-decoded -->

where ρ &gt; 0 and ∥ · ∥F represents the Frobenius norm of matrix. This problem can be solved in an ADMM framework. In particular, we update T , Z, and U iteratively. In the l-th iteration, we solve the following three subproblems.

1) Update T : Fixing current Z (l) and U (l) , we achieve this step via solving K independent optimal transport problems with ℓ2 regularizers: for k = 1, ...,K, we have

<!-- formula-not-decoded -->

Ignoring the term irrelevant to T , we can equivalently formulate (6) as a quadratically-regularized optimal transport (QROT) problem, i.e.,

<!-- formula-not-decoded -->

where C = -P + U (l) k - Z(l) k . We solve the problem based on the smooth OT algorithm proposed in [Blondel et al., 2018]. Here, we consider the following two situations.

- Balanced supply-demand cases: When ∥µ∥1 = ∥νk∥1, we impose the doubly-stochastic constraint on the coupling T , i.e., the feasible domain is Π(µ, νk) = {T ≥ 0|T1N = µ, T T 1M = νk}. In this situation,

we can rewrite the QROT problem in (7) in its smooth relaxed dual formulation:

<!-- formula-not-decoded -->

where [·]+ means setting negative elements to be zeros, and a and b are dual variables. This dual problem is unconstrained and can be solved efficiently by the L-BFGS algorithm [Liu and Nocedal, 1989]. Given the optimal solution of (8), denoted as a∗ and b∗ , we can derive the optimal coupling T (l+1) k = 1ρ [a∗1TN +1M(b∗ )T -C]+.

- Unbalanced supply-demand cases: When ∥µ∥1 ̸= ∥νk∥1, the feasible domain of the coupling becomes Π(µ, νk) = {T ≥ 0|T1N ≤ µ, T T 1M ≤ νk, and 1TM T1N = min{∥µ∥1, ∥νk∥1}}. In this situation, we relax the problem in (7) to a semi-relaxed smooth form: without the loss of generality, we assume that ∥µ∥1 &lt; ∥νk∥1, and the problem becomes

<!-- formula-not-decoded -->

This problem can be solved by the FISTA algorithm [Beck and Teboulle, 2009] or conditional gradient.

2) Update Z: Given fixed T (l+1) and U (l) , we achieve this step via solving a group sparse problem.

<!-- formula-not-decoded -->

Applying the soft-thresholding method [Bengio et al., 2009; Deng et al., 2013], the solution of this problem corresponds to M × N independent updates:

<!-- formula-not-decoded -->

where tmn = [tmnk] ∈ RK , umn = [umnk] ∈ RK and zmn = [zmnk] ∈ RK are vectors in T , U and Z.

3) Update U: In the framework of ADMM, the update of

<!-- formula-not-decoded -->

Convergence and efficiency. As aforementioned, the GSOT problem is a convex optimization problem. Applying our ADMM framework, each subproblem is convex as well. Therefore, the convergence of our optimization algorithm is guaranteed in theory. More specifically, the convergence rate is O(ϵ-1 ) in general because we apply first-order optimization algorithms to solve the subproblems [Nishihara et al., 2015]. 2 Regarding time complexity, the main bottleneck is solving K QROT problems that involve iterative updates. Fortunately, these problems are independent. We can solve them in parallel and accelerate the whole process.

2 The L-BFGS is a quasi-Newton method, which applies a firstorder strategy to approximate the second-order strategy.

## Algorithm 1 Oue GSOT method for the SPFD problem

- 1: Input: A supply vector µ and a set of demand vectors {νk}Kk=1 . The profit matrix P. The hyperparameters α

```
and ρ. The maximum number of edges L. 2: Output: Top-L edges {(mℓ, nℓ)}Lℓ=1 . 3: Initialize T = {Tk = µν T k }, Z = T , U = 0M×N×K. 4: while Not converge do 5: for k = 1 : K do 6: C = -P +Uk -Zk. 7: if ∥µ∥1 = ∥νk∥1 (Balanced case) then 8: Solve (8) via L-BFGS and get a∗ and b∗ . 9: Update Tk as 1ρ [a∗1TN +1M(b∗ )T -C]+. 10: else 11: Solve (9) via FISTA or conditional gradient and get Tk accordingly. 12: end if 13: end for 14: Update Z via (11). 15: Update U via (12). 16: end while 17: T ∗ ← T , and derive {(mℓ, nℓ)}Lℓ=1 via (13).
```

## 3.3 Sparse Network Topology Design

Repeating the above optimization steps till the objective function converges, we obtain the optimal variable, denoted as T ∗ = [T ∗ k ]. As aforementioned, the aggregation of the T ∗ k 's (i.e., T ∗ := 1K ∑Kk=1 T ∗ k ) is sparse and indicates the proposed network topology.

Besides automatically designing a sparse manufacturing network, our GSOT method can also design a network with a predefined number of edges (e.g., the L in (1)). In particular, from a statistical viewpoint, each element of T ∗ , i.e., t∗mn , approximates the expected amount of product n made by plant m. Accordingly, the expected profit created by the network is approximated as P ⊙ T ∗ , where each element pmnt∗mn corresponds to the expected profit generated by the edge (m, n). We can take P ⊙ T ∗ as the significance of the edges. As a result, we can sort the edges according to their significance and select the top-L significant edges, i.e.,

<!-- formula-not-decoded -->

The scheme of our GSOT method is shown in Algorithm 1.

## 4 Experiments

## 4.1 Implementation Details

Methods. To demonstrate the feasibility and effectiveness of our GSOT method, we compare it with representative SPFD methods on both synthetic and real-world datasets. According to the methodological categories proposed by [Chou et al., 2008], we consider the following three methods as our baselines, which represent three different strategies.

- Sampling method in [Lien et al., 2011]: This method samples edges and updates the significance of the remaining edges in an iterative and alternating manner. In each step, the sampling probability of each edge is estimated based on the solution of a max-flow problem.
- Expander method in [Chou et al., 2011]: This method applies the theory of graph expander, designing a sparse process structure by constructing a graph with a good node expansion ratio. It adds edges iteratively to improve the node expansion ratio in a greedy manner.
- Chaining method in [Jordan and Graves, 1995]: When the number of plants is equal to that of products and the supplies and demands are balanced, we further consider the classic long chain as a baseline.

All the methods, including the proposed method and the baselines, are implemented in Python. The Python Optimal Transport (POT) toolbox [Flamary et al., 2021] is used to solve the QROT problems in our method. Additionally, to verify the usefulness of our ADMM-based optimization algorithm (i.e., GSOT (Ours)), we apply the commercial optimization software Gurobi [Gurobi Optimization, 2021] to solve the proposed GSOT problem in (4) (i.e., GSOT (Gurobi)). For our GSOT method, the comparison between our ADMM-based algorithm and the Gurobi-based implementation is provided.

Datasets. We consider the symmetry between plants and products (i.e., whether M = N or not) and the balance between supplies and demands (i.e., whether ∥µ∥1 = ∥νk∥1 or not), constructing four synthetic datasets accordingly. The symmetric and balanced case corresponds to an ideal scenario used in many existing work [Simchi-Levi and Wei, 2012; Chou et al., 2008]: We set M = N; The supplies provided by the M plants are fixed and formulated as a uniform vector, while the demands of each product are sampled from independently a uniform distribution, leading to K Ndimensional vectors; we then normalize the supply vector and the demand vectors (i.e., ∥µ∥1 = ∥νk∥1 = 1). In the asymmetric cases, we set N = round(1.5M). In the unbalanced cases, we use unnormalized supply and demand vectors. For each synthetic dataset, we construct a grounding profit matrix randomly. We set M ∈ {5, 10, 15, 20} and K ∈ {10, 20, 50}, testing the scalability of different methods accordingly.

Besides the synthetic datasets, we consider the real-world fashion manufacturing dataset [Chou et al., 2014], which contains ten plants and ten products in the fashion industry. There exists a profit matrix with size 10 × 10, whose element indicates the profit per product generated by a plant. 3 The supplies and the mean and variance of demands are predefined, so we set one normalized supply vector and simulate 100 unnormalized demand vectors. Fifty demand vectors are used for training, while the remaining vectors are used for testing.

Evaluation criteria. For each dataset, we split the demand vectors into two sets. A training set is used to design the manufacturing network, and a testing set is used to simulate and evaluate the performance of the designed network. Given a dataset, we evaluate the performance of various methods in two ways: i) setting the number of edges from max{M,N} to round(2.5max{M,N}) and comparing the profits achieved by the networks designed by different methods; ii) designing a network with 2max{M,N} edges by

3 In [Chou et al., 2014], all the plants have the same profit per product, which is impractical in general. We add random noise to the rank-1 profit matrix and get a full-rank profit matrix.

![Image](https://ijcai.org/proceedings/2023/fulltext/679/image_000001_c495e1d02f5dee982f6661c35d02860bbea743700d0c79fc7229c296dfd81697.png)

- (a) Balanced symmetric case
- (b) Balanced asymmetric case (c) Unbalanced symmetric case (d) Unbalanced asymmetric case

Figure 2: The comparisons for various methods on their expected profits achieved under different network topologies. In this experiment, we set M = 15 and K = 20. The fully-flexible scenario corresponds to a complete bipartite manufacturing network.

![Image](https://ijcai.org/proceedings/2023/fulltext/679/image_000002_e6c715afeb3a4bbb25769c87b91ce3a3d6521a96a6d3559b057b328a0f1281f7.png)

Figure 3: The comparisons for various methods on their scalability. In each figure, we change the number of plants (M) and design a manufacturing network with 2 max{M, N} edges.

different methods and evaluating the performance with the increase of N (or M). In each experiment, we test each method in 10 trials and record the average and standard deviation of its performance. For our GSOT method and its ADMM-based algorithm, we i) analyze the robustness of the hyperparameters α and ρ, ii) analyze its empirical sample complexity, and iii) compare it with the Gurobi-based implementation.

## 4.2 Evaluation on Synthetic Data

Figure 2 shows the performance of different methods on four synthetic datasets. We can find that with the increase of edge numbers, the networks designed by different methods can always increase expected profits. However, our GSOT method consistently works better than the baselines in all four datasets. In particular, when designing a network with L edges, our GSOT method obtains higher expected profits than the baselines, which quickly approaches the completebipartite network's performance (with full flexibility).

Figure 3 further shows the comparisons for various methods on their scalability. In this experiment, we increase M from 5 to 20 and design a network with 2 max{M, N} edges in each scenario. Similar to Figure 2, we can find that our GSOT method consistently outperforms the baselines. Note that our GSOT method works better than the long chain in the balanced symmetric case. A potential reason for this phenomenon is that we introduce a random price matrix to assign the edges with different significance. This leads to a complicated scenario ignored by theoretical analysis [D´esir et al., 2016] - even if the supplies and demands are balanced and the plants and the products have the same amount, the long chain may not be optimal when the edges are weighted.

## 4.3 Feasibility of ADMM-Based Algorithm

According to the results in Figures 2 and 3, we can find that compared to the Gurobi-based implementation, our ADMMbased algorithm does not perform so well in the balanced cases. However, in the challenging unbalanced cases, it is comparable to even superior to the Gurobi-based implementation - as shown in Figures 2(c, d), our GSOT method converges to the full flexibility situation more quickly when it is implemented based on our ADMM framework.

Additionally, we find that in the unbalanced cases, the performance of Gurobi and our algorithm is different. Denote the expected profits achieved by Gurobi and our ADMM algorithm as pG and pA, respectively. Figure 4(a) shows that when supplies are insufficient (i.e., ∥µ∥1 &lt; ∥νk∥1) Gurobi often works better (i.e., pG ≥ pA). On the contrary, when supplies are sufficient (i.e., ∥µ∥1 &gt; ∥νk∥1), our ADMM algorithm often outperforms Gurobi. This result gives us useful insight into the selection of algorithms in practice.

Figure 4: (a) The influence of unbalanced supplies and demands on our algorithm. (b) The performance of our GSOT method with respect to K (in the balanced and symmetric case M = N = 10).

![Image](https://ijcai.org/proceedings/2023/fulltext/679/image_000003_c5af146461fe1f0fd87620ae5454ae0cb58590e2bd026e38f3d6e2dc98a3084a.png)

Table 1: The influences of hyperparameters on expected profit

| α (ρ = 1)   | α (ρ = 1)   |   0.01 |   0.1 |    1 |   10 |
|-------------|-------------|--------|-------|------|------|
| M M         | M = 5       |   1.12 |  1.14 | 1.11 | 1.13 |
| Balanced    | = 10        |   1.10 |  1.10 | 1.10 | 1.06 |
| Symmetric   | = 15        |   1.11 |  1.11 | 1.11 | 0.41 |
| M =         | 20          |   1.09 |  1.10 | 1.08 | 0.34 |
| ρ (α =      | 1)          |   0.01 |   0.1 |    1 |   10 |
| M =         | M = 5       |   0.97 |  1.13 | 1.13 | 1.12 |
| Balanced    | 10          |   0.78 |  1.13 | 1.13 | 1.11 |
| Asymmetric  | M = 15      |   0.68 |  0.47 | 1.06 | 1.10 |
| Asymmetric  | M = 20      |   0.63 |  0.56 | 1.09 | 1.07 |

Remark. Gurobi is a sophisticated commercial optimization software, while our ADMM algorithm is implemented based on Python. Refactoring our code with professional acceleration and optimization techniques can further improve the performance of our algorithm, which is left as our future work. Compared with Gurobi, our work is more friendly to developing countries because we will provide open-source code and make it free for commercial applications.

Additionally, our ADMM-based algorithm is robust to the data size and hyperparameter settings. As shown in Figure 4(b), the performance of our algorithm becomes stable given K = 20 supply-demand pairs. In other words, applying our algorithm, we merely need to sample a small set of supply-demand pairs when designing manufacturing networks. Table 1 shows the influences of two key hyperparameters (α and ρ) on the performance of our algorithm. We can find that when α ∈ [0.01, 1] and ρ ∈ [1, 10], our algorithm works well, whose performance is with good stability. It means that our algorithm is robust to the hyperparameters, and we can set them in a wide range. Empirically, we set α = ρ = 1 in our experiments.

## 4.4 Evaluation on Real-World Data

The real-world data is generated by a fashion manufacturing company with M = 10 manufacturing plants. The company aims to design a flexible and sparse manufacturing network to produce N = 10 different styles of parkas. For the plant m, the profit of making one product n is denoted as pmn, and all the profits are stored in a profit matrix. Here, the supplies of the plants are fixed, and the demands of the products obey different independent distributions. This problem corresponds to an unbalanced symmetric SPFD problem.

Table 2: Comparison for various methods on real-world data

| Method           | Expected Profit - 20 Edges   | Expected Profit - 23 Edges   | Expected Profit - 25 Edges   |
|------------------|------------------------------|------------------------------|------------------------------|
| Sampling         | 0.501±0.009                  | 0.529±0.001                  | 0.520±0.001                  |
| Expander         | 0.522±0.010                  | 0.532±0.008                  | 0.535±0.007                  |
| GSOT (Ours)      | 0.502±0.008                  | 0.535±0.011                  | 0.541±0.005                  |
| GSOT (Gurobi)    | 0.528±0.004                  | 0.543±0.000                  | 0.544±0.002                  |
| Full Flexibility | 0.548                        | 0.548                        | 0.548                        |

Table 2 shows the expected profits obtained by different methods. Following the experiments on synthetic data, for each method, we set the number of the corresponding network in the range [2M, round(2.5M)], which corresponds to the networks with edges equal to and more than a long chain, respectively. The results show that when the network is as sparse as a long chain, our GSOT method is comparable to the baselines. When further increasing the number of edges, our GSOT method outperforms the baselines consistently and approaches the performance of a complete bipartite network. Our ADMM-based algorithm is comparable to Gurobi - by na¨ıve Python implementation, its performance is just slightly worse than Gurobi, demonstrating its potential for practical scenarios. In the future, we plan to improve our algorithm by optimizing and adjusting its hyperparameter setting.

## 5 Conclusion

In this study, we have proposed a GSOT method for sparse process flexibility design, which provides a new solution to this significant OR problem. Instead of applying the traditional two-stage design strategy, our method constructs and evaluates the designed network in a joint optimization framework, which applies to unbalanced and asymmetric SPFD problems. We develop an effective optimization algorithm to achieve our GSOT method, leading to encouraging experimental results. In the future, we will study the theory of our method and consider improving its computational efficiency further. Moreover, we would like to find partners in the industry and test our method in real-world cases.

## Acknowledgements

Dr. Dixin Luo was supported by the National Natural Science Foundation of China (No. 62102031) and the Young Scholar Program (No. XSQD-202107001) from Beijing Institute of Technology. Dr. Hongteng Xu was supported by the National Natural Science Foundation of China (No. 62106271, 92270110), CAAI-Huawei MindSpore Open Fund, the Fundamental Research Funds for the Central Universities, and the Research Funds of Renmin University of China. He thanks the support from the Beijing Key Laboratory of Big Data Management and Analysis Methods, the Intelligent Social Governance Platform, Major Innovation &amp; Planning Interdisciplinary Platform for the "Double-First Class" Initiative.

## References

1. [Agueh and Carlier, 2011] Martial Agueh and Guillaume Carlier. Barycenters in the Wasserstein space. SIAM Jour- nal on Mathematical Analysis, 43(2):904-924, 2011.
2. [Altschuler et al., 2017] Jason Altschuler, Jonathan Weed, and Philippe Rigollet. Near-linear time approximation algorithms for optimal transport via sinkhorn iteration. In Advances in Neural Information Processing Systems, pages 1964-1974, 2017.
3. [Arjovsky et al., 2017] Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasserstein generative adversarial net- works. In International Conference on Machine Learning, pages 214-223, 2017.
4. [Barrett et al., 2020] Thomas Barrett, William Clements, Jakob Foerster, and Alex Lvovsky. Exploratory combi- natorial optimization with reinforcement learning. In Pro- ceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 3243-3250, 2020.
5. [Beck and Teboulle, 2009] Amir Beck and Marc Teboulle. A fast iterative shrinkage-thresholding algorithm for lin- ear inverse problems. SIAM journal on imaging sciences, 2(1):183-202, 2009.
6. [Benamou et al., 2015] Jean-David Benamou, Guillaume Carlier, Marco Cuturi, Luca Nenna, and Gabriel Peyré. Iterative Bregman projections for regularized transporta- tion problems. SIAM Journal on Scientific Computing, 37(2):A1111-A1138, 2015.
7. [Bengio et al., 2009] Samy Bengio, Fernando Pereira, Yoram Singer, and Dennis Strelow. Group sparse coding. Advances in neural information processing systems, 22, 2009.
8. [Blondel et al., 2018] Mathieu Blondel, Vivien Seguy, and Antoine Rolet. Smooth and sparse optimal transport. In In- ternational conference on artificial intelligence and statis- tics, pages 880-889. PMLR, 2018.
9. [Boyd et al., 2011] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, Jonathan Eckstein, et al. Distributed op- timization and statistical learning via the alternating direc- tion method of multipliers. Foundations and Trends® in Machine learning, 3(1):1-122, 2011.
10. [Chan et al., 2022] Timothy CY Chan, Daniel Letourneau, and Benjamin G Potter. Sparse flexible design: a machine learning approach. Flexible Services and Manufacturing Journal, pages 1-51, 2022.
11. [Chen et al., 2015] Xi Chen, Jiawei Zhang, and Yuan Zhou. Optimal sparse designs for process flexibility via prob- abilistic expanders. Operations Research, 63(5):1159- 1176, 2015.
12. [Chen et al., 2019] Xi Chen, Tengyu Ma, Jiawei Zhang, and Yuan Zhou. Optimal design of process flexibility for gen- eral production systems. Operations Research, 67(2):516- 531, 2019.
13. [Chizat et al., 2018] Lenaic Chizat, Gabriel Peyré, Bernhard Schmitzer, and Franc ¸ois-Xavier Vialard. Scaling algo- rithms for unbalanced optimal transport problems. Math- ematics of Computation, 87(314):2563-2609, 2018.
14. [Chou et al., 2008] Mabel C Chou, Chung-Piaw Teo, and Huan Zheng. Process flexibility: design, evaluation, and applications. Flexible Services and Manufacturing Jour- nal, 20:59-94, 2008.
15. [Chou et al., 2010] Mabel C Chou, Geoffrey A Chua, Chung-Piaw Teo, and Huan Zheng. Design for process flexibility: Efficiency of the long chain and sparse struc- ture. Operations research, 58(1):43-58, 2010.
16. [Chou et al., 2011] Mabel C Chou, Geoffrey A Chua, Chung-Piaw Teo, and Huan Zheng. Process flexibility re- visited: The graph expander and its applications. Opera- tions research, 59(5):1090-1105, 2011.
17. [Chou et al., 2014] Mabel C Chou, Geoffrey A Chua, and Huan Zheng. On the performance of sparse process struc- tures in partial postponement production systems. Opera- tions research, 62(2):348-365, 2014.
18. [Cuturi and Doucet, 2014] Marco Cuturi and Arnaud Doucet. Fast computation of Wasserstein barycenters. In International Conference on Machine Learning, pages 685-693, 2014.
19. [Cuturi, 2013] Marco Cuturi. Sinkhorn distances: Light- speed computation of optimal transport. In Advances in neural information processing systems, pages 2292-2300, 2013.
20. [Deng et al., 2013] Wei Deng, Wotao Yin, and Yin Zhang. Group sparse optimization by alternating direction method. In Wavelets and Sparsity XV, volume 8858, pages 242-256. SPIE, 2013.
21. [Désir et al., 2016] Antoine Désir, Vineet Goyal, Yehua Wei, and Jiawei Zhang. Sparse process flexibility designs: Is the long chain really optimal? Operations Research, 64(2):416-431, 2016.
22. [Feng et al., 2017] Wancheng Feng, Chen Wang, and Zuo- Jun Max Shen. Process flexibility design in heterogeneous and unbalanced networks: A stochastic programming ap- proach. IISE Transactions, 49(8):781-799, 2017.
23. [Flamary et al., 2021] Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Z. Alaya, Aurélie Bois- bunon, Stanislas Chambon, Laetitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, Léo Gautheron, Nathalie T.H. Gayraud, Hicham Janati, Alain Rakotoma- monjy, Ievgen Redko, Antoine Rolet, Antony Schutz, Vivien Seguy, Danica J. Sutherland, Romain Tavenard, Alexander Tong, and Titouan Vayer. Pot: Python opti- mal transport. Journal of Machine Learning Research, 22(78):1-8, 2021.
24. [Frogner et al., 2015] Charlie Frogner, Chiyuan Zhang, Hos- sein Mobahi, Mauricio Araya, and Tomaso A Poggio. Learning with a wasserstein loss. In Advances in Neural Information Processing Systems, pages 2053-2061, 2015.
25. [Graves and Tomlin, 2003] Stephen C Graves and Brian T Tomlin. Process flexibility in supply chains. Management Science, 49(7):907-919, 2003.
26. [Gurobi Optimization, 2021] LLC Gurobi Optimization. Gurobi optimizer reference manual, 2021.
27. [Harris and Ross, 1956] TE Harris and FS Ross. Fundamen- tals of a method for evaluating rail net capacities. Research Memorandum, 1573, 1956.
28. [Hughes and Chen, 2021] Jason Hughes and Juntao Chen. Fair and distributed dynamic optimal transport for re- source allocation over networks. In 2021 55th Annual Conference on Information Sciences and Systems (CISS), pages 1-6. IEEE, 2021.
29. [Jordan and Graves, 1995] William C Jordan and Stephen C Graves. Principles on the benefits of manufacturing pro- cess flexibility. Management science, 41(4):577-594, 1995.
30. [Khalil et al., 2017] Elias Khalil, Hanjun Dai, Yuyu Zhang, Bistra Dilkina, and Le Song. Learning combinatorial op- timization algorithms over graphs. Advances in neural in- formation processing systems, 30, 2017.
31. [Kusner et al., 2015] Matt Kusner, Yu Sun, Nicholas Kolkin, and Kilian Weinberger. From word embeddings to docu- ment distances. In International conference on machine learning, pages 957-966. PMLR, 2015.
32. [Lien et al., 2011] Robert W Lien, Seyed MR Iravani, Karen Smilowitz, and Michal Tzur. An efficient and robust de- sign for transshipment networks. Production and Opera- tions Management, 20(5):699-713, 2011.
33. [Liu and Nocedal, 1989] Dong C Liu and Jorge Nocedal. On the limited memory bfgs method for large scale optimiza- tion. Mathematical programming, 45(1-3):503-528, 1989.
34. [Maretic et al., 2022] Hermina Petric Maretic, Mireille El Gheche, Matthias Minder, Giovanni Chierchia, and Pascal Frossard. Wasserstein-based graph alignment. IEEE Transactions on Signal and Information Processing over Networks, 8:353-363, 2022.
35. [Nishihara et al., 2015] Robert Nishihara, Laurent Lessard, Ben Recht, Andrew Packard, and Michael Jordan. A gen- eral analysis of the convergence of admm. In International conference on machine learning, pages 343-352. PMLR, 2015.
36. [Peyré et al., 2019] Gabriel Peyré, Marco Cuturi, et al. Com- putational optimal transport. Foundations and Trends® in Machine Learning, 11(5-6):355-607, 2019.
37. [Schölkopf et al., 2001] Bernhard Schölkopf, Ralf Herbrich, and Alex J Smola. A generalized representer theorem. In Proceedings of the 14th Annual Conference on Computa- tional Learning Theory and and 5th European Conference on Computational Learning Theory, pages 416-426, 2001.
38. [Simchi-Levi and Wei, 2012] David Simchi-Levi and Yehua Wei. Understanding the performance of the long chain and sparse designs in process flexibility. Operations research, 60(5):1125-1141, 2012.
39. [Simchi-Levi et al., 2018] David Simchi-Levi, He Wang, and Yehua Wei. Increasing supply chain robustness through process flexibility and inventory. Production and Operations Management, 27(8):1476-1491, 2018.
40. [Simon et al., 2013] Noah Simon, Jerome Friedman, Trevor Hastie, and Robert Tibshirani. A sparse-group lasso. Jour- nal of computational and graphical statistics, 22(2):231- 245, 2013.
41. [Sinkhorn and Knopp, 1967] Richard Sinkhorn and Paul Knopp. Concerning nonnegative matrices and doubly stochastic matrices. Pacific Journal of Mathematics, 21(2):343-348, 1967.
42. [Sultana et al., 2022] Nasrin Sultana, Jeffrey Chan, Tabinda Sarwar, and AK Qin. Learning to optimise general tsp instances. International Journal of Machine Learning and Cybernetics, 13(8):2213-2228, 2022.
43. [Titouan et al., 2019] Vayer Titouan, Nicolas Courty, Ro- main Tavenard, and Rémi Flamary. Optimal transport for structured data with application on graphs. In Inter- national Conference on Machine Learning, pages 6275- 6284. PMLR, 2019.
44. [Villani, 2008] Cédric Villani. Optimal transport: Old and new, volume 338. Springer Science & Business Media, 2008.
45. [Wang and Banerjee, 2014] Huahua Wang and Arindam Banerjee. Bregman alternating direction method of mul- tipliers. In Advances in Neural Information Processing Systems, pages 2816-2824, 2014.
46. [Wei et al., 2021] Yehua Wei, Lei Zhang, Ruiyi Zhang, Shi- jing Si, Hao Zhang, and Lawrence Carin. Reinforcement learning for flexibility design problems. arXiv preprint arXiv:2101.00355, 2021.
47. [Xie et al., 2020] Yujia Xie, Xiangfeng Wang, Ruijia Wang, and Hongyuan Zha. A fast proximal point method for com- puting exact wasserstein distance. In Uncertainty in artifi- cial intelligence, pages 433-453. PMLR, 2020.
48. [Xu and Cheng, 2022] Hongteng Xu and Minjie Cheng. Regularized optimal transport layers for gener- alized global pooling operations. arXiv preprint arXiv:2212.06339, 2022.
49. [Yan et al., 2018] Zhenzhen Yan, Sarah Yini Gao, and Chung Piaw Teo. On the design of sparse but ef- ficient structures in operations. Management Science, 64(7):3421-3445, 2018.
50. [Ye et al., 2017] Jianbo Ye, Panruo Wu, James Z Wang, and Jia Li. Fast discrete distribution clustering using wasser- stein barycenter with sparse support. IEEE Transactions on Signal Processing, 65(9):2317-2332, 2017.
51. [Zhang and Zhu, 2019] Rui Zhang and Quanyan Zhu. Consensus-based distributed discrete optimal transport for decentralized resource matching. IEEE Transactions on Signal and Information Processing over Networks, 5(3):511-524, 2019.