Solvability and stability analysis of a coupled system involving generalized fractional derivatives

1.   Introduction

We mainly consider the linear multiplicative problem of the following form:

where $c_{i}$ is an $n$-dimensional column vector, $d_{i}$ is a real number and $\alpha_{i}$ is a real number different from zero. $T$ represents the transpose of a vector, $A\in \mathbb{R}^{m\times n}$, $b\in \mathbb{R}^{m}$ and $\mathcal{Y}$ is a non-empty bounded set. In this article, for any $y\in \mathcal{Y}$, we suppose that $c_{i}^{T}y+d_{i} > 0,$ $i = 1, ..., p$.

LMP is a typical non-convex optimization problem with important applications in real life. It and its variants have been applied in various fields such as robust optimization ^[1], financial optimization ^[2], VLSI chip design ^[3], decision tree optimization ^[4], network flow optimization ^[5], supply chain problem ^[6], investment portfolio ^[7], etc. Moreover, LMP is an NP-hard problem ^[8] with multiple local solutions rather than global solutions, which increases the computational complexity. Hence, researching this problem holds great significance.

Various solution algorithms for LMP and its special forms have been proposed by numerous experts and scholars. These algorithms can be broadly categorized into the following groups: branch-and-bound algorithms ^{[9,10,11,12,13,14]}, finite algorithm ^[15], heuristic method ^[16], approximation algorithm ^[17], polynomial time algorithm ^[18], parameterized simplex method ^[19], cutting-plane method ^[20], level set algorithm ^[21], etc. Despite the advancements made by these approaches, solving high-dimensional LMP remains a challenging task. In the past 20 years, searching for global solutions of LMP in the outer space using different relaxation methods has become a hot topic among scholars ^{[22,23,24,25,26,27,28,29]}. For example, the authors in references ^[22,27] used new convex relaxation techniques to put forward different outer space branch-and-bound algorithms for solving LMP. References ^{[23,25,26,28,29,30]} adopted various linear relaxation programming problems and proposed novel outer space branch-and-bound algorithms, respectively.

In this paper, an outer space branch-and-bound algorithm is designed to solve large-scale LMP. The major characteristics of the algorithm are as follows: First, $p$ auxiliary variables are introduced to transform LMP into an equivalent problem, where $p$ is the number of linear functions. Second, based on the properties of exponential and logarithmic functions, the second equivalent problem (ELMP1) is established. Third, a novel linear relaxation approach is employed to derive the linear relaxation programming problem for ELMP1. Moreover, the branching rule in $p$-dimensional outer space is given, and the corresponding outer space branch-and-bound algorithm is developed by embedding the rectangular compression technique into the branch-and-bound framework. Finally, the computational complexity of the algorithm is analyzed to estimate the number of iterations in the worst case, which also implies that our algorithm is convergent.

Compared with existing methods ^{[10,23,26,28]}, the proposed algorithm has the following advantages:

(1) The solved LMP is universal, and the exponents of its objective function are real numbers rather than being limited to just 1 or only positive values.

(2) After the first equivalent conversion, the exponents of the objective function are all positive. Therefore, when linearly approximating the objective function of equivalent problems, there is no need to consider the case where the coefficient is negative, which simplifies the implementation of linear relaxation.

(3) The branching process takes place only in the $p$-dimensional outer space of the linear function. This leads to cost savings compared to the branching operation in the $n$-dimensional decision variable space, as $p$ is often much smaller than $n$ in practical problems.

(4) To demonstrate the efficiency of the algorithm, we compared it with the methods in references ^{[10,23,26,28]}, our algorithm is suitable for solving LMP with large-scale decision variables.

The rest of this paper is organized as follows: Section 2 presents the equivalent problems of LMP and establishes its linear relaxation problem. In Section 3, the branching operation and compression technology are discussed in detail. Moreover, the outer space branch-and-bound algorithm and its properties are described. Section 4 provides numerical comparison results for some problems. Finally, a brief summary of this paper is presented in Section 5.

2.   Equivalent problem and its linear relaxation problem

In order to search the global optimal solution of LMP, we transform it into the equivalent problem (ELMP). For convenience, we first define the following sets:

For any $i\in\{1, ..., p\}$, denote

Since $c_{i}^{T}y+d_{i}$ is a bounded linear function, by solving the above $2p$ linear programs, we can easily get that $\underline{t}_{i}^{0}$ and $\overline{t}_{i}^{0}$. Simultaneously, the initial hyper-rectangle

is constructed. Thus, let us consider the following equivalent problem:

We denote the feasible region of ELMP as $\mathcal{V} = \{t_{i} = c_{i}^{T}y+d_{i}, i\in I^{+}, \; \; t_{i} = \frac{1}{z_{i}}, \; \; z_{i} = c_{i}^{T}y+d_{i}, \; \; i\in I^{-}, \; \; y\in \mathcal{Y}, \; \; t\in \mathcal{H}^{0}\}.$ It is evident that $\mathcal{V}$ is a nonempty bounded compact set if and only if $\mathcal{Y}\neq \emptyset$. Theorem 1 below explains the equivalence between LMP and ELMP.

Theorem 1. $(y^*, t^*, z^*)$ is a global optimal solution for ELMP if and only if $y^*$ is a global optimal solution for LMP, where $t_{i}^* = c_{i}^{T}y^*+d_{i}$ when $i\in I^{+}$, $t_{i}^* = \frac{1}{z_{i}^*}$ and $z_{i}^* = c_{i}^{T}y^*+d_{i}$ when $i\in I^{-}$.

Proof. We will develop our proof in two aspects. On one hand, for any $y\in \mathcal{Y}$, let $t_{i} = c_{i}^{T}y+d_{i}$ for $i\in I^{+}$, $t_{i} = \frac{1}{z_{i}}$ and $z_{i} = c_{i}^{T}y+d_{i}$ for $i\in I^{-}$, thus $(y, t, z)\in \mathcal{V}$. If $y^*$ is a global optimal solution for LMP, then $t_{i}^* = c_{i}^{T}y^*+d_{i}$ for $i\in I^{+}$, $t_{i}^* = \frac{1}{z_{i}^*}$ and $z_{i}^* = c_{i}^{T}y^*+d_{i}$ for $i\in I^{-}$, so $(y^*, t^*, z^*)\in \mathcal{V}$, which shows that $(y^*, t^*, z^*)$ is a feasible solution to ELMP, and by the optimality of $y^*$, the following inequalities hold:

Thus, $(y^*, t^*, z^*)$ is a global optimal solution for ELMP.

On the other hand, if $(y^*, t^*, z^*)$ is a global optimal solution of ELMP, where $t^*, z^*$ are satisfied: if $i\in I^{+}$, then $t_{i}^* = c_{i}^{T}y^*+d_{i}$; if $i\in I^{-}$, then $t_{i}^* = \frac{1}{z_{i}^*}$ and $z_{i}^* = c_{i}^{T}y^*+d_{i}$. Suppose $y$ is a global optimal solution of LMP such that $\prod_{i = 1}^{p}(c_{i}^{T}y+d_{i})^{\alpha_{i}} < \prod_{i = 1}^{p}(c_{i}^{T}y^*+d_{i})^{\alpha_{i}}$ holds, for $i\in I^{+}$, let $t_{i} = c_{i}^{T}y+d_{i}$, for $i\in I^{-}$, let $t_{i} = \frac{1}{z_{i}}$ and $z_{i} = c_{i}^{T}y+d_{i}$, it follows that

This contradicts the optimality of $(y^*, t^*, z^*)$, thus $y^*$ is a global optimal solution of LMP. □

For the sake of convenience, we denote $\beta_{i} = -\alpha_{i} > 0$ for $i\in I^{-}$. In the meantime, ${\rm ln}(\bullet)$ and ${\rm exp}(\bullet)$ represent the logarithmic and the exponential functions, respectively. The objective function of ELMP is further equivalently transformed according to the properties of the exponential and logarithmic functions, namely,

where $\zeta\in \mathbb{R}^{p}$ and $\zeta = [\alpha_{1}, \alpha_{2}, \cdots, \alpha_{\kappa}, \beta_{\kappa+1}, \cdots, \beta_{p}]$. Hence, ELMP is reformulated to the following form:

where $\mathcal{H}^{k}$ represents $\mathcal{H}^{0}$ or the sub-rectangle of $\mathcal{H}^{0}$. Obviously, the optimal solution for ELMP1 is the same as that for ELMP. Therefore, we shift our focus to solving ELMP1, but ELMP1 cannot be solved directly due to the nonlinearity of the objective function and the constraints $z_{i} = \frac{1}{t_{i}}$ for $i\in I^{-}$. To address this, we propose a linear relaxation technique to obtain the lower bound problem of ELMP1.

Theorem 2. For $i = 1, ..., p$, $t_{i}\in [\underline{t}_{i}, \overline{t}_{i}]$, define

where $k_{i} = \frac{{\rm ln}\overline{t}_{i}-{\rm ln}\underline{t}_{i}}{\overline{t}_{i}-\underline{t}_{i}}$. Then we have the following conclusions:

(i) $\underline{g}(t_{i})\leq g(t_{i})$; (ii) $\underline{\mathcal{L}}(y, t, z)\leq \mathcal{L}(y, t, z).$

Proof. Since the function $g(t_{i})$ is a monotonically increasing concave function on $[\underline{t}_{i}, \overline{t}_{i}]$ with respect to $t_{i}$, $\underline{g}(t_{i})$ is its secant approximation, so (ⅰ) and (ⅱ) obviously hold. □

Theorem 3. For each $i\in I^{-}$, define

Then the functions $\underline{\psi}(t_{i})$ and $\overline{\psi}(t_{i})$ satisfy the conclusions below:

(i) $\underline{\psi}(t_{i})\leq \frac{1}{t_{i}}\leq \overline{\psi}(t_{i})$;

(ii) Denote $\Delta t_{i} = \overline{t}_{i}-\underline{t}_{i}$, then $\lim\limits_{\Delta t_{i}\rightarrow 0}(\frac{1}{t_{i}}-\underline{\psi}(t_{i})) = 0$, $\lim\limits_{\Delta t_{i}\rightarrow 0}(\overline{\psi}(t_{i})-\frac{1}{t_{i}}) = 0$.

Proof. (ⅰ) For each $i\in I^{-}$, since $t_{i}\in [\underline{t}_{i}, \overline{t}_{i}]$ and $t_{i} > 0$, it follows from the definitions of $\underline{\psi}(t_{i})$ and $\overline{\psi}(t_{i})$ that

and

(ⅱ) From (ⅰ), when $\Delta t_{i}\rightarrow0$ for each $i\in I^{-}$, the following inequalities are valid:

and

□

Consequently, we obtain the linear relaxation program of ELMP1:

In the constraint of LRP, we substitute $z_{i} = c_{i}^{T}y+d_{i}$ into the inequalities $z_{i}\leq \overline{\psi}(t_{i})$ and $\underline{\psi}(t_{i})\leq z_{i}$, respectively, that is

For each $i = 1, ..., p$, $\zeta_{i} > 0$, LRP is reformulated as

We have derived the linear relaxation problem for ELMP1 through a series of relaxation processes. This relaxation enables us to simplify the problem by loosening certain constraints while providing a reliable lower bound for the global optimal value of ELMP1, facilitating informed decision-making in subsequent optimization steps.

3.   Branch-and-bound algorithm and its property analysis

In this section, we present an efficient deterministic algorithm for solving ELMP1 by combining the linear relaxation problem proposed in Section 2 with subsequent branching operation in Section 3.1 and rectangle compression technique in Section 3.2. The algorithm requires solving a series of linear relaxation problems on the subdivided rectangles of $\mathcal{H}^{0}$. Furthermore, we provide rigorous proofs for the convergence and complexity of the algorithm based on the employed branching operation.

3.1.   Branching operation

For any selected sub-rectangle $\mathcal{H}^{k} = \{t\in \mathbb{R}^{p}|\underline{t}_{i}\leq t_{i}\leq \overline{t}_{i}\}\subseteq \mathcal{H}^{0}$, the following branching rules are given:

(ⅰ) Let $\tau = {\rm argmax}\{\overline{t}_{i}-\underline{t}_{i}, i = 1, ..., p\}$;

(ⅱ) $\mathcal{H}^{k}$ is divided into two sub-rectangles

where $\mathcal{H}_{i} = [t_{i}\in \mathbb{R}|\underline{t}_{i}\leq t_{i} \leq \overline{t}_{i}, i = 1, ..., p, i\neq \tau].$

3.2.   Compression technology

We introduce a compression technique to enhance the convergence speed of the algorithm. When the algorithm iterates to the $k$th time, multiple sub-rectangles are obtained through rectangular subdivision. For any sub-rectangle $\widetilde{\mathcal{H}}^{k}\subseteq \mathcal{H}^{k}$, we further investigate whether ELMP1 over $\widetilde{\mathcal{H}}^{k}$ has a global minimum, where $\widetilde{\mathcal{H}}^{k} = \widetilde{\mathcal{H}}_{1}\times \widetilde{\mathcal{H}}_{2} \times \cdots \times \widetilde{\mathcal{H}}_{p}$ and $\widetilde{\mathcal{H}}_{i} = \{t_{i}\in \mathbb{R}|\underline{t}_{i}\leq t_{i} \leq \overline{t}_{i}, i = 1, ..., p\}.$ The embedded compression technology in the algorithm involves replacing sub-rectangles with smaller rectangles $\widetilde{\mathcal{H}}^{k}$, while ensuring that the global optimal solution of ELMP1 remains intact and unaffected.

Theorem 4. When the algorithm iterates to the $k$th time, let $\widehat{UB}$ be the current best upper bound of the global optimum of ELMP1. Denote

For any sub-rectangle $\widetilde{\mathcal{H}}^{k}\subseteq \mathcal{H}^{k}$, it can be inferred that

(i) If $\Xi > \widehat{UB}$, then there is no global optimal solution over $\widetilde{\mathcal{H}}^{k}$ for ELMP1;

(ii) If $\Xi < \widehat{UB}$, then ELMP1 has no global optimal solution over $\ddot{\mathcal{H}}$, where

with ${\ddot{\mathcal{H}}}_{\iota} = \{t_{\iota}\in \mathbb{R}|\pi_{\iota}\leq t_{\iota}\leq \overline{t}_{\iota}\}\cap \widetilde{\mathcal{H}}_{\iota}$.

Proof. (ⅰ) If $\Xi > \widehat{UB}$, then

Therefore, $\widetilde{\mathcal{H}}^{k}$ does not contain a global optimal solution for ELMP1.

(ⅱ) If $\Xi < \widehat{UB}$, for any $t\in \ddot{\mathcal{H}}$,

Therefore, $\ddot{\mathcal{H}}$ does not contain a global optimal solution for ELMP1. □

3.3.   Branch-and-bound algorithm

The branching operation proposed in Section 3.1 partitions the initial rectangle $\mathcal{H}^{0}$ into smaller sub-rectangles, enabling the algorithm to search for local optimal solutions of ELMP1 over $\mathcal{V}$ that necessarily include the global minimal solution of ELMP1. During the $k$th iteration of the algorithm, we provide some relevant notations. $Q_{k}$ denotes the list of active nodes. For each branching node $\mathcal{H}\in Q_{k}$, $(y(\mathcal{H}), t(\mathcal{H}))$ and $LB(\mathcal{H})$ represent the optimal solution and the optimal value of ${\rm LRP}(\mathcal{H})$, respectively. The current best lower bound for ELMP1 is noted as $LB_{k} = \min\{LB(\mathcal{H}), \mathcal{H}\in Q_{k}\}$. $v_{k}$ represents the objective function value corresponding to the best feasible solution of ELMP1, and the current best upper bound of $v_{k}$ is denoted as $UB$. We choose a divided rectangle $\mathcal{H}^{k}$ from $Q_{k}$ that satisfies $LB(\mathcal{H}) = LB_{k}$, and segment it into two sub-rectangles $\mathcal{H}^{k1}$ and $\mathcal{H}^{k2}$ by branching operation. These sub-rectangles are then added to the set $\mathcal{T}$, and the set $\mathcal{T}$ is updated as $\mathcal{T} = \{\mathcal{T}\backslash\mathcal{H}^{k}\}\bigcup\{\mathcal{H}^{k1}, \mathcal{H}^{k2}\}$. Let $\mathcal{F}$ be the set of feasible points, and $\epsilon$ denotes algorithmic tolerance. In a more precise manner, we can describe the presented outer space branch-and-bound algorithm as follows:

Step 0. Initialize the best known solution as null and the best known upper bound as infinity. Create a root node and initialize the list of active nodes with this root node. Set the algorithmic tolerance to the desired value.

Step 1. Solve a relaxation problem ${\rm LRP(\mathcal{H}^{0})}$ in order to get a lower bound (or prove infeasibility). If problem is feasible, update the incumbent if necessary. Let

If $(y^{0}, t^{0}, z^{0})$ is a feasible solution of ELMP1, then let $UB = \mathcal{L}(y^{0}, t^0, z^{0})$, $\mathcal{F} = \mathcal{F}\bigcup\{y^0, t^0, z^{0}\}$. If $UB-LB_{0}\leq \epsilon$, the algorithm terminates and obtains the global $\epsilon$-optimal solution $(y^{0}, t^{0}, z^{0})$ for ELMP1. Otherwise, denote $\mathcal{T} = \{\mathcal{H}^{0}\}$.

Step 2. Split the current node $\mathcal{H}^{k}$ into two new nodes $\mathcal{H}^{kj}(j = 1, 2)$ and reduce them by using the compression technique, the reduced rectangle is still denoted as $\mathcal{H}^{kj}(j = 1, 2)$ and set $\mathcal{T} = \{\mathcal{T}\backslash \mathcal{H}^{k}\}\bigcup\{\mathcal{H}^{k1}, \mathcal{H}^{k2}\}.$

Step 3. For each child node $\mathcal{H}^{kj}\in \mathcal{T}(j = 1, 2)$, the corresponding optimal value $LB(\mathcal{H}^{kj})$ and optimal solution $(y(\mathcal{H}^{kj}), t(\mathcal{H}^{kj}))$ are obtained by solving ${\rm LRP}(\mathcal{H}^{kj})$. Set $\mathcal{F} = \mathcal{F} \bigcup (y(\mathcal{H}^{kj}), t(\mathcal{H}^{kj}), z(\mathcal{H}^{kj}))$, $(\hat{y}^{k}, \hat{t}^{k}, \hat{z}^{k}) = {\rm arg min}_{(y, t, z)\in \mathcal{F}} \mathcal{L}(y, t, z)$, set $v_{k} = \mathcal{L}(\hat{y}^{k}, \hat{t}^{k}, \hat{z}^{k})$. If $v_{k} < UB$, then update the upper bound $UB = v_{k}$, the current best solution for ELMP1 is updated as $(y^k, t^k, z^{k}) = (\hat{y}^{k}, \hat{t}^{k}, \hat{z}^{k})$, and set $\mathcal{F} = \emptyset$. If $LB(\mathcal{H}^{kj}) > UB$, then remove it from $\mathcal{T}$, i.e. $\mathcal{T} = \mathcal{T}\backslash \{\mathcal{H}^{kj}\}$. Set $Q_{k} = (Q_{k}\backslash \mathcal{H}^{k})\bigcup \mathcal{T}$ and update the lower bound $LB_{k} = \min\{LB(\mathcal{H})|\mathcal{H}\in Q_{k}\}.$

Step 4. Let the list of active nodes $Q_{k+1} = \{\mathcal{H}|UB-LB(\mathcal{H}) > \epsilon, \mathcal{H}\in Q_{k}\}$. If $Q_{k+1}$ is empty, return the best known solution as a global $\epsilon$-optimal solution. Otherwise select an active node $\mathcal{H}^{k+1}\in {\rm arg min}\{LB(\mathcal{H}), \mathcal{H}\in Q_{k+1}\}$. Set $k = k+1$ and go back to Step 2.

Definition 1 provides the concept of the global $\epsilon$-optimal solution involved in the proposed algorithm.

Definition 1. Given $\epsilon > 0$, the feasible solution $(\hat{y}, \hat{t}, \hat{z})$ is considered a global $\epsilon$-optimal solution for ELMP1, if $\mathcal{L}(\hat{y}, \hat{t}, \hat{z})\leq v+\varepsilon$, where $v$ represents the global optimal value of ELMP1.

3.4.   Convergence analysis

This subsection discusses the convergence analysis of the algorithm. Supposing the algorithm is infinite, according to the branching operation, there exists an infinite rectangular sequence $\{\mathcal{H}^{k}\}_{k = 1}^{\infty}$ such that for each $k = 1, 2, ...,$ we have $\mathcal{H}^{k+1}\subseteq\mathcal{H}^{k}\subseteq\cdots\subseteq\mathcal{H}^{0}$, where $\mathcal{H}^{k}\in \mathbb{R}^{p}$. The following Lemma 1 provides a theoretical basis for Theorem 5.

Lemma 1. For any $t\in \mathcal{H}$, when $\overline{t}_{i}-\underline{t}_{i}\rightarrow0$, $i = 1, ..., p$, for which we have $\mathcal{L}(y, t, z)-\underline{\mathcal{L}}(y, t, z)\rightarrow0$.

Proof. It follows from Theorem 2 that

Therefore, for any $t\in \mathcal{H}$, $\mathcal{L}(y, t, z)-\underline{\mathcal{L}}(y, t, z)\rightarrow0$ holds while $\overline{t}_{i}-\underline{t}_{i}\rightarrow0$, $i = 1, ..., p$. □

Theorem 5. Given $\epsilon > 0$, assuming that the feasible domain of ELMP1 is non-empty, the algorithm either obtains a global $\epsilon$-optimal solution of ELMP1 at a finite number of iterations, or produces a sequence of feasible solutions $\{(y^k, t^k, z^{k})\}$, each of whose accumulation points is a global $\epsilon$-optimal solution of ELMP1.

Proof. Assuming that the algorithm terminates within a finite number of iterations, without loss of generality, let us assume that the algorithm terminates at the $k$th iteration, which gets $UB-LB_{k}\leq\epsilon.$ According to Step 3 in the algorithm, we have $UB = v_{k} = \mathcal{L}(\hat{y}^{k}, \hat{t}^{k}, \hat{z}^{k})$, thus

It follows from (3.1) that

Thereby, $(\hat{y}^{k}, \hat{t}^{k}, \hat{z}^{k})$ is a global $\epsilon$-optimal solution to ELMP1.

If the algorithm iterates an infinite number of times and produces an infinite sequence of rectangles $\{\mathcal{H}^{k}\}_{k = 1}^{\infty}$, where $\mathcal{H}^{k} = \prod_{i = 1}^{p}[\underline{t}_{i}^{k}, \overline{t}_{i}^{k}]\in \mathbb{R}^{p}$. Without losing generality, suppose that $\lim\limits_{k\rightarrow \infty}\mathcal{H}^{k} = H^{\infty}$, then for the subsequences $K$ of sequence $\{1, 2, ... \}$ we have

For any $k\in K$, depending on Step 3 of the algorithm, the lower bound is updated to

For any $k\in K$, it follows that $t^{k}\in \mathcal{H}^{k}\subseteq \mathcal{H}$. Therefore, $\{{t^k}\}_{k = 1}^{\infty}$ exists a convergent subsequence $\{{t^k}\}_{k\in K}$, by formula (3.2), the limit of $\{{t^k}\}_{k\in K}$ is in $\mathcal{H}^{\infty}$, let

where $\hat{t}$ is a accumulation point of $\{{t^k}\}_{k\in K}$. Since $\mathcal{L}(y, t, z)$ is continuous, combining with (3.3) we have

For any $k\in K$, $t^k\in \mathcal{H}^{k}$, it follows from Lemma 1 that

Hence, we have

Integrating (3.4)–(3.6), we obtain $\mathcal{L}(\hat{y}, \hat{t}, \hat{z}) = \lim_{k\in K}\mathcal{L}(y^{k}, t^k, z^{k}) = \lim_{k\in K} \underline{\mathcal{L}}(y^{k}, t^k, z^{k}).$ For each $k\in K$ we get $\mathcal{L}(\hat{y}^k, \hat{t}^k, \hat{z}^k) = \sum_{i = 1}^{p}\zeta_{i}{\rm ln}(t_{i}^{k})$, therefore

Because $\{\sum_{i = 1}^{p}\zeta_{i}{\rm ln}(t_{i}^{k})\}_{k\in K}$ is a subsequence of $\{\sum_{i = 1}^{p}\zeta_{i}{\rm ln}(t_{i}^{k})\}_{k = 1}^{\infty}$, following from (3.7), we then obtain $\lim_{k\in K}\sum_{i = 1}^{p}\zeta_{i}{\rm ln}(t_{i}^{k}) = v.$ From the continuity of $\mathcal{L}(y, t, z)$ and formula (3.5), we have

So we get

Combining Eq (3.8), we conclude that $(\hat{y}^{k}, \hat{t}^{k}, \hat{z}^{k})$ is a globally optimal solution to ELMP1. □

3.5.   Complexity analysis

We use $\Omega(\mathcal{H})$ to define the size of the sub-rectangle $\mathcal{H} = \{t\in \mathbb{R}^{p}, \underline{t}_{i}\leq t_{i}\leq \overline{t}_{i}, i = 1, ..., p\}$, i.e.,

In addition,

Lemma 2. Given any $\epsilon > 0$, for a sub-rectangle $\mathcal{H}\subseteq\mathcal{H}^{0}$, if $\Omega(\mathcal{H}) < \epsilon$, then for any optimal solution $(y, t)$ of $LRP(\mathcal{H})$, we have

Proof. Clearly, $(y, t, z)$ is a feasible solution to ELMP1, it follows from Lemma 1 that

□

The following Theorem 6 illustrates an estimation of the maximum number of iterations for the proposed algorithm in the worst case, indirectly indicating that the algorithm will eventually find the global optimal solution for LMP.

Theorem 6. For any $\epsilon_{0}\in (0, 1)$, the algorithm iterates at most

times to obtain a global $\epsilon_{0}$-optimal solution in a worst case.

Proof. Let $\epsilon = \epsilon_{0}$, suppose that during the $k$th iteration, when the algorithm reaches Step 3, $\mathcal{H}\subseteq \mathcal{H}^{0}$ is a rectangle selected by the branching rule to be dissected, which satisfies

then, the algorithm is terminated. In the worst case, when the algorithm terminates at Step 4 on the $k$th iteration, splitting the initial rectangle $\mathcal{H}^{0}$ results in $k+1$ sub-rectangles, with the assumption that $\mathcal{H}^{\iota}$ is any one of these sub-rectangles and satisfies

here, $\iota_{i}$ denotes the initial interval $[\underline{t}_{i}^{0}, \overline{t}_{i}^{0}]$ after $\iota_{i}$ subdivision to produce $[\underline{t}_{i}^{\iota}, \overline{t}_{i}^{\iota}]$. From Lemma 2, it follows that

Since the subdivision $\mathcal{H}^{0}$ yields no more than $\prod_{i = 1}^{p}2^{\iota_{i}}$ sub-rectangles, if every sub-rectangle satisfies (3.9), the algorithm must terminate. By Eq (3.9) we have

Let $\chi_{i} = \lceil\log_{2}\frac{p\theta \lambda\Omega(\mathcal{H}^{0})}{\epsilon_{0}}\rceil, i = 1, ..., p.$ The initial rectangle is split into $k+1$ sub-rectangles and $k+1 = \prod_{i = 1}^{p}2^{\chi_{i}}.$ At this point the algorithm terminates. Thus, the algorithm iterates at most $2^{\sum_{i = 1}^{p}\lceil\log_{2}\frac{p\theta\lambda\Omega(\mathcal{H}^{0})}{\epsilon_{0}}\rceil}-1$ times. Furthermore, it can be derived that

where $(y^*, t^*, z^*)$ is a global optimal solution of ELMP1. For $i\in I^{+}$, $t_{i}^* = c_{i}^{T}y^*+d_{i}$. For $i\in I^{-}$, $z_{i}^* = \frac{1}{t_{i}^*}$ and $z_{i}^* = c_{i}^{T}y^*+d_{i}$. Based on the bounding process, let $\hat{y}^{k}$ be the best feasible solution obtained so far, and denote that $\{f(\hat{y}^{k})\}$ is the decreasing sequence such that $f(\hat{y}^{k})\leq f(y)$. Combined with (3.10) can be obtained

Therefore, the algorithm terminates, and $\hat{y}^{k}$ is a global $\epsilon_{0}$-optimal solution to LMP. □

4.   Numerical experiments

In this section, several problems are employed to demonstrate the feasibility and effectiveness of the proposed algorithm in this paper. All linear programming problems are solved by the dual-simplex method, all tests of the algorithm are carried out using MATLAB9.2.0.538062 (R2017a) on an Inter(R)Core(TM) i5-8250U, CPU@1.60GHz, 4GB memory and 64 bit Windows10 operating system.

Initially, we utilize the existing branch and bound algorithms ^{[10,23,26,28]} and the proposed algorithm to compute the deterministic Problems 1–8 with a predefined convergence tolerance. This is to demonstrate the feasibility of our algorithm. To validate the efficiency of the proposed algorithm, we conduct tests on random Problems 9–11 with a tolerance of $10^{-6}$.

Problem 1 ^[10,26,28]

Problem 2 ^[10,26,28]

Problem 3 ^[10,26,28]

Problem 4 ^[26,28]

Problem 5 ^[10,26,28]

Problem 6 ^[26,28]

Problem 7 ^[23,26,28]

Problem 8 ^[23,26,28]

Problem 9 ^[10,23]

where $\hat{c}_{ij}$ is generated randomly in the interval [0, 1], $i = 1, 2, j = 1, ..., n$. All elements $\hat{a}_{ij}$ of the matrix $\hat{A}$ are randomly generated in the interval [-1, 1], i.e., $\hat{a}_{ij} = 2\hat{\lambda}-1$, where $\hat{\lambda}$ is randomly generated in [0, 1]. The components of $\hat{b}$ is set to $\sum_{j = 1}^{n}\hat{a}_{ij}+2\hat{t}$, where $\hat{t}$ is randomly generated at [0, 1].

Problem 10 ^[23,26]

where $\widetilde{c}_{ij}$ is randomly generated in [0, 1]. $\widetilde{a}_{mj}$ is randomly generated at [-1, 1], $m = 1, ..., M, j = 1, ..., n$. $\widetilde{b}_{m} = \sum_{j = 1}^{n}\widetilde{a}_{mj}+2\tilde{t}$, where $\tilde{t}$ is randomly generated at [0, 1].

Problem 11 ^[26,28]

where $\overline{c}_{ij}$, and $\overline{d}_{i}$ are randomly generated in [0, 1], $i = 1, ..., p, j = 1, ..., n$. Each element of the matrix $\overline{A}$ and $\overline{\alpha}^{i}$ are randomly generated in [-1, 1]. The components of the vector $\overline{b}$ are generated by $\sum_{j = 1}^{n}\bar{a}_{ij}+2\bar{t}$, where $\bar{t}$ is generated randomly at [0, 1].

Table 1 shows the numerical comparison between some algorithms and our algorithm on Problems 1–8.

For stochastic Problems 9–11, we solve 10 randomly generated problems for each set of parameters $(p, m, n)$ and place their average number of iterations and average CPU time in Tables 2–6. Specifically, Problem 9 represents an LMP with only two linear functions and exponents of 1. Table 2 shows the results of numerical comparisons between our algorithm and the algorithms proposed in ^[10,23]. Problem 10 is an LMP with multiple linear functions and exponents of 1. Tables 3 and 4 display the numerical results of our algorithm compared with the algorithms in ^[23,26]. Additionally, Figures 1 and 2 plot some of the data results in Table 3. Problem 11 is an LMP with real exponents and multiple linear functions. Tables 5 and 6 show the numerical results of our algorithm compared with the algorithms in ^[26,28]. Figures 3–6 depict some of the data results from Tables 5 and 6.

For convenience, the symbols in the table headers in Tables 1–6 are specified as follows: Opt.val: the global optimum of the tested problem; Iter: the number of iterations of the algorithm; Time: the CPU time in seconds; Avg.Iter: the average number of iterations of the 10 randomly generated problems; Std.Iter: the standard deviation of the number of iterations; Avg.Time: the average CPU time of the 10 randomly generated problems; Std.Time: the standard deviation of the average CPU time; $p$: the number of linear functions; $m$: the number of constraints; $n$: the dimensionality of decision variables.

As can be seen from the numerical results in Table 1, our algorithm effectively calculates the global optimal solutions and optimal values for low-dimensional Problems 1–8. In comparison to the algorithms proposed in ^[10,26,28], our algorithm demonstrates shorter computation time and fewer iterations. Most deterministic problems require only one iteration, with a maximum of nine iterations, indicating the feasibility of our algorithm.

Upon observing Table 2, it is evident that our algorithm exhibits a lower average number of iterations and shorter average CPU time compared to the algorithm proposed in ^[10] for medium-scale Problem 9. The primary reason for this disparity is that our algorithm solves the problem in the $p$-dimensional space, whereas the algorithm in ^[10] tackles it in the $n$-dimensional space. In comparison to the algorithm presented in ^[23], it is apparent that the iterations of the algorithm in ^[23] is generally lower than our algorithm. However, when considering the average CPU time, our algorithm outperforms in terms of efficiency.

From Table 3, it is easy to see that the proposed algorithm is more efficient in solving Problem 10. First, when compared to the algorithm in ^[26], our algorithm consumes less time and requires fewer iterations. Second, in comparison to the algorithm in ^[23], our algorithm exhibits a lower number of iterations for parameter sets where $p$ is fixed at 7. Although the algorithm in ^[23] may have less iterations than ours for some parameter groups, our average CPU time is still lower than that of ^[23]. For instance, when $p = 4$, our algorithm demonstrates higher iterations than ^[23] as $m$ and $n$ vary, yet our CPU time is shorter than theirs.

To visualize the effectiveness of the proposed algorithm in this paper, Figures 1 and 2 depict the trends of the average number of iterations and average CPU time for fixing $(m, n)$ to (10, 20), (100,200) and $p$ to change from 4 to 7 in Problem 10, respectively. From Figure 1, when $(m, n) = (10, 20)$ (indicated by the solid line), the green solid line is always above the blue solid line, indicating that the number of iterations of the algorithm in ^[26] is higher than ours. On the other hand, the red solid line is above the blue solid line only after $p = 6$, which means that the iterations of the algorithm in ^[23] are higher than ours after $p > 6$. When $(m, n) = (100,200)$ (indicated by the dashed line), both the red dashed line and the green dashed line are above the blue dashed line, indicating that the number of iterations of our algorithm is lower than the other two algorithms. In the vertical direction, the blue dashed line is always higher than the blue solid line, but the vertical distance between them is shorter than the corresponding vertical distance of the other two algorithms. This implies that as $(m, n)$ increases, the number of iterations of our algorithm exhibits a slight increase, but the magnitude of the increase is not significant. Based on Figure 2, we observe that when $(m, n) = (10, 20)$, $p = 4, 5, 6$, the three solid lines approximately coincide and when $p = 7$, it is obvious that our algorithm takes less time. When $(m, n) = (100,200)$, the distance between the red dashed line and the blue dashed line becomes larger as $p$ increases. The time taken by our algorithm increases as $(m, n)$ increases, but not significantly.

Table 4 gives the numerical results for solving the large-scale Problem 10. It can be observed that our algorithm is more time-saving compared to the algorithm in ^[26]. In contrast to the algorithm in ^[23], the first two sets of parameters yield better results than ours. However, for the case of $(m, n) = (3, 10, 1000)$, the number of iterations in our algorithm is comparable to that of ^[23], but our algorithm exhibits lower time consumption. When $(p, m, m) = (3, 10, 2000)$, although our algorithm requires fewer iterations than that of ^[23], the CPU time is slightly longer. When $(p, m, n) = (4, 10, 1000), (4, 10, 2000)$, respectively, our algorithm demonstrates clear advantages. Specifically, the algorithm in ^[23] requires twice as many iterations and three times as much CPU time as ours.

Table 5 shows the results of a medium-scale test for Problem 11. Compared to the algorithms in ^[26,28], our algorithm exhibits a distinct advantage. Figures 3 and 4 show the average number of iterations and the average CPU time among the algorithms in ^[26,28] and our algorithm when $(m, n)$ is fixed (10, 20), (100,200) on Problem 11, respectively. From Figure 3, it can be observed that the average number of iterations in the algorithm proposed by ^[26] is higher than our algorithm when $(m, n) = (10, 20)$ and lower than our algorithm only when $p = 5$. When $(m, n) = (100,200)$, the average number of iterations in our algorithm is lower than those of the algorithms in ^[26,28]. Moreover, our algorithm significantly outperforms the algorithms in ^[26,28] when $(m, n) = (100,200)$. From Figure 4, we find that while $(m, n) = (10, 20), (100,200)$, the average CPU time of the proposed algorithm is less than the other two algorithms, and the trend of time change is not obvious as $p$ increases.

Table 6 shows the numerical results of a large-scale Problem 11, we can see that the number of linear functions $p$ is much smaller than $n$. Fixing $(p, m) = (2, 10)$, the number of iterations of our algorithm decreases as the decision variable $n$ increases, while the number of iterations of the algorithms in the references ^[26,28] decreases and then increases as the decision variable $n$ increases, as reflected in Figure 5. Moreover, we find that the algorithm in ^[26,28] consume more CPU time than our algorithm. Furthermore, we conclude that CPU time increases with an increase in the decision variable or linear function. Figure 6 provides a clearer picture of this conclusion.

5.   Conclusions

We propose a new branch-and-bound algorithm for solving LMP in outer space. First, the original problem is reduced to an equivalent problem by introducing auxiliary variables. Then, the equivalent problem is further simplified by leveraging the properties of exponential and logarithmic functions.The focus switches to resolving the second equivalent problem. Afterwards, the nonlinear constraints are linearly relaxed, and the objective function of the equivalent problem is linearly approximated to obtain the linear relaxation programming problem. Consequently, the proposed rectangular compression technique is embedded into the branch-and-bound framework, leading to the development of the outer space branch-and-bound algorithm. Furthermore, we demonstrate the computational complexity to estimate the maximum number of iterations of the algorithm in the worst case. Finally, in order to verify the feasibility and efficiency of the algorithm, some deterministic and random problems are tested. The experimental results show that the algorithm exhibits good computational performance, particularly for large-scale random LMP where the number of linear functions $p$ is significantly smaller than the decision variable $n$.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (11961001), the Construction Project of first-class subjects in Ningxia higher Education(NXYLXK2017B09), the Basic discipline research projects supported by Nanjing Securities(NJZQJCXK202201).

Conflict of interest

The authors declare no conflicts of interest.