Simplicial decomposition of variational inequalities with multiple nonlinear column generation

1.   Introduction

Column generation (CG) is a method to solve large-scale mathematical programming problems such as optimization models and variational inequalities (Ⅵ) involving a vast number of variables. These models can be found in multi-cloud systems, multi-commodity economic equilibrium models, and transportation models. These models would become too large to be solved due to the rise of big data. For example, large-scale VI problems can be easily found in stochastic VI settings with large-scale scenario analysis. Even a well-established VI solver (PATH) has been found to terminate when attempting to solve a large-scale VI problem due to a lack of memory ^[1]. Using column generation algorithms, such large-scale problems can be decomposed into smaller subproblems to be solved on several computers. Another motivation for using CG is to improve model development and maintenance efficiency by joining separately well-developed submodels when a converged solution is needed. Other research groups, such as Murphy et al. ^[2], have made similar observations. These smaller submodels can even reside on different computers connected by a network. On top of improving model development and maintenance, CG may improve computational efficiency. For example, Chung ^[3] reported that some computational instances of CG were faster than the reference method, the PATH solver. Chung ^[3] also reported that the PATH solver terminated when attempting to solve a large-scale VI problem because of insufficient memory.

However, CG has a major drawback which is the long-tail convergence property, called long-tail effect. CG involves solving a subproblem and a master problem iteratively. An example of CG is simplicial decomposition (SD) method. Simplicial decomposition was recently employed and modified for different application aspects, such as Bettiol et al. ^[4], Uciński ^[5], Guignard and Ahlatcioglu ^[6], and Delle Site ^[7]. For the multiple subproblems aspect, see Morabit et al. ^[8]. It is assumed that there are efficient solution methods to solve the subproblem and the master problem. However, it is well known that the SD converge rapidly at first but subsequently slowly, with a long tail of near-optimal solutions. As the final optimality gap cannot be closed, this long-tail convergence resulting in poor computational performance becomes a drawback. Other drawbacks are dual oscillations and primal degeneracy, and alternative dual optimal solutions in optimization models. In my experience, CG methods for VI problems inherit this long-tail convergence property. That is, CG-VI is likely to approach the equilibrium solution rapidly at first but then tail off, making only a tiny amount of progress per CG iteration. Consequently, the time savings afforded by the rapid initial convergence tend to be offset by the tailing off. I develop a new method to alleviate the long-tail effect of CG-VI in the current paper.

1.1.   Literature review of the methods of resolving the tailing-off effect

1.1.1.   Existing methods of resolving the tailing-off effect for LP

Implementation enhancement: Nazareth ^[9] used a set of LP problem examples consisting of a few constraints and variables to identify the numerical difficulties that can generally occur when applying the DW method. Nazareth showed that even when stable techniques, such as the stable refactorization of the corresponding basis matrix, are used, the columns of the computed master program can differ substantially from those of the true one. Later, Nazareth ^[10] noted that because changes in primal and dual variables are applied iteratively, a certain amount of "cleaning up" is to be expected. According to Ho ^[11] explanation of long-tail convergence, implicit error bounds may imply a tolerance for convergence below which further apparent improvements should be considered as noise. However, Lübbecke and Desrosiers ^[12] remarked that tailing off also occurs when columns are computed exactly, e.g., through the use of combinatorial algorithms.

Column dropping: One may anticipate that dropping columns from the master problem will improve computational efficiency. Dantzig ^[13] proposed some methods for deciding how long proposals should be kept in the master problem, such as deleting all non-basic proposals immediately (apart from those just generated). Indeed, it is not entirely clear when to delete proposals from the master problem. Beale et al. ^[14] reported that when using Dantzig's ^[13] above-mentioned methods, for a problem with 450 constraints, the DW method performed marginally better than the simplex method. For NLP, Murphy ^[15] provided two conditions under which columns may be dropped from the restricted master problem. However, O'Neill's ^[16] computational results revealed virtually no differences in the measures of computational efficiency when all columns were retained and when only basic columns were retained.

Stabilized column generation methods: Another approach is to use stabilized column generation methods to ensure that the dual variable values of the linking constraints smoothly converge to their respective optima without vehement oscillation. From the literature, three stabilization techniques of optimization models can be found. These techniques are the proximity of a stability center, smoothing techniques, and centralized prizes for stabilizing the iterative dual solutions from the master problems. There are some well-known stabilization principles, such as the Boxtep method developed by Marsten et al. ^[17]. These stabilization methods are widely used in branch-and-bound algorithms for mixed integer programs. It should be noted that these techniques rely on solving dual master optimization problems. It would become a big challenge when one tries to use the dual master VI problems.

1.1.2.   Existing methods of resolving the tailing-off effect for VI

To resolve the tailing-off problem afflicting simplicial decomposition of VI (SD-VI), Larsson and Patriksson ^[18] considered the use of subproblem column generations of nonlinear simplicial decomposition (NSD) to reduce the number of decomposition iterations. The generalized NSD method is obtained from the restricted simplicial decomposition method by replacing the linear column generation subproblem with a nonlinear column generation (NCG) subproblem. Generating columns based on the NCG subproblem requires fewer columns to describe an optimal solution, resulting in fewer decomposition iterations. There are three kinds of NSD column generation subproblems addressed by adopting Newton's method, the diagonalized Newton algorithm, and the projection method.

Chung et al. ^[19] extended the application of the DW procedure from linear programming (LP) to VI problems. The master problem and the subproblems of the decomposed VI are VIs. To shorten the computation time, Chung and Fuller ^[20] derived an approximation method to solve the subproblem VIs. Chung ^[3] integrated the results of Çelebi and Fuller ^[1] and Chung and Fuller ^[20] to develop an approximation method to solve DW-VI, in which both master VI and subproblem VI can be solved approximately as DW for LP. Chung ^[21] derived a truncated DW-VI method, in which the subproblem VI was approximately solved by one iteration of DW. Although they did not discuss the tailing-off effect, their computational results showed that the tailing-off effect for VI was alleviated.

In recent years, in the context of optimization problems, some research works concentrated on using machine learning to predict the stabilization technique for column generation and the patterns of different existing columns for re-optimization, see Morabit et al. ^[22] and Kraul et al. ^[23]. Although these approaches do not focus on column generation, their results can help to overcome the tailing-off effect in the re-optimization processes.

1.2.   Approaches and contributions

In the context of Dantzig-Wolfe of LP, Dirickx and Jennergren ^[24] reported that, in general, the more proposals are used to initialize the algorithm, the faster the solution can be found. Therefore, in the current paper, I propose to solve multiple column generation subproblems in each SD-VI iteration instead of solving only one subproblem. By solving multiple column generation subproblems in each SD-VI iteration, I obtain two advantages, as follows: (1) Different columns will be generated and stored in the master problem to improve the representativeness of the corresponding convex hull. (2) Generating multiple column generation subproblem solutions in each SD-VI iteration will enable the corresponding convex hull to be rapidly enlarged. Consequently, the number of SD-VI iterations can be greatly reduced. The idea of generating multiple columns at a given iteration is not a new strategy in nonlinear programming, see Morabit et al. ^[8]. In short, I propose a new method, simplicial decomposition with multiple subproblems for VI, that solves multiple column generation subproblem in each SD-VI iteration to reduce the number of decomposition steps and thereby ease the tailing-off effect. One of the challenges of the new method is how to derive "multiple" subproblems. According to my knowledge, there is no research on solving multiple subproblems in SD-VI. It is a research gap in the literature. Moreover, the proposed method can also apply to NLP since the VI framework includes NLP.

While employing multiple subproblems (multiple column generation) in SD of optimization models, convergence properties and theories may be easily derived based on the existing results of the monotonic convergence properties. However, the convergence property of the SD-VI is non-monotone. Hence, the convergence properties and theories of multiple column generation of SD-VI is derived, and empirical tests are given.

In brief, the main contribution of the current paper is to derive a new method, simplicial decomposition with multiple subproblems for VI. The major challenges are (1) to derive "multiple subproblems" in SD-VI in the Subsection 3.2 and (2) to prove the convergence of the new method.

Section 2 provides the background of the VI problems, SD-VI methods, and SD of NLP with nonlinear column generation (NCG). Section 3 derives SD-VI with NCG (SD-VI-NCG) and with multiple NCGs (SD-VI-mNCG). SD-VI-NCG and SD-VI-mNCG are SD with one subproblem and multiple subproblems for VI, respectively. Convergence properties of the SD-VI-NCG and SD-VI-mNCG are also included. Section 4 is used to report computational performance of SD-NLP-mNCG, SD-VI-mNCG, and their re-optimization processes. Section 5 provides the conclusion with further research topics.

2.   Background

The class of the VI problems I will study can be described as follows. $VI\left(G, K\right)$: Find a vector ${x}^{*}\in K\subseteq {R}^{n}$ such that $G\left({x}^{*}\right)^{T}\left(x-{x}^{*}\right)\ge 0 \;\forall x\in K$, where G is given continuous mapping from $Kto{R}^{n}$, superscript T denotes the transpose, and all vectors are considered to be column vectors. K is a nonempty, closed, and convex set. Applications of $VI\left(G, K\right)$ can be found in user equilibrium traffic assignment problems and energy equilibrium problems. For large-scale VI, see Murphy et al. ^[2]. Harker and Pang ^[25] provided standard conditions for the existence and uniqueness of solutions to this class of VI problems. There are different methods of solving $VI\left(G, K\right)$, such as PIES-like methods or Newtonian methods.

2.1.   Column generation of VI methods (CG-VI)

There are two popular column generation of VI methods, simplicial decomposition of VI (see Lawphongpanich and Hearn ^[26]) and Dantzig-Wolfe of VI (see Fuller and Chung ^[27]). Since the my results are derived from simplicial decomposition of VI (SD-VI), a brief description of SD-VI is given below. Moreover, multiple subproblems are generated from SD with nonlinear column generation (NCG) of nonlinear programming (Larsson et al. ^[28]), NCG is given in this subsection.

2.2.   Simplicial decomposition of VI (SD-VI)

Lawphongpanich and Hearn ^[26] presented a type of simplicial decomposition of VI method for the user equilibrium traffic assignment problem, in which $K$ is a closed convex set of feasible flow patterns, and $G$ is a cost mapping. The following SD-VI is not exactly the one presented in Lawphongpanich and Hearn ^[26]. The major differences are that there is no pre-set convergence sequence ${(\varepsilon }^{k})$ for convergent of their SD-VI and no column dropping consideration for improving the computational performance of the master-VI problems. It is because my results of employing multiple subproblems do not require the pre-set convergence sequence and column dropping procedure.

Note that SD-VI consists of linear programming subproblem, $Sub-{LP}^{k}$, and an equilibrium master problem, ${\mathrm{M}\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}-\mathrm{V}\mathrm{I}}^{k}$, where superscript $k$ represents the iteration number of CG-VI. $Sub-{LP}^{k}$ is defined as follows:

where ${x}_{M}^{k}$ equals the solution of master problem (${\mathrm{M}\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}-\mathrm{V}\mathrm{I}}^{\mathrm{k}}$), described below.

The feasible set for the master problem at iteration $k$ is restricted to all convex combinations of the $k$ proposals (solutions of the subproblem). Let ${X}^{k}$ be the $n\times k$ matrix whose columns are the solutions ${x}_{S}^{i}$ of the subproblems solved at iterations $i = 1, \dots, k$; i.e., ${X}^{k} = \left[{x}_{S}^{1}, {x}_{S}^{2}, \dots, {x}_{S}^{k}\right]$. The weights on the proposals in the convex combination are contained in the vector $\lambda \in {R}^{k}$. The feasible set for the master problem is defined as

where ${e}^{k}\in {R}^{k}$ is a vector with 1 for every entry. The master problem at iteration $k$ is defined as

Alternatively, let $conv\left({X}^{k}\right)$ represent the set of all convex combinations of the columns of ${X}^{k}$. Then, I have the following alternative master problem at iteration k, defined as

With the convergence gap, ${CG}^{k} = {G\left({x}_{M}^{k}\right)}^{T}\left({x}_{S}^{k}-{x}_{M}^{k}\right)$, for the stopping condition, I can define the standard SD-VI algorithm as follows:

[SD-VI]

Step 0: Set $k = 0$. Choose $\varepsilon > 0$, ${x}_{M}^{1}\in K$, and ${X}^{0}$ is a null matrix.

Step 1: Increment $k\leftarrow k+1$.

Solve $Sub-{LP}^{k}\left(G\left({x}_{M}^{k}\right), K\right)$ and place the solution ${x}_{S}^{k}$ in the matrix ${X}^{k} = \left[{X}^{k-1}, {x}_{S}^{k}\right]$.

If $k = 1$ then go to Step 2; if ${CG}^{k}\ge -\varepsilon$ then STOP; else go to Step 2.

Step 2: Solve ${\mathrm{M}\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}-\mathrm{V}\mathrm{I}}^{k}\left(G, {K}^{k}\right).$ Record $G\left({x}_{M}^{k}\right)$. Go to Step 1.

In the next subsection, I describe the result of Larsson et al. ^[28] concerning how the nonlinear column generation (NCG) is employed in the context of SD of NLP. Based on their results, I first derive a SD of VI with NCG, then I derive a SD of VI with multiple NCG.

2.3.   Simplical decomposition of NLP with nonlinear column generation (NCG)

According to Larsson et al. ^[28], SD of NLP with NCG (SD-NLP-NCG) is the combination of the SD principle and the nonlinear search direction finding subproblem of a primal descent algorithm for NLP, not VI.

For NLP, the master problem and the subproblem of SD are shown below:

For NCG, the above subproblem, $Sub-{LP}^{k}$, becomes the combination of the SD principle and the nonlinear search direction finding subproblem of a primal descent algorithm for NLP. Larsson et al. ^[28] modified the subproblem by adding a continuous regularization function, $\varphi \left(x, x\right):K\times K⟼R$, with the following properties:

(1) $\varphi \left(\cdot, y\right)$ is strictly convex and continuously differentiable on $K$ for every $y\in K$;

(2) ${\nabla }_{x}\varphi \left(x, y\right) = 0$ holds if and only if $x = y$.

Then, given an iterate ${x}_{M}^{k}$, the subproblem of SD-NLP-NCG becomes an NLP.

It is noted that the solution of ${Sub}_{NCG}-{NLP}^{k}$, ${x}_{S}^{k}$, is unique due to the compactness of $K$ and the strict convexity of $\varphi \left(\cdot, y\right)$. Moreover, ${x}_{M}^{k-1}$ is optimal in NLP if and only if ${x}_{M}^{k-1}$ solves ${Sub}_{NCG}-{NLP}^{k}$, that is, ${x}_{S}^{k} = {x}_{M}^{k-1}$. In the numerical experiments of Larsson et al. ^[28], the function $\varphi$ had the form:

where $Q\left(y\right)$ is a positive definite and symmetric matrix for every $y$ in $K$. Three matrixes were proposed: $Q\left(y\right)\equiv \gamma I, \gamma > 0$, $Q\left(y\right) = \mathit{\nabla} G\left({x}_{M}^{k}\right)$, and $Q\left(y\right) = \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\mathit{\nabla }G\left({x}_{M}^{k}\right)$.

3.   Simplicial decomposition of VI and NCG

Based on the results of Larrsson et al. ^[28], I first define the method of SD-VI with a NCG (SD-VI-NCG), then SD-VI with multiple NCG (SD-VI-mNCG).

3.1.   Simplicial decomposition of VI with NCG (SD-VI-NCG)

To derive a SD-VI with NCG, the $\boldsymbol{S}\boldsymbol{u}\boldsymbol{b}-{\boldsymbol{L}\boldsymbol{P}}^{\boldsymbol{k}}\left(\boldsymbol{G}\left({\boldsymbol{x}}_{\boldsymbol{M}}^{\boldsymbol{k}}\right), \boldsymbol{K}\right)$ of SD-VI is replaced with ${\boldsymbol{S}\boldsymbol{u}\boldsymbol{b}}_{\boldsymbol{N}\boldsymbol{C}\boldsymbol{G}}-{\boldsymbol{N}\boldsymbol{L}\boldsymbol{P}}^{\boldsymbol{k}}\left(\boldsymbol{G}\left({\boldsymbol{x}}_{\boldsymbol{M}}^{\boldsymbol{k}}\right), \boldsymbol{K}\right)$. Then, I can define the standard SD-VI-NCG algorithm as follows:

[SD-VI-NCG]

Step 0: Set $k = 0$. Choose $\varepsilon > 0$, ${x}_{M}^{1}\in K$, and ${X}^{0}$ is a null matrix.

Step 1: Increment $k\leftarrow k+1$.

Solve $\widetilde{S u b}-{NLP}^{k}\left(\tilde{G}\left(x; {x}_{M}^{k}\right), K\right)$ and place the solution ${x}_{S}^{k}$ in the matrix ${X}^{k} = \left[{X}^{k-1}, {x}_{S}^{k}\right]$.

If $k = 1$ then go to Step 2; if ${CG}^{k}\ge -\varepsilon$ then STOP; else go to Step 2.

Step 2: Solve ${\boldsymbol{M}\boldsymbol{a}\boldsymbol{s}\boldsymbol{t}\boldsymbol{e}\boldsymbol{r}-\boldsymbol{V}\boldsymbol{I}}^{\boldsymbol{k}}\left(\boldsymbol{G}, {\boldsymbol{K}}^{\boldsymbol{k}}\right).$ Record $G\left({x}_{M}^{k}\right)$. Go to Step 1.

Noted that if I use $\tilde{G}\left(x; {x}_{M}^{k}\right) = G\left({x}_{M}^{k}\right)+Q\left({x}_{M}^{k}\right)\left(x-{x}_{M}^{k}\right)$, the subproblem is equal to the one in SD-NLP-NCG, ${\boldsymbol{S}\boldsymbol{u}\boldsymbol{b}}_{\boldsymbol{N}\boldsymbol{C}\boldsymbol{G}}-{\boldsymbol{N}\boldsymbol{L}\boldsymbol{P}}^{\boldsymbol{k}}\left(\boldsymbol{G}\left({\boldsymbol{x}}_{\boldsymbol{M}}^{\boldsymbol{k}}\right), \boldsymbol{K}\right)$. That is,

3.1.1.   Convergence of SD-VI-NCG

With the assumption of $\tilde{G}\left(x; {x}_{M}^{k}\right)$ is strictly monotone in $x$ and $\tilde{G}\left({x}_{M}^{k}; {x}_{M}^{k}\right) = G\left({x}_{M}^{k}\right)$, I can extend the results of Chung and Fuller ^[20] to have the following properties.

Theorem 1a. ${\lambda }^{k}$ solves ${Master-VI}^{k}\left(G, {\mathrm{\Lambda }}^{k}\right)$ iff there exist ${\lambda }^{k}\in {R}_{+}^{k}$ and ${\theta }^{k}\in R$ such that the following relations hold:

Proof. This is a standard result for VI problems. See, e.g., Proposition 2.2 in Harker and Pang ^[25].

Theorem 2a. Given the property $\tilde{G}\left({x}_{M}^{k}; {\xi }^{k}\right) = G\left({x}_{M}^{k}\right)$. If ${x}_{M}^{k}$ solves $\widetilde{S u b}-{NLP}^{k}\left(\tilde{G}\left(x; {x}_{M}^{k}\right), K\right)$, then ${x}_{M}^{k}$ solves $VI\left(G, K\right)$.

Proof. Ref. new Theorem 4 in Chung and Fuller ^[20].

Suppose that ${x}_{M}^{k} = {X}^{k}{\lambda }^{k}$ solves $\widetilde{S u b}-{NLP}^{k}\left(\tilde{G}\left(x; {x}_{M}^{k}\right), K\right)$. It follows that ${\tilde{G}\left({x}_{M}^{k}; {x}_{M}^{k}\right)}^{T}\left(x-{x}_{M}^{k}\right)\ge 0 \;\forall x\in K$. As $\tilde{G}\left({x}_{M}^{k}; {x}_{M}^{k}\right) = G\left({x}_{M}^{k}\right)$, I have ${G\left({x}_{M}^{k}\right)}^{T}\left(x-{x}_{M}^{k}\right)\ge 0 \;\forall x\in K$ and I may conclude that ${x}_{M}^{k}$ solves $VI\left(G, K\right)$.

Theorem 3a. Assume that $\tilde{G}\left(x; {\xi }^{k}\right)$ is strictly monotone in $x$. Given the property $\tilde{G}\left({x}_{M}^{k}; {\xi }^{k}\right) = G\left({x}_{M}^{k}\right)$. If ${CG}^{k} = {G\left({x}_{M}^{k}\right)}^{T}\left({x}_{S}^{k}-{x}_{M}^{k}\right)\ge 0$, then ${x}_{M}^{k-1}$ solves $VI\left(G, K\right)$.

Proof. Ref. new Theorem 6(a) of Chung and Fuller ^[20].

We shall show that if ${x}_{M}^{k}$ does not solve $VI\left(G, K\right)$, then ${CG}^{k} < 0$. By Theorem 2a, ${x}_{M}^{k}\ne {x}_{S}^{k}$, and since $\tilde{G}\left({x}_{M}^{k}; {\xi }^{k}\right) = G\left({x}_{M}^{k}\right)$, strict monotonicity of $\tilde{G}\left(x; {\xi }^{k}\right)$ implies that ${\left(G\left({x}_{M}^{k}\right)-\tilde{G}\left({x}_{S}^{k}; {\xi }^{k}\right)\right)}^{T}\left({x}_{M}^{k}-{x}_{S}^{k}\right) > 0$. Since ${x}_{S}^{k}$ solves $\widetilde{S u b}-{NLP}^{k}\left(\tilde{G}\left(x; {x}_{M}^{k}\right), K\right)$, it follows that ${\tilde{G}\left({x}_{S}^{k}; {\xi }^{k}\right)}^{T}\left({x}_{M}^{k}-{x}_{S}^{k}\right)\ge 0$. Adding this last inequality to the strict inequality and multiply by -1 yields ${G\left({x}_{M}^{k}\right)}^{T}\left({x}_{S}^{k}-{x}_{M}^{k}\right) < 0$.

Theorem 4. Either ${CG}^{k}\ge 0$ at a finite iteration number $k$, or ${CG}^{k} < 0$ for all iterations $k$. In the latter case, any infinite subsequence of ${\left\{\left({x}_{M}^{k}, {x}_{S}^{k}\right)\right\}}_{k = 1}^{\infty }$ has at least one limit point, and if $G$ is continuous then ${lim}_{k\to \infty }{CG}^{k} = 0$.

Proof. Ref. new Theorem 8 of Chung and Fuller ^[20].

Suppose that ${CG}^{k}{ = G\left({x}_{M}^{k}\right)}^{T}\left({x}_{S}^{k}-{x}_{M}^{k}\right) < 0$ for all $k$ and suppose, contrary to my desired conclusion, that there exists an $\varepsilon > 0$ and infinite set of iteration numbers, $\mathcal{T}$, such that ${G\left({x}_{M}^{k}\right)}^{T}\left({x}_{S}^{k}-{x}_{M}^{k}\right) < -\varepsilon$ for all $k\in \mathcal{T}$. For any $k$ and $r$ with $r > k$, ${x}_{S}^{k}$ is one of the proposals available to ${Master-VI}^{r}\left(G, {\mathrm{\Lambda }}^{r}\right)$; thus, I may use the complementarity conditions in Theorem 1 to derive an inequality. I do this by examining the dual feasibility constraint in ${Master-VI}^{r}\left(G, {\mathrm{\Lambda }}^{r}\right)$ associated with the primal variable ${\lambda }_{k}^{r}$ which is the weight associate with the proposal ${x}_{S}^{k}$: ${{x}_{S}^{k}}^{T}G\left({x}_{M}^{r}\right)+{\theta }^{r}\ge 0$. I may eliminate the variable ${\theta }^{r}$ using the complementarity slackness condition ${\sum }_{i = 1}^{r}{\left({{x}_{S}^{i}}^{T}G\left({x}_{M}^{r}\right)+{\theta }^{r}\right)}^{T}{\lambda }_{i}^{r} = 0$, and using the constraint ${\sum }_{i = 1}^{r}{\lambda }_{i}^{r} = 1$, and the fact that ${x}_{M}^{r} = {\sum }_{i = 1}^{r}{\lambda }_{i}^{r}{x}_{S}^{i}$: ${\theta }^{r} = -{{x}_{M}^{r}}^{T}G\left({x}_{M}^{r}\right)$. This allow us to rewrite the constraint associated with ${\lambda }_{k}^{r}$:${G\left({x}_{M}^{r}\right)}^{T}\left({x}_{S}^{k+1}-{x}_{M}^{r}\right)\ge 0$, for all $k$ and $r$ with $r > k$. Subtracting this from the strict inequality derived earlier, yields ${G\left({x}_{M}^{k}\right)}^{T}\left({x}_{S}^{k}-{x}_{M}^{k}\right)-{G\left({x}_{M}^{r}\right)}^{T}\left({x}_{S}^{k}-{x}_{M}^{r}\right) < -\varepsilon$, for all $k, r\in \mathcal{T}$ with $r > k$. Note that continuity of $G$ makes the left of the inequality continuous in (${x}_{M}, {x}_{S})$. By the property that any infinite subsequence ${\left\{\left({x}_{M}^{k}, {x}_{S}^{k}\right)\right\}}_{k = 1}^{\infty }$ has at least one limit point, there exists a subset $\hat{\mathcal{T}}\subset \mathcal{T}$ such that ${\mathrm{lim}}_{k\to \infty, k\in \hat{\mathcal{T}}}\left({x}_{M}^{k}, {x}_{S}^{k}\right) = \left({\widehat{x}}_{M}, {\widehat{x}}_{S}\right)$, a limit point. Finally, I let $r\to \infty$ through values $r\in \hat{\mathcal{T}}$, in the inequality, and then let $k\to \infty$ through values $k\in \hat{\mathcal{T}}$ (this order of limits ensures that $r > k$ throughout the limiting process) to derive the contradiction $0 = {G\left({\widehat{x}}_{M}\right)}^{T}\left({\widehat{x}}_{S}-{\widehat{x}}_{M}\right)-{G\left({\widehat{x}}_{M}\right)}^{T}\left({\widehat{x}}_{S}-{\widehat{x}}_{M}\right) < -\varepsilon < 0$. That is, ${CG}^{k} < 0$, contradicting the assumption that ${CG}^{k}\ge 0$.

After deriving a SD-VI-NCG, I derive SD-VI with multiple NCG (SD-VI-mNCG) in the next subsection.

3.2.   Simplicial decomposition of VI with multiple NCG (SD-VI-mNCG)

For having multiple NCGs, I introduce a parameter, ${\alpha }_{j}$, to the additional component of SD-VI subproblem approximation. Then, I define the ${j}^{th}$ subproblem approximation mapping:

The corresponding subproblem of SD-VI-mNCG:

Assuming that I have $J$ subproblems, then, for each SD iteration, I have $J$ subproblem proposals, $\left({x}_{S1}^{k}, \dots, {x}_{Sj}^{k}, \dots, {x}_{SJ}^{k}\right)$. For the ${j}^{th}$ subproblem, the proposal, ${x}_{Sj}^{k}$, is stored in the matrix ${X}_{Sj}^{k} = \left[{X}_{Sj}^{k-1}, {x}_{Sj}^{k}\right] = \left[{x}_{Sj}^{1}, {x}_{Sj}^{2}, \dots, {x}_{Sj}^{k}\right]$. Then, let ${X}_{mNCG}^{k} = \left[{X}_{S1}^{k}, \dots, {X}_{Sj}^{k}, \dots, {X}_{SJ}^{k}\right]$.

The feasible set for the master problem at iteration $k$ is restricted to all convex combinations of the $k$ proposals (solutions of all subproblems). Let ${X}_{mNCG}^{k}$ be the $n\times \left(k\times j\right)$ matrix whose columns are solutions ${x}_{Sj}^{i}$ of all subproblems solved at iterations $i = 1, \dots, k$. The weights on the proposals in the convex combination are contained in the vector ${\lambda }_{j}^{k}\in {R}^{k\times j}$. The feasible set for the master problem is defined as ${\mathrm{\Lambda }}_{mNCG}^{k} = \left\{{\lambda }_{j}^{k}\in {R}^{k\times j}|{\left({e}^{k\times j}\right)}^{T}{\lambda }_{j}^{k} = 1, {\lambda }_{j}^{k}\ge 0\right\}$, where ${e}^{k\times j}\in {R}^{k\times j}$ is a vector with 1 for every entry. The master problem at iteration $k$ is defined as

Alternatively, let $conv\left({X}^{k}\right)$ represent the set of all convex combinations of the columns of ${X}^{k}$. Then, I have the following alternative master problem at iteration k, defined as

Since I have $J$ subproblems, I can have $J$ convergence gaps for the stopping condition, according to Theorems 3 and 4, ${CG}_{j}^{k}{ = G\left({x}_{M}^{k}\right)}^{T}\left({x}_{Sj}^{k}-{x}_{M}^{k}\right)$ for all $j$.

We first state the algorithm of SD-VI-mNCG and then the corresponding new theorems about the convergence and the stopping conditions.

[SD-VI-mNCG]

Step 0: Set $k = 0$. Choose $\varepsilon > 0$, ${\alpha }_{j} >$ ${\alpha }_{j-1} > \dots > {\alpha }_{1} > 0, {x}_{M}^{1}\in K$, and ${X}_{i}^{0}$is a null matrix i = 1, $\dots, j$.

Step 1: Increment $k\leftarrow k+1$.

Solve ${\widetilde{S u b}-NLP}_{j}^{k}\left({\tilde{G}}_{j}\left(x; {x}_{M}^{k}; {\alpha }_{j}\right), K\right)$ for all $j$, and place the solution ${x}_{Si}^{k}$ in the matrix $\left[{X}_{Si}^{k-1}, {x}_{Si}^{k}\right], i = 1, \dots, j;$ and update ${X}_{mNCG}^{k} = \left[{X}_{S1}^{k}, \dots, {X}_{Sj}^{k}\right]$. If $k = 1$ then go to Step 2; else if ${CG}_{j}^{k}{ = G\left({x}_{M}^{k}\right)}^{T}\left({x}_{Sj}^{k}-{x}_{M}^{k}\right)\ge -\varepsilon$ then STOP; else go to Step 2.

Step 2: Solve ${\mathrm{M}\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}-\mathrm{V}\mathrm{I}}_{mNCG}^{k}\left(G, {K}_{mNCG}^{k}\right).$ Record $G\left({x}_{M}^{k}\right)$. Go to Step 1.

3.2.1.   The computational sequence of SD-VI-mNCG

The SD-VI-mNCG algorithm first solves ${\widetilde{S u b}-NLP}_{1}^{1}$ with ${\alpha }_{1} = 0$, to generate a column proposal $\left({x}_{Sj}^{1}\right)$ for all $j$ to be placed in the matrix ${X}_{mNCG}^{1}$ of ${Master-VI}_{mNCG}^{1}$ to enlarge the set ${\mathrm{\Lambda }}_{mNCG}^{1}$. ${Master-VI}_{mNCG}^{1}$ is solved to obtain a new ${x}_{M}^{1}$ for new $G\left({x}_{M}^{1}\right)$ and $Q\left({x}_{M}^{1}\right)$. The algorithm proceeds with $k = k+1$, and the new column proposals from ${\widetilde{S u b}-NLP}_{j}^{k}$ is added to the matrix ${X}_{Sj}^{k}$ for all $j$ and update ${X}_{mNCG}^{k}$ which ${Master-VI}_{mNCG}^{k}$ uses to define its convex combinations, and the next solution of the ${Master-VI}^{k}$ produces a new $G\left({x}_{M}^{k}\right)$ and $Q\left({x}_{M}^{k}\right)$ for ${\widetilde{S u b}-NLP}_{j}^{k}$. The algorithm proceeds in this manner until the stopping criterion is satisfied (see Theorem 5 in the following subsection).

3.2.2.   Convergence of SD-VI-mNCG

Theorem 1b. ${\lambda }_{j}^{k}$ solves ${\mathrm{M}\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}-\mathrm{V}\mathrm{I}}_{mNCG}^{k}\left(G, {\mathrm{\Lambda }}_{mNCG}^{k}\right)$ iff there exist ${\lambda }_{j}^{k}\in {R}_{+}^{k\times j}$ and ${\theta }_{j}^{k}\in R$ such that all of the following relations hold:

Theorem 2b. Given the property ${\tilde{G}}_{j}\left({x}_{M}^{k}; {x}_{M}^{k}; {\alpha }_{j}\right) = G\left({x}_{M}^{k}\right)$. If ${x}_{M}^{k}$ solves any one of $J$ subproblem: ${\widetilde{S u b}-NLP}_{j}^{k}\left({\tilde{G}}_{j}\left(x; {x}_{M}^{k}; {\alpha }_{j}\right), K\right)$, then ${x}_{M}^{k}$ solves $VI\left(G, K\right)$.

Proof. Ref. new Theorem 4 in Chung and Fuller ^[20].

Suppose that ${x}_{M}^{k}$ solves the ${j}^{th}$ subproblem, ${\widetilde{S u b}-NLP}_{j}^{k}$. It follows that ${\tilde{G}\left({x}_{M}^{k}; {x}_{M}^{k}; {\alpha }_{j}\right)}^{T}\left(x-{x}_{M}^{k}\right)\ge 0 \;\forall x\in K$. As ${\tilde{G}}_{j}\left({x}_{M}^{k}, {x}_{M}^{k}; {\alpha }_{j}\right) = G\left({x}_{M}^{k}\right)$, I have ${G\left({x}_{M}^{k}\right)}^{T}\left(x-{x}_{M}^{k}\right)\ge 0 \;\forall x\in K$ and I may conclude that ${x}_{M}^{k}$ solves $VI\left(G, K\right)$.

Theorem 3b. Assume that ${\tilde{G}}_{j}\left(x; {x}_{M}^{k}; {\alpha }_{j}\right)$ is strictly monotone in $x$. Given the property ${\tilde{G}}_{j}\left({x}_{M}^{k}; {x}_{M}^{k}; {\alpha }_{j}\right) = G\left({x}_{M}^{k}\right)$. If any ${CG}_{j}^{k}{ = G\left({x}_{M}^{k}\right)}^{T}\left({x}_{Sj}^{k}-{x}_{M}^{k}\right)\ge 0$, then ${x}_{M}^{k}$ solves $VI\left(G, K\right)$.

Proof. Ref. new Theorem 6(a) of Chung and Fuller ^[20].

We shall show that if ${x}_{M}^{k}$ does not solve $VI\left(G, K\right)$, then ${CG}^{k} < 0$. By Theorem 2b, ${x}_{M}^{k}\ne {x}_{Sj}^{k}$, and since ${\tilde{G}}_{j}\left({x}_{M}^{k}; {x}_{M}^{k}; {\alpha }_{j}\right) = G\left({x}_{M}^{k}\right)$, strict monotonicity of ${\tilde{G}}_{j}\left({x}_{M}^{k}; {x}_{M}^{k}; {\alpha }_{j}\right)$ implies that ${\left(G\left({x}_{M}^{k}\right)-{\tilde{G}}_{j}\left({x}_{M}^{k}; {x}_{M}^{k}; {\alpha }_{j}\right)\right)}^{T}\left({x}_{M}^{k}-{x}_{Sj}^{k+1}\right) > 0$. Since ${x}_{Sj}^{k}$ solves ${\widetilde{S u b}-NLP}_{j}^{k}\left({\tilde{G}}_{j}\left(x; {x}_{M}^{k}; {\alpha }_{j}\right), K\right)$, it follows that ${{\tilde{G}}_{j}\left({x}_{Sj}^{k}; {x}_{M}^{k}; {\alpha }_{j}\right)}^{T}\left({x}_{M}^{k}-{x}_{Sj}^{k}\right)\ge 0$. Adding this last inequality to the strict inequality and multiplying by -1 yields ${G\left({x}_{M}^{k}\right)}^{T}\left({x}_{Sj}^{k}-{x}_{M}^{k}\right) = {CG}_{j}^{k} < 0$.

Theorem 5. Let ${\alpha }_{j} > {\alpha }_{j-1} > \dots > {\alpha }_{1} > 0$ and let ${z}_{j}^{k}$ be the optimal value of objective function of the subproblem ${\widetilde{S u b}-NLP}_{j}^{k}$. Then ${z}_{j}^{k}\ge {z}_{j-1}^{k}\ge \dots \ge {z}_{1}^{k}$, and ${CG}_{j}^{k}\ge {CG}_{j-1}^{k}\ge \dots \ge {CG}_{1}^{k}$, where ${CG}_{j}^{k}{ = G\left({x}_{M}^{k}\right)}^{T}\left({x}_{Sj}^{k}-{x}_{M}^{k}\right)$.

Proof. For any subproblem $i$, let ${x}_{Si}^{k}$ solve ${\widetilde{S u b}-NLP}_{i}^{k}\left({\tilde{G}}_{j}\left(x; {x}_{M}^{k}; {\alpha }_{i}\right), K\right)$, and let ${z}_{i}^{k}\left({x}_{Si}^{k}\right)$ be the corresponding objective value. That is, ${z}_{i}^{k} = {G\left({x}_{M}^{k}\right)}^{T}{x}_{Si}^{k}+\frac{{\alpha }_{i}}{2}{\left({x}_{Si}^{k}-{x}_{M}^{k}\right)}^{T}Q\left({x}_{M}^{k}\right)\left({x}_{Si}^{k}-{x}_{M}^{k}\right).$ As ${x}_{Si+1}^{k}$ is another feasible solution, it follows that ${G\left({x}_{M}^{k}\right)}^{T}{x}_{Si+1}^{k}+\frac{{\alpha }_{i}}{2}{\left({x}_{Si+1}^{k}-{x}_{M}^{k}\right)}^{T}Q\left({x}_{M}^{k}\right)\left({x}_{Si+1}^{k}-{x}_{M}^{k}\right)\ge {G\left({x}_{M}^{k}\right)}^{T}{x}_{Si}^{k}+\frac{{\alpha }_{i}}{2}{\left({x}_{Si}^{k}-{x}_{M}^{k}\right)}^{T}Q\left({x}_{M}^{k}\right)\left({x}_{Si}^{k}-{x}_{M}^{k}\right)$. Since ${\alpha }_{i+1} > {\alpha }_{i}$ for all $j$, I have ${G\left({x}_{M}^{k}\right)}^{T}{x}_{Si+1}^{k}+\frac{{\alpha }_{i+1}}{2}{\left({x}_{Si+1}^{k}-{x}_{M}^{k}\right)}^{T}Q\left({x}_{M}^{k}\right)\left({x}_{Si+1}^{k}-{x}_{M}^{k}\right)\ge {G\left({x}_{M}^{k}\right)}^{T}{x}_{Si}^{k}+\frac{{\alpha }_{i}}{2}{\left({x}_{Si}^{k}-{x}_{M}^{k}\right)}^{T}Q\left({x}_{M}^{k}\right)\left({x}_{Si}^{k}-{x}_{M}^{k}\right)$. That is ${z}_{i+1}^{k}\ge {z}_{i}^{k}$. Adding ${(-G\left({x}_{M}^{k}\right)}^{T}\left({x}_{M}^{k}\right))$ in both sides, I have ${G\left({x}_{M}^{k}\right)}^{T}\left({x}_{Si+1}^{k}-{x}_{M}^{k}\right)+\frac{{\alpha }_{i+1}}{2}{\left({x}_{Si+1}^{k}-{x}_{M}^{k}\right)}^{T}Q\left({x}_{M}^{k}\right)\left({x}_{Si+1}^{k}-{x}_{M}^{k}\right)\ge {G\left({x}_{M}^{k}\right)}^{T}\left({x}_{Si}^{k}-{x}_{M}^{k}\right)+\frac{{\alpha }_{i}}{2}{\left({x}_{Si}^{k}-{x}_{M}^{k}\right)}^{T}Q\left({x}_{M}^{k}\right)\left({x}_{Si}^{k}-{x}_{M}^{k}\right)$. It implies ${CG}_{i+1}^{k}+\frac{{\alpha }_{i+1}}{2}{\left({x}_{Si+1}^{k}-{x}_{M}^{k}\right)}^{T}Q\left({x}_{M}^{k}\right)\left({x}_{Si+1}^{k}-{x}_{M}^{k}\right)\ge {CG}_{i}^{k}+\frac{{\alpha }_{i}}{2}{\left({x}_{Si}^{k}-{x}_{M}^{k}\right)}^{T}Q\left({x}_{M}^{k}\right)\left({x}_{Si}^{k}-{x}_{M}^{k}\right)$. When ${x}_{M}^{k}$ solve ${\widetilde{S u b}-NLP}_{i+1}^{k}$, ${x}_{M}^{k} = {x}_{Si+1}^{k}$ and $\frac{{\alpha }_{i+1}}{2}{\left({x}_{Si+1}^{k}-{x}_{M}^{k}\right)}^{T}Q\left({x}_{M}^{k}\right)\left({x}_{Si+1}^{k}-{x}_{M}^{k}\right) = 0.$Then, ${CG}_{i+1}^{k}\ge {CG}_{i}^{k}+\frac{{\alpha }_{i}}{2}{\left({x}_{Si}^{k}-{x}_{M}^{k}\right)}^{T}Q\left({x}_{M}^{k}\right)\left({x}_{Si}^{k}-{x}_{M}^{k}\right)\ge {CG}_{i}^{k}$.

With Theorem 4, the stopping condition can be defined by ${CG}_{j}^{k}{ = G\left({x}_{M}^{k}\right)}^{T}\left({x}_{Sj}^{k}-{x}_{M}^{k}\right)$. In addition, it should be noted that the SD-VI-mNCG can be easily reduced to SD-NLP-mNCG for nonlinear programming by replacing the ${\mathrm{M}\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}-\mathrm{V}\mathrm{I}}_{mNCG}^{k}\left(G, {K}_{mNCG}^{k}\right)$ with ${\mathrm{M}\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}-\mathrm{N}\mathrm{L}\mathrm{P}}^{k}\left(G, {K}_{mNCG}^{k}\right)$ in the Subsection 2.3.

4.   Empirical results of SD with multiple NCG

In my implementation, the test models, the column generation algorithms, and the reference algorithm without column generation are coded into GAMS programs, executed on a PC with Intel Processor (12^th Gen Core i7-12700, 2100 Mhz, 12 Core(s), 20 Logical Processors and 16 GB RAM. The reference algorithm is a relaxation algorithm ^[29] by which the test models are solved so that I can have reference results for evaluating the accuracy and speed of the column generation algorithms. The master problem is also solved by the same relaxation algorithm. Each iteration of the relaxation algorithm is a nonlinear programming (NLP) calculation. For all models solved by SD-VI method, in the first subproblem at step $k$ = 1, the initial values of $G\left({x}_{M}^{1}\right)$ and ${\alpha }_{j}$ are arbitrarily set to 10 and 0, respectively, which produces first proposals that ignore the nonlinear column generation component. In the relaxation method for the reference calculations, the first NLP also has the initial values of all ${x}_{M}$ are set to 10.

For all column generation calculations, I set the convergence tolerance $\varepsilon = 0.000001$, except where noted. The relaxation iterations of the restricted master problems and the reference method, and all subproblems, except subproblems with ${\alpha }_{j} = 0,$ are solved by CONOPT 3 called from GAMS. The subproblems with ${\alpha }_{j} = 0$ are linear programs that are solved by CPLEX. It should be noted that the SD-VI-mNCG with ${\alpha }_{j} = 0$ reduces to the SD-VI method of ^[26]. The SD-VI-mNCG with one NCG and ${\alpha }_{j} = 0.5$ reduces to the SD-VI-NCG method of ^[28].

In the calculations, many similar optimization problems are solved repeatedly: The subproblem, a sequence of optimization problems related to the restricted master problem, and another related to the reference algorithm; except for the first time, each time an optimization calculation is done, it is started from the last solution in order to reduce computation time. This process is not used in the multiple column generation method since I assume that each nonlinear column generation subproblem is implemented in a server individually and all these subproblems are solved in parallel.

To summarize, even if a column generation algorithm takes more time than no column generation algorithm (the relaxation algorithm), there can be an advantage to a column generation approach, in model development. However, if column generation takes a very large amount of time, then such an advantage might not be worthwhile. Therefore, I am interested in tests that measure the time of column generation, compared with no column generation, and in variants of multiple column generation algorithms that may run in shorter times.

4.1.   A simple example (NEW)

An example is Example 2 of Nagurney and Dhanda ^[30], a VI problem of a multiproduct multipollutant oligopolistic market model (MMOM-VI) with ambient-based pollution permits and transaction costs. In this model, $m$ firms or sources of pollution, $r$ pollutants emitted by the firms, and $n$ receptor points exist. I also let ${e}_{i}^{t}$ be the amount of pollutant $t$ emitted by firm $i$, ${l}_{ij}^{t}$ be the number of licenses for pollutant $t$ at receptor point $j$ held by firm $i$, ${l}_{ij}^{t0}$ be the initial allocation of licenses made by a regulatory agency and ${p}_{j}^{t}$ be the price of the licenses for pollutant $t$ that affects receptor point $j$.

The cost of purchasing licenses for specific pollutant $t$ that affects receptor point $j$ for source $i$ is given by $\sum _{j = 1}^{n}{p}_{j}^{t*}\left({l}_{ij}^{t}-{l}_{ij}^{t0}\right)$, where ${p}_{j}^{t*}$ is the market clearing price determined by the VI.

Each firm $i$ in the oligopoly is faced with cost ${f}_{i}$ for producing the vector of quantities ${q}_{i}$, where

and ${q}_{id}$ = the quantity of product $d$ produced by firm $i$.

Each firm i in the region is also faced with joint-cost ${g}_{i}$, where ${g}_{i} = {g}_{i}\left({e}_{i}, {q}_{i}\right) = \sum _{t = 1}^{r}\left[{g1}_{it}{\left({e}_{i}^{t}\right)}^{2}+{g2}_{it}{e}_{i}^{t}+{g4}_{it}\right]+\sum _{d = 1}^{s}{g3}_{it}{q}_{id}$.

The transaction cost function employed by firm $i$ for pollutant t at receptor point $j$ was of the form ${c}_{ij}^{t} = {c}_{ij}^{t}\left({l}_{ij}^{t}\right) = {\varnothing 1}_{ijt}{\left({l}_{ij}^{t}\right)}^{2}+{\varnothing 2}_{ijt}{l}_{ij}^{t}+{\alpha }_{ijt}$.

The firms are oligopolistic in their product markets, and they affect the prices of the outputs. The price of product $d$ is denoted by ${\rho }_{d} = {\rho }_{d}\left(\sum _{i = 1}^{m}{q}_{id}\right)$.

Then, I can have the MMOM-VI model shown below.

MMOM-VI: Find $({q}_{id}^{\mathrm{*}}, {e}_{i}^{\mathrm{t}\mathrm{*}}, {l}_{ij}^{\mathrm{t}\mathrm{*}})\in {K}_{q, e, l}$ such that

where $\sum _{i = 1}^{m}{l0}_{ij}^{t}$ is the environmental quality standard, and ${h}_{ij}^{t}$ denotes the contribution of one unit of emission by source $i$ to the average pollutant concentration of type $t$ at receptor point $j$. All the parameters and the reference results can be found in the online appendix of ^[30]. For convenience, they are given in the Appendix of the current paper.

Similarly, I adopt $Q\left({x}_{M}\right) = \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\mathit{\nabla} G\left({x}_{M}\right)$. That is, let ${x}_{M} = ({q}_{idM}, {e}_{iM}^{t}, {l}_{ijM}^{t})$, and I have

which is added to the subproblem and the $\boldsymbol{K} = {K}_{q, e, l}$.

Table 1 summarises the performance of SD-VI-NCG(x) and SD-VI-mNCG of MMOM-VI, in which the first column provides the names of the models. Ref_method is the reference method, which is the relaxation algorithm. SD-VI-NCG($\mathrm{\alpha }/2$) is the method of simplicial decomposition with one NCG of $\mathrm{\alpha }/2$. For instance, SD-VI-NCG(0.5) is the SD-VI-NCG with $\mathrm{\alpha }/2 = 0.5$. The last two rows of Table 1 present the performance of SD-VI-mNCG. I use m = 3, 5 to represent how many NCG subproblems are used. For example, SD-VI-3NCG(0.1, 0.3, 0.5) consists of three subproblems with $\mathrm{\alpha }/2$ = 0.1, 0.3, and 0.5. It is interesting that all two SD-VI-mNCG methods use the smallest number of decomposition steps, 7. However, there is no conclusive result in terms of computational performance. The last column of Table 1 shows the stopping condition mentioned in Theorem 5 that the SD-VI was terminated by the subproblem with the largest value of $\mathrm{\alpha }/2$.

4.2.   The test model (was 4.1)

Since the SD method is usually used for solving transportation network equilibrium problems, I use this kind of problems to develop test model. A large-scale real transportation network equilibrium problem is used¹. The network of this problem represents the extra-urban area of the city of Arezzo (Italy). It consists of 213 nodes, 598 arcs and 2423 O/D pairs. The form of the link $\left(i, j\right)$ cost functions is asymmetric:

¹ https://pages.di.unipi.it/passacantando/test_networks.html.

where

The objective function of the subproblem, ${\widetilde{S u b}-NLP}_{j}^{k}\left({\tilde{G}}_{j}\left(x; {x}_{M}^{k-1};{\alpha }_{j}\right), K\right)$, is

where I adopt $Q\left({x}_{M}^{k}\right) = \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\mathit{\nabla} G\left({x}_{M}^{k}\right)$. That is,

with different ${\alpha }_{j}$, I can have different nonlinear column generation subproblems.

We first provide the computational performance results of SD-NLP with multiple NCG, and then the results of SD-VI with multiple NCG. For having SD-NLP test model, I set the asymmetric factor d = 0.

4.3.   Performance of SD with multiple NCG for NLP (was 4.2)

In this subsection, I report the solution time of simplicial decomposition with one NCG subproblem of NLP. Table 2 summarises the computational performance of SD-NLP-NCG, in which the first column provides the names of the models. SD-NLP-NCG($\mathrm{\alpha }/2$) is the method of simplicial decomposition with one NCG of $\mathrm{\alpha }/2$. For instance, SD-NLP-NCG(0.5) is the SD-NLP-NCG with $\mathrm{\alpha }/2 = 0.5$. Hence, SD-NLP-NCG(0) is the simple SD-NLP with LP subproblem, and SD-NLP-NSD(0.5) is the one in Larsson et al. ^[28]. From Table 2, the empirical result shows that the decomposition steps taken by SD-NLP is 87, which is expected and much higher than any SD-NLP-NCG methods with $\mathrm{\alpha }/2 = 0.1, \mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }0.2, \mathrm{ }0.3, \mathrm{ }0.4, \mathrm{ }0.5, \mathrm{ }\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{ }1.0.$ The range of the decomposition steps is 8 to 22. The SD-NLP-NCG with $\mathrm{\alpha }/2$ = 0.4 provides the best performance in terms of the number of decomposition steps, 8. However, there is no clear relationship between the number of decomposition steps and the total solution time while changing the value of $\mathrm{\alpha }$. On the other hand, the fastest method is SD-NLP-NCG(0) because its subproblem is an LP solved by CPLEX, which is one of the fastest LP solvers. It is noted that SD-NLP-NCG(0) is a simple SD-NLP, and it takes more time to solve its master NLP problem comparing with other methods with different values of $\mathrm{\alpha }$. It implies that SD-NLP-NCG(0) may not be the fastest if a larger model incurs a larger number of generated columns in the master problem.

The last three rows of Table 2 present the performance of SD-NLP-mNCG. I use m = 3, 5, and 10 to represent how many NCG subproblems are used. For example, SD-NLP-3NCG(0.1, 0.3, 0.5) consists of three subproblems with $\mathrm{\alpha }/2$ = 0.1, 0.3, and 0.5. All three SD-NLP-mNCG methods use fewer decomposition steps, as expected, than SD-NLP-NCG. It is interesting that all three SD-NLP-mNCG methods use 6 decomposition steps, which may imply that computational performance of SD-NLP-mNCG does not rely on the number of subproblems.

4.4.   Performance of SD with multiple NCG for VI (was 4.3)

In this subsection, I report the solution time of simplicial decomposition of VI with one NCG subproblem, which is similar to the one in Larsson et al. ^[28]. Table 3 summarizes the computational performance of SD-VI-NCG(x). This empirical result shows that the SD-VI-NSD with $\mathrm{\alpha }/2$ = 0.3 provides the best performance in terms of the number of decomposition steps, 13. It is expected that the number of decomposition steps of SD-VI-NSD($\mathrm{\alpha }/2)$is much smaller than that of SD-VI, 152, as discussed in Larsson et al. ^[28]. Like the SD-NLP-NCG results, there is no clear relationship between the number of decomposition steps and the total solution time while changing the value of $\mathrm{\alpha }$. On the other hand, the fastest method is SD-VI because its subproblem is an LP solved by CPLEX. It is noted that SD-VI takes more time to solve its master VI problem. It implies that SD-VI may not be the fastest if a larger model incurs a larger number of generated columns in the master problem.

Table 4 reports the solution time of simplicial decomposition with multiple NCG subproblems (SD-VI-mNCG). That is, in each iteration of the SD-VI, there are a number of NCG subproblems with different values of $\mathrm{\alpha }/2$. I first use three NCG subproblems for my test model. Table 4 summarises the computational performance of SD-VI-3NCG with nine sets of $\mathrm{\alpha }/2$. For example, SD-VI-3NCG (0.1, 0.2, 0.3) consists of three subproblems with $\mathrm{\alpha }/2$ = 0.1, 0.2, and 0.3. As expected, on average, SD-VI-3NCG uses fewer decomposition steps than SD-VI-NSD. The range of decomposition steps of SD-VI-3NCG is 11 to 13, while SD-VI-NCG is 13 to 23 shown in Table 3. Moreover, it is observed that the number of decomposition steps exploited in SD-VI-3NCG is more stable. I further use five NCG subproblems, with $\mathrm{\alpha }/2$ = 0.1, 0.2, 0.3, 0.4, and 0.5, (SD-VI-5NCG(0.1, …, 0.5)), the resulting decomposition steps is 13, in the range of 11–13. 11 decomposition steps are also used in ten NCG subproblems with $\mathrm{\alpha }/2$ = 0.1, 0.2, …, and 1.0, (SD-VI-5NCG(0.1, …, 0.5)).

Remarks on using the CONOPT solver: Solvers, the MINOS and the PATHNLP cannot find the solution after 8 hours. The CONOPT solver can obtain solutions within an hour in all cases. When the CONOPT solver generates a "feasible" solution to the subproblem, not an optimal solution, in particular, in the first few iterations, the corresponding CG cannot be used for the stopping criterion.

In implementing the multiple subproblem cases, there are two phases. The first phase is to allow a few subproblems to generate infeasible, feasible, or optimal solutions due to the capacity of CONOPT. When any subproblem generates an infeasible solution in an iteration, I use its previous solution obtained from the previous iteration to substitute this infeasible solution. By doing so, a set of feasible solutions can be obtained in each iteration until all subproblems obtain their optimal solutions and thereafter, which is the second phase and follows my theoretical results.

4.5.   Re-optimization performance (was 4.4)

For studying the re-optimization performance of the SD-VI-mNCG, the asymmetric parameter d of the test model (see Subsection 4.1) is changed to 0.6 from 0.5. There are three approaches, SD-VI_(0.6), Final_soln(0.5)_(0.6), and All_soln(0.5)_(o.6), to re-optimization analysis of different multi-subproblem methods.

(1) SD-VI_(0.6): Use the exiting coding of SD-VI-NCG(x), where x = 0, 0.1, 0.3, 0.5, and 1.0, with new d = 0.6, as reference result.

(2) Final_soln(0.5)_0.6: Use the previous final solution obtained from d = 0.5 as the feasible starting point to run the re-optimization process for the model with new d = 0.6.

(3) All_soln(0.5)_(0.6): Use all previous proposals with d = 5 to start the re-optimization process with new d = 0.6.

Table 5 summarizes the computational results. As expected, Approach 3, All_soln(0.5)_(0.6), used the smallest decomposition steps in all methods. With Approach 3, the decomposition steps of SD-VI-5NSD is 8, which is smaller than the decomposition steps of SD-VI-NSD(x), ranging from 10 to 23.

4.6.   A comparison summary of the method SD-VI-mNCG, SD-VI, and SD-VI-NCG

Table 6 summarized the performance of the method SD-VI-mNCG, the SD-VI of Lawphongpanich and Hearn ^[26], and the SD-VI-NCG(0.5) of Larsson et al. ^[28]. In terms of column generation steps, SD-VI-mNCG is always smaller than SD-VI and SD-NCG. In terms of computational performance, for the model MMOM, the methods with NCG are faster than the SD-VI. However, SD-VI is the fastest method for the model Arezzo. As mentioned above, SD-VI takes more time to solve its master VI problem. It implies that SD-VI may not be the fastest if a larger model incurs a larger number of generated columns in the master problem. In short, solving multiple subproblems in each SD step can reduce the number of iterations. There is no conclusive result for the computational performance. Similar performance results are found in the re-optimization processes.

4.7.   Other practical examples

The simple example of a multiproduct multipollutant oligopolistic market model in Subsection 4.1 and the application of transportation network equilibrium problems in Subsection 4.2 are practical examples of directly applying the SD-VI-mNCG to the problem. Indeed, the SD-VI-mNCG can be used in some combined methods that solve part of the problems. For example, Sharma et al. ^[31] developed a method combining the simplicial decomposition, gauss-seidel, and the augmented Lagrangian to solve the coordination problem of the optimal power and traffic flows with EVs. Wang et al. ^[32] employed SD in their combined method to have a parallel decentralized solution for multi-regional unit commitment with convex AC power flow constraints. It should be noted that the constraint sets of SD-VI-mNCG can be convex constraint sets. In addition, the SD-VI-mNCG can also be employed to solve stochastic transportation problems ^[28].

5.   Conclusions and further research

In this paper, I derive a simplicial decomposition with multiple nonlinear column generation subproblems, which can be applied to nonlinear programming and variational inequalities. I provide convergence properties and derive new stopping condition.

From the computational performance, I generally conclude that for solution methods, simplicial decomposition of variational inequalities with multiple nonlinear column generation subproblems (SD-VI-mNCG), can perform better than the simplicial decomposition of variational inequalities (SD-VI) and the SD-VI with a nonlinear column generation subproblem (SD-VI-NCG) in terms of decomposition steps. The number of column generation steps is greatly reduced, and the long-tail convergence property of SD-VI is alleviated. The same conclusion can be applied to nonlinear programming.

Concerning the selection of the number of subproblems, m, the empirical results show that 3 to 5 is fine.

Concerning the selection of the values ${\mathrm{\alpha }}_{\mathrm{j}}$, I use $\left({\mathrm{\alpha }}_{\mathrm{j}}/2\right)$ = 1/2 as the reference value by which the subproblem becomes diagonalized Newton algorithm. This approach may be beneficial in other column generation schemes, like Dantzig-Wolfe decomposition of VI, which is one of further research directions. Moreover, employing SD-VI-mNCG for different kind of VI problems, like large-scale energy equilibrium problems, can be another research direction.

On the other hand, since there is no conclusive result for the computational performance, it is worthwhile to consider the computational performance as one of the further research topics. For instance, I may not need to wait for all subproblems to be solved to move on to the next SD step. The calculation framework can be one master problem to m subproblems until one of the subproblems solves the VI problem.

Use of AI tools declaration

The author declares that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Financial support for William Chung's work came from the Research Grants Council of Hong Kong S.A.R., China (CityU 11500022).

Conflict of interest

The author declares no conflicts of interest.

Appendix

Parameters and reference results of MMOM-VI.

For cost function, ${f}_{i} = {f}_{i}\left({q}_{i}\right) = \sum _{d = 1}^{s}\left[{c}_{id}{q}_{id}+\frac{{\beta }_{id}}{{\beta }_{id}+1}{K}_{id}^{-1/{\beta }_{id}}{q}_{id}^{({\beta }_{id}+1)/{\beta }_{id}}\right]$:

For joint-cost function, ${g}_{i} = {g}_{i}\left({e}_{i}, {q}_{i}\right) = \sum _{t = 1}^{r}\left[{g1}_{it}{\left({e}_{i}^{t}\right)}^{2}+{g2}_{it}{e}_{i}^{t}+{g4}_{it}\right]+\sum _{d = 1}^{s}{g3}_{it}{q}_{id}$:

For transaction cost function, ${c}_{ij}^{t} = {c}_{ij}^{t}\left({l}_{ij}^{t}\right) = {\varnothing 1}_{ijt}{\left({l}_{ij}^{t}\right)}^{2}+{\varnothing 2}_{ijt}{l}_{ij}^{t}+{\alpha }_{ijt}$:

For the price functions of product $d$: ${\rho }_{d} = {\rho }_{d}\left(\sum _{i = 1}^{m}{q}_{id}\right)$, ${\rho }_{d} = 5000$.

For constraints ${{l}_{ij}^{t}-{h}_{ij}^{t}e}_{i}^{t}\ge 0$:

For $\sum _{i = 1}^{m}{l0}_{ij}^{t}-\sum _{i = 1}^{m}{l}_{ij}^{t}\ge 0$, ${l0}_{ij}^{t} = 3$.

Reference results are generated by PIES without a decomposition algorithm: