A generalized trapezoidal fuzzy number A ̃ is characterized by a quadruple (a, b, c, d) with 0 < a ≤ b ≤ c ≤ d and a height μ ∈ (0,1]. It is denoted as A ̃ = ⟨(a, b, c, d); μ⟩. Its membership function is x − a μ. a ≤ x < b, b − a μ b ≤ x ≤ c, μ (x) = (14) d − x μ. c < x ≤ d, d − c { 0 otherwise.

In the mTSP under study, the cost of traversing any arc (i, j) is a positive generalized trapezoidal fuzzy number

c ̃ij = 〈(aij, bij, cij, dij); μij〉.

Definition 1 (parametric defuzzification function ([27])). For a generalized trapezoidal fuzzy number A ̃ = ⟨(a, b, c, d); μ⟩ and a parameter α ∈ [0, μ], the function M (A ̃) is defined as:

μ − α M (A ̃) = (a + d)(μ − α) + (b + c − a − d) . (15)

The parameter α acts as a fuzzy cut level. By increasing α, the contribution of the core region of the fuzzy number (where membership is high) is emphasised, while the influence of the boundaries (where membership is low) diminishes.

Theorem 1 (continuity and differentiability). As a function of α, M (A ̃) is continuous on [0, μ] and infinitely differentiable on (0, μ).

Proof: Expanding Eq. (15) yields

b + c − a − d ̃ 2 2 Mα(A) = (a + d)μ − (a + d)α + (μ − α )

b + c − a − d b + c − a − d = (a + d)μ + μ − (a + d)α − α .

2 2μ

This expression is a second-degree polynomial in α. Since polynomials are continuous and possess derivatives of all orders on R, the claim follows. The first and second derivatives are

dMα b + c − a − d = −(a + d) − α, dα μ d Mα b + c − a − d = − .

dα μ

Theorem 2 (strict monotonicity). For any positive generalized trapezoidal fuzzy number A ̃ (a > 0), the function M (A ̃) is strictly decreasing with respect to α on [0, μ].

dMα

Proof: We show that < 0 for all α ∈ (0, μ).

dα

dMα b + c − a − d = −(a + d) − α.

dα μ This derivative is an affine function of α. Its value at the endpoints of the interval is

dMα |α=0 = −(a + d) < 0, dα dMα |α=μ = −(a + d) − (b + c − a − d) = −(b + c) < 0.

dα (Since a + d > 0 and b + c > 0), both endpoint values are strictly negative. Because the derivative is linear in α, its sign on the whole interval [0, μ] is determined by its values at the endpoints and its slope.

b+c−a−d

I. If b + c − a − d ≥ 0, the slope − ≤ 0; thus, the derivative is a non-increasing function of α. Its maximum on [0, μ] occurs at α = 0, which is −(a + d) < 0. Hence the derivative remains negative.

II. If b + c − a − d < 0, the slope is positive, so the derivative is increasing. Its maximum on [0, μ] occurs at α = μ, which equals −(b + c) < 0. Therefore, the derivative is negative throughout.

dMα

In both cases, < 0 for all α ∈ (0, μ). Consequently, Mα is strictly decreasing on [0, μ].

dα

Theorem 3 (positivity and boundedness). Let A ̃ be a positive generalized trapezoidal fuzzy number. Then

0 = M (A ̃) < M (A ̃) ≤ M (A ̃) for all α ∈ [0, μ).

μ α 0

Moreover, M (A ̃) > 0 holds for every α < μ without any additional condition .

Proof: The strictly decreasing property from Theorem 2 immediately gives M (A ̃) < M (A ̃) < M (A ̃) for α ∈

μ α 0

(0, μ). Substituting α = μ into Factor (15) yields M (A ̃) = 0. To prove that M (A ̃) > 0 for α < μ, Factor (15) as

μ + α M (A ̃) = (μ − α) [a + d + (b + c − a − d) ].

The factor μ − α is strictly positive for α < μ. Let B(α) denote the bracketed expression. B(α) is a linear

b+c−a−d

function of α with slope .

I. If b + c − a − d ≥ 0, B is non-decreasing. Its minimum on [0, μ] occurs at α = 0:

b + c − a − d a + d b + c B(0) = a + d + = + > 0.

2 2 2 II. If b + c − a − d < 0, B is decreasing. Its minimum occurs at α = μ:

a + d + (b + c − a − d). 1 = b + c > 0.

In both cases, B(α) > 0 for every α ∈ [0, μ]. Hence M (A ̃) = (μ − α)B(α) > 0 for all α < μ. The upper bound

b+c−a−d

is M (A ̃) = (a + d)μ + μ.

It should be noted that M (A ̃) = 0 is a limiting value that is not attainable in practice, because the parameter α must be chosen strictly less than min μij to prevent any crisp cost from becoming zero or negative.

ij

Theorem 4. The function M (A ̃) is I. Strictly concave on [0, μ] if b + c > a + d, II. Linear if b + c = a + d, III. Strictly convex if b + c < a + d.

d Mα b+c−a−d

Proof: The second derivative is = − . The sign of this constant determines the curvature:

dα2 μ

I. If b + c − a − d > 0, the second derivative is negative, implying strict concavity .

II. If b + c − a − d = 0, the second derivative is zero, so the function is linear.

III. If b + c − a − d < 0, the second derivative is positive, giving strict convexity.

To illustrate the theoretical properties, we evaluate Mα for three arcs taken from the numerical example (Section 4). The results are summarized in Table 1.

For c ̃14 = ⟨(5.7,7.3,9.6,11.2); 0.8⟩, we have a + d = 16.9 = b + c, hence the function is linear. For c ̃15 = ⟨(4.2,6.1,8.4,10.7); 0.5⟩, a + d = 14.9 > 14.5 = b + c, so the function is convex. For c ̃12 = ⟨(3.5,5.2,7.8,9.1); 0.7⟩, a + d = 12.6 < 13.0 = b + c, making the function concave. In all cases, monotonic decrease with α is evident. Positivity is maintained for α < μ; the zero value at α = 0.5 for c ̃15 is consistent with α = μ15, confirming that the cut level must be kept below the minimum height across all arcs .

The established properties endow the defuzzification function Mα with a clear decision-theoretic interpretation. The parameter α controls the degree of optimism/pessimism embedded in the crisp equivalent cost:

I. At α = 0, the entire support of the fuzzy number contributes to the crisp value, yielding the most conservative (i.e., highest) cost estimate.

II. As α increases, the influence of low-membership regions diminishes, and the crisp cost decreases monotonically.

III. As α → μ, the crisp cost approaches zero, which represents an extreme, unattainable optimistic scenario. Consequently, α serves as a continuous risk-adjustment lever: A risk-averse decision-maker would choose a small α, while a risk-seeking decision-maker might select a larger α. The convexity/concavity of the function also indicates that the rate of cost reduction is not constant but depends on the asymmetry of the trapezoidal number.

These theoretical guarantees ensure that the defuzzification stage does not introduce artifacts (e.g., negative costs) that could compromise the MILP solution. In the following section, the complete two-stage algorithm is applied to a 7-city mTSP instance, and the impact of α on both the cost and the structure of optimal tours is thoroughly analyzed .

The instance consists of n=7 cities, with city 1 serving as the common depot. There are m=2 salesmen. The minimum and maximum number of non‑depot cities that may be assigned to a single salesman are set to k=2 and w=6, respectively. With only six non‑depot cities (2, … ,7), these bounds ensure feasibility while still permitting imbalanced allocations (e.g., a tour of length 4 together with a tour of length 2, or the symmetric 3 + 3 split).

The travel costs are modelled as generalized trapezoidal fuzzy numbers c ̃ij = 〈(aij, bij, cij, dij); μij〉 with aij > 0. The complete set of directed‑arc costs is given below in a row‑by‑row fashion. For completeness, the diagonal elements (self‑loops) are set to a dummy value ⟨(100,100,100,100); 1.0⟩ and are never selected because of the constraint xii = 0.

Row 1 (from city 1 – depot):

c ̃12 = 〈(3.5, 5.2, 7.8, 9.1); 0.7〉.

c ̃13 = 〈(2.8, 4.1, 6.3, 8.5); 0.6 〉.

c ̃14 = 〈(5.7, 7.3, 9.6, 11.2); 0.8 〉.

c ̃15 = 〈(4.2, 6.1, 8.4, 10.7); 0.5〉.

c ̃16 = 〈(7.9, 9.5, 12.3, 14.8); 0.9〉.

c ̃17 = 〈(6.4, 8.2, 10.9, 13.5); 0.7〉.

Row 2 (from city 2):

c ̃21 = 〈(2.7, 4.3, 6.8, 8.9); 0.8〉.

c ̃23 = 〈(5.9, 7.6, 10.1, 12.3); 0.6〉.

c ̃24 = 〈(4.8, 6.7, 9.2, 11.8); 0.5〉.

c ̃25 = 〈(8.6, 10.9, 13.7, 16.2); 0.8〉.

c ̃26 = 〈(7.2, 9.4, 12.1, 14.9); 0.7〉.

c ̃27 = 〈(10.8, 13.2, 16.3, 19.1); 0.9〉.

Row 3 (from city 3):

c ̃31 = 〈(4.1, 5.9, 8.3, 10.6); 0.6〉.

c ̃32 = 〈(6.8, 8.7, 11.4, 13.9); 0.7〉.

c ̃34 = 〈(8.9, 11.2, 14.1, 16.8); 0.8〉.

c ̃35 = 〈(7.5, 9.8, 12.6, 15.3); 0.5〉.

c ̃36 = 〈(11.7, 14.3, 17.5, 20.8); 0.9〉.

c ̃37 = 〈(10.2, 12.9, 16.1, 19.3); 0.7〉.

Row 4 (from city 4):

c ̃41 = 〈(5.6, 7.4, 10.2, 12.7); 0.7〉.

c ̃42 = 〈(8.3, 10.6, 13.7, 16.4); 0.8〉.

c ̃43 = 〈(10.1, 12.8, 16.2, 19.5); 0.9〉.

c ̃45 = 〈(12.9, 15.7, 19.2, 22.8); 0.6〉.

c ̃46 = 〈(11.3, 14.1, 17.6, 21.2); 0.5〉.

c ̃47 = 〈(15.8, 18.9, 22.7, 26.4); 0.8〉.

Row 5 (from city 5):

c ̃51 = 〈(7.2, 9.5, 12.8, 15.6); 0.8〉.

c ̃52 = 〈(10.4, 13.2, 16.7, 19.9); 0.9〉.

c ̃53 = 〈(12.7, 15.6, 19.3, 22.8); 0.5〉.

c ̃54 = 〈(15.3, 18.4, 22.3, 26.1); 0.7〉.

c ̃56 = 〈(17.9, 21.2, 25.1, 29.3); 0.8〉.

c ̃57 = 〈(16.4, 19.7, 23.6, 27.5); 0.6〉.

Row 6 (from city 6):

c ̃61 = 〈 (9.1, 11.8, 15.3, 18.6); 0.9〉.

c ̃62 = 〈(12.6, 15.7, 19.6, 23.2); 0.5〉.

c ̃63 = 〈(15.2, 18.6, 22.7, 26.5); 0.7〉.

c ̃64 = 〈(18.1, 21.5, 25.8, 29.7); 0.8〉.

c ̃65 = 〈(20.7, 24.3, 28.6, 32.9); 0.6〉.

c ̃67 = 〈(22.9, 26.7, 30.9, 35.2); 0.9〉.

Row 7 (from city 7):

c ̃71 = 〈 (11.3, 14.2, 18.1, 21.7); 0.5〉.

c ̃72 = 〈(14.9, 18.3, 22.8, 26.9); 0.7〉.

c ̃73 = 〈(17.8, 21.4, 25.9, 30.2); 0.8〉.

c ̃74 = 〈(20.9, 24.8, 29.5, 34.1); 0.9〉.

c ̃75 = 〈(23.8, 27.9, 32.7, 37.4); 0.6〉.

c ̃76 = 〈(26.7, 30.9, 35.8, 40.6); 0.5〉.

The parametric defuzzification function Mα, defined in Eq. (15) of Section 3, was applied to each fuzzy cost c ̃ij. The resulting crisp cost matrices are displayed in Tables 2–7 for α = 0.0, 0.1, 0.2, 0.3, 0.4, 0.5. All off‑diagonal entries are rounded to two decimal places. A clear monotonic decrease with increasing α is observed; at α = 0.5, several costs become zero, reflecting the most optimistic attitude (see Theorem 3). The diagonal entries, though irrelevant, are left at the dummy value of the crisp equivalent of the self‑loop fuzzy cost, which is 200 for α =0.0 and gradually decreases to 100 at α =0.5.

As α increases, all off‑diagonal costs strictly decrease (or remain zero). The zero entries at α = 0.5 correspond to arcs whose fuzzy numbers have μ = 0.5 and for which the crisp value vanishes exactly at the boundary α = μ. In practice, α should be chosen strictly smaller than min μij to avoid zero costs; here α = 0.5 is the theoretical limit and is included only for illustration.

Each crisp matrix was supplied to the corrected MILP model of Section 2. The solver returned an optimal solution (exit flag = 1) for every α. The solutions were verified to ensure that all non‑depot cities are visited exactly once, tour lengths lie within [3], [13], and no subtour exists. The results are summarized in Table 8.

The total cost decreases strictly as α rises, empirically confirming Theorem 2. The solver times are all under one second, reflecting the small size of the instance and the efficiency of the MILP formulation.

A notable observation is the structural change of the optimal tour set at α = 0.3. For conservative attitudes (α = 0.0,0.1,0.2), the solution consists of one longer tour (1–4–6–2–3–1, 4 cities) and one shorter tour (1–5– 7–1, 2 cities). From α = 0.3 onward, the pattern reverses: the shorter tour becomes 1–5–3–1 (2 cities), and the longer tour covers the remaining cities (1–7–6–2–4–1, 4 cities). This switch is driven by the differential reduction of arc costs as α increases. For instance, the fuzzy cost c ̃15 has height μ = 0.5 and a relatively broad support; its defuzzified value drops from 7.35 at α = 0.0 to 2.92 at α = 0.3 and eventually to zero at α = 0.5. Simultaneously, c ̃53 declines to zero as well. Consequently, the cycle 1–5–3–1 becomes increasingly attractive, freeing salesman 2 to serve the remaining four cities.

The auxiliary variables ui from the MILP solution confirm the strict sequencing. For the solution at α = 0.0, the values are:

Tour 1: u4 = 1, u6 = 2, u2 = 3, u3 = 4, Tour 2: u5 = 1, u7 = 2.

Both sequences satisfy uj = ui + 1 for each traversed arc, ensuring that the number of non‑depot cities in a tour equals the u-value of its last city. The satisfaction of these sequencing constraints demonstrates that the strict sequencing constraints correctly enforce the minimum tour length k=2 (since u3 = 4 ≥ 2 and u7 = 2 ≥ 2) and maximum length w=6.

The parameter α acts as an intuitive, continuous dial that translates the decision‑maker’s attitude toward risk into concrete routing plans. A small α (e.g., 0.0–0.2) corresponds to a conservative stance: The full spread of the fuzzy numbers is considered, resulting in higher total cost estimates but routes that are likely to remain feasible under worst‑case conditions. Larger α values (0.3–0.5) represent an increasingly optimistic outlook, producing lower cost estimates and more aggressive routing choices that exploit the core of the fuzzy numbers.

The abrupt change in the optimal tour configuration at α= 0.3 carries an important practical message: Even a moderate shift in the perceived level of uncertainty can fundamentally alter the optimal operational plan. A logistics manager who slightly relaxes the conservatism of cost estimates may find that a completely different assignment of cities to vehicles becomes preferable. The proposed methodology enables such managers to quantify and visualize the trade‑off between cost and robustness, supporting informed decision‑making under uncertainty.

Moreover, the analysis highlights the necessity of selecting α strictly below the minimum fuzzy height of all arcs. In the present instance, min μij = 0.5; choosing α = 0.5 yields several zero costs, which may be unrealistic. In practice, one would select α values in the interval [0, μmin), for example 0.4 or lower, to avoid degenerate costs while still exploring the full spectrum of risk attitudes.

The numerical experiments confirm the correctness and practical utility of the proposed two‑stage method. The corrected MILP formulation reliably produces feasible, subtour‑free solutions that respect all cardinality constraints. The defuzzification function Mα behaves precisely as predicted by the theory: It is continuous, strictly decreasing in α, and produces positive crisp costs for α < μ. The sensitivity analysis reveals that α is not merely a numerical scaling factor but can induce qualitative changes in optimal routing, thereby serving as an effective risk‑management lever. The next section concludes the paper and outlines promising directions for future research.