Mohamed Hebiri

Johannes C Lederer

How Correlations Influence Lasso Prediction

Keywords: Correlations, Lars Algorithm, Lasso, Restricted Eigenvalue, Tuning Parameter

We study how correlations in the design matrix influence Lasso prediction. First, we argue that the higher the correlations are, the smaller the optimal tuning parameter is. This implies in particular that the standard tuning parameters, that do not depend on the design matrix, are not favorable. Furthermore, we argue that Lasso prediction works well for any degree of correlations if suitable tuning parameters are chosen. We study these two subjects theoretically as well as with simulations.

Introduction

Although the Lasso estimator is very popular and correlations are present in many of its diverse applications, the influence of these correlations is still not entirely understood. Correlations are surely problematic for parameter estimation and variable selection. The influence of correlations on prediction, however, is far less clear.

Let us first set the framework for our study. We consider the linear regression model

Y = Xβ 0 + σǫ, (1)
where Y ∈ R n is the response vector, X ∈ R n×p is the design matrix, ǫ ∈ R n is the noise and σ ∈ R + is the noise level. We assume in the following that the noise level σ is known and that the noise ǫ obeys an n dimensional normal distribution with covariance matrix equal to the identity. Moreover, we assume that the design matrix X is normalized, that is, X T X jj = n for 1 ≤ j ≤ p. Three main tasks are then usually considered: estimating β 0 (parameter estimation), selecting the nonzero components of β 0 (variable selection), and estimating Xβ 0 (prediction). Many applications of the above regression model are high dimensional, that is, the number of variables p is larger than the number of observations n but are also sparse, that is, the true solution β 0 has only few nonzero entries. A computationally feasible method for the mentioned tasks is, for instance, the widely used Lasso estimator introduced in [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]:

β := arg min β∈R p Y -Xβ 2 2 + λ β 1 .
In this paper, we focus on the prediction error of this estimator for different degrees of correlations. The literature on the Lasso estimator has become very large, we refer the reader to the well written books [BC11, [START_REF] Bühlmann | Statistics for High Dimensional Data[END_REF][START_REF] Hastie | The elements of statistical learning[END_REF] and the references therein.

Two types of bounds for the prediction error are known in the theory for the Lasso estimator. On the one hand, there are the so called fast rate bounds (see [BRT09, [START_REF] Bunea | Sparsity oracle inequalities for the Lasso[END_REF][START_REF] Van De Geer | On the conditions used to prove oracle results for the lasso[END_REF] and references therein). These bounds are nearly optimal but imply restricted eigenvalues or similar conditions and therefore only apply for weakly correlated designs. On the other hand, there are the so called slow rate bounds (see [START_REF] Huang | Risk of penalized least squares, greedy selection and L1 penalization for flexible function libraries[END_REF][START_REF] Koltchinskii | Nuclear norm penalization and optimal rates for noisy low rank matrix completion[END_REF][START_REF] Massart | The Lasso as an ℓ 1 -ball model selection procedure[END_REF][START_REF] Rigollet | Exponential Screening and optimal rates of sparse estimation[END_REF]). These bounds are valid for any degree of correlations but -as their name suggests -are usually thought of as unfavorable.

Regarding the mentioned bounds, one could claim that correlations lead in general to large prediction errors. However, recent results in [START_REF] Van De Geer | The lasso, correlated design, and improved oracle inequalities[END_REF] suggest that this is not true. It is argued in [START_REF] Van De Geer | The lasso, correlated design, and improved oracle inequalities[END_REF] that for (very) highly correlated designs, small tuning parameters can be chosen and favorable slow rate bounds are obtained. In the present paper, we provide more insight into the relation between Lasso prediction and correlations. We find that the larger the correlations are, the smaller the optimal tuning parameter is. Moreover, we find both in theory and simulations that Lasso performs well for any degree of correlations if the tuning parameter is chosen suitably.

We finally give a short outline of this paper: we first discuss the known bounds on the Lasso prediction error. Then, after some illustrating numerical results, we study the subject theoretically. We then present several simulations, interpret our results and finally close with a discussion.

Known Bounds for Lasso Prediction

To set the context of our contribution, we first discuss briefly the known bounds for the prediction error of the Lasso estimator. We refer to the books [BC11] and [START_REF] Bühlmann | Statistics for High Dimensional Data[END_REF] for a detailed introduction to the theory of the Lasso.

Fast rate bounds, on the one hand, are bounds proportional to the square of the tuning parameter λ. These bounds are only valid for weakly correlated design matrices. We first recall the corresponding assumption. Let a be a vector in R p , J a subset of {1, . . . , p}, and finally a J the vector in R p that has the same coordinates as a on J and zero coordinates on the complement J c . Denote the cardinality of a given set by | • |. For a given integer s, the Restricted Eigenvalues (RE) assumption introduced in [BRT09] reads then Assumption RE(s):

φ(s) := min J 0 ⊂{1,...,p}:|J 0 |≤s min ∆ =0: ∆ J c 0 1 ≤3 ∆ J 0 1 X∆ 2 √ n ∆ J 0 2 > 0.
The integer s plays the role of a sparsity index and is usually comparable to the number of nonzero entries of β 0 . More precisely, to obtain the following fast rates, it is assumed that s ≥ s, where s := |{j : (β 0 ) j = 0}|. Also, we notice that φ(s) ≈ 0 corresponds to correlations. Under the above assumption it holds (see for example Bickel et al.

[BRT09] and more recently Koltchinskii et al. [START_REF] Koltchinskii | Nuclear norm penalization and optimal rates for noisy low rank matrix completion[END_REF]):

X( β -β 0 ) 2 2 ≤ λ 2 s nφ 2 (s) (2)
on the set T := sup β 2σ|ǫ T Xβ| β 1 ≤ λ . Similar bounds, under slightly different assumptions, can be found in [START_REF] Van De Geer | On the conditions used to prove oracle results for the lasso[END_REF]. Usually, the tuning parameter λ is chosen proportional to σ n log(p). For fixed φ, the above rate then is optimal up to a logarithmic term (see [START_REF] Bunea | Aggregation for Gaussian regression[END_REF]Theorem 5.1]) and the set T has a high probability (see Section 3.2.2). For correlated designs, however, this choice of the tuning parameter is not suitable. This is detailed in the following section.

Slow rate bounds, on the other hand, are bounds only proportional to the tuning parameter λ. These bounds are valid for arbitrary designs, in particular, they are valid for highly correlated designs. The result [KTL11, Eq. (2.3) in Theorem 1] yields in our setting

X( β -β 0 ) 2 2 ≤ 2λ β 0 1 (3)
on the set T . Similar bounds can be found in [START_REF] Huang | Risk of penalized least squares, greedy selection and L1 penalization for flexible function libraries[END_REF] (for a related work on a truncated version of the Lasso), in [MM11, Theorem 3.1] (for estimation of a general function in a Banach space), and in [RT11, Theorem 4.1] (which also applies to the non-parametric setting; note that the corresponding bound can be written in the form above with arbitrary tuning parameter λ). We note that these bounds depend on β 0 1 instead of s. Moreover, they depend on λ to the first power, and these bounds are therefore considered unfavorable compared to the fast rate bounds.

The mentioned bounds are only useful for sufficiently large tuning parameters such that the set T has a high probability. This is crucial for the following. We show that the higher the correlations, the larger the probability of T is. Correlations thus allow for small tuning parameters; this implies for correlated designs, via the factor λ in the slow rate bounds, favorable bounds even though no fast rate bounds are available.

Remark 2.1. The slow rate bound (3) can be improved if β 0 1 is large. Indeed, we proof in the Appendix that

X( β -β 0 ) 2 2 ≤ 2λ min β 0 1 , ( β -β 0 ) J 0 1
on the set T . That is, the prediction error can be bounded both with β 0 1 and with the ℓ 1 estimation error restricted on the sparsity pattern of β 0 . The latter term can be considerably smaller than β 0 1 (in particular for weakly correlated designs). A detailed analysis of this observation, however, is not within the scope of the present paper.

The Lasso and Correlations

We show in this section that correlations strongly influence the optimal tuning parameters. Moreover, we show that -for suitably chosen tuning parameters -Lasso performs well in prediction for different levels of correlations. For this, we first present simulations where we compare Lasso prediction for an initial design with Lasso prediction for an expanded design with additional impertinent variables. Then, we discuss the theoretical aspects of correlations. We introduce, in particular, a simple and illustrating notion about correlations. Further simulations finally confirm our analysis.

The Lasso on Expanded Design Matrices

Is Lasso prediction becoming worse when many impertinent variables are added to the design? Regarding the bounds and the usual value of the tuning parameter λ described in the last section, one may expect that many additional variables lead to notably larger optimal tuning parameters and prediction errors. However, as we see in the following, this is not true in general.

Let us first describe the experiments.

Algorithm 1 We simulate from the linear regression model (1) and take as input the number of observations n, the number of variables p, the noise level σ, the number of nonzero entries of the true solution s := {j : (β 0 ) j = 0} and finally a correlation factor ρ ∈ [0, 1). We then sample the n independent rows of the design matrix X from a normal distribution with mean zero and covariance matrix with diagonal entries equal to 1 and off-diagonal entries equal to ρ, and we normalize X such that (X ⊤ X) jj = n for 1 ≤ j ≤ n. Then, we define (β 0 ) i := 1 for 1 ≤ i ≤ s and (β 0 ) i := 0 otherwise, sample the error ǫ from a standard normal distribution and compute the response vector Y according to (1). After calculating the Lasso solution β, we finally compute the prediction error X( β -β 0 ) 2 2 for different tuning parameters λ and find the optimal tuning parameter, that is, the tuning parameter that leads to the smallest prediction error.

Algorithm 2 This algorithm only differs from the above algorithm in one point. In an additional step after the initial design matrix X is sampled, we add for each column X (j) of the initial design matrix p -1 columns sampled according to X (j) + ηN. We finally normalize the resulting matrix. The parameter η controls the correlation among the added columns and the initial columns and N is a standard normally distributed random vector. Compared to the initial design, we have now a design with p 2p additional impertinent variables.

Several algorithms for computing a Lasso solution have been proposed: For example, using interior point methods [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], using homotopy parameters [EHJT04, OPT00, Tur05], or using a so-called shooting algorithm [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF][START_REF] Fu | Penalized regressions: the bridge versus the lasso[END_REF][START_REF] Friedman | Pathwise coordinate optimization[END_REF]. We use the LARS algorithm introduced in [EHJT04], since, among others, Bach et al. [BJMO11, Section 1.7.1] have confirmed the good behavior of this algorithm when the variables are correlated. 2 for 1000 iterations as a function of the tuning parameter λ. The blue, dashed line corresponds to Algorithm 1, where A stands for the initial design matrices. The blue, dotted lines give the confidence bounds. The red, solid line corresponds to Algorithm 2, where A represents the extended matrices. The faint, red lines give the confidence bounds. The parameters for the algorithms are given in the header. The mean of the optimal tuning parameters is 3.62 ± 0.01 for Algorithm 1 and 3.63 ± 0.01 for Algorithm 2. These values are represented by the blue and red vertical lines.

Results

We did 1000 iterations of the above algorithms for different λ and with n = 20, p = 40, s = 4, σ = 1, ρ = 0 and with η = 0.001 (Figure 1) and η = 0.1 (Figure 2). We plot the means P E of the prediction errors as a function of λ. The blue, dashed curves correspond to the initial designs (Algorithm 1), the red, solid 2 for 1000 iterations as a function of the tuning parameter λ. The blue, dashed line corresponds to Algorithm 1, where A stands for the initial design matrices. The blue, dotted lines give the confidence bounds. The red, solid line corresponds to Algorithm 2, where A stands for the extended matrices. The faint, red lines give the confidence bounds. The parameters for the algorithms are given in the header. The mean of the optimal tuning parameters is 3.63 ± 0.01 for Algorithm 1 and 4.77 ± 0.01 for Algorithm 2. These values are represented by the blue and red vertical lines.

curves correspond to the extended designs (Algorithm 2). The confidence bounds are plotted with faint lines in the according color and finally the mean values of the optimal tuning parameters are plotted with vertical lines. We find in both examples that the minimal prediction errors, that is, the minima of the red and blue curves, do not differ significantly. Additionally, in the first example, corresponding to highly correlated added variables (η = 0.001, see Figure 1), also the optimal tuning parameters do not differ significantly. However, in the second example (η = 0.1, see Figure 2), the optimal tuning parameter is considerably larger for the extended designs.

First Conclusions Our results indicate that tuning parameters proportional to √ n log p (cf.

[BRT09] and most other contributions on the subject) independent of the degree of correlations are not favorable. Indeed, for Algorithm 2, this would lead to a tuning parameter proportional to n log p 2 = √ 2n log p, whereas for Algorithm 1 to √ n log p. But regarding Figure 1, the two optimal tuning parameters are nearly equal and hence these choices are not favorable. In contrast, the results illustrate that the optimal tuning parameters depend strongly on the level of correlations: for Algorithm 2 (red, solid curves), the means of the optimal tuning parameters corresponding to the highly correlated case (3.63 ±0.01, see Figure 1) are be considerably smaller than the ones corresponding to the weakly correlated case (4.77 ± 0.01, see Figure 2).

Our results indicate additionally that the minimal mean prediction errors are comparable for all cases. This implies, that a suitable tuning parameters lead to good prediction even with additional impertinent parameters. We only give two examples here but made these observations for any values of n, p and s.

Theoretical Evidence

We provide in this section theoretical explanations for the above observations. For this, we first discuss results derived in [START_REF] Van De Geer | The lasso, correlated design, and improved oracle inequalities[END_REF]. They find that high correlations allow for small tuning parameters and that this can lead to bounds for Lasso prediction that are even more favorable than the fast rate bounds. Then, we introduce and apply new correlation measures that provide some insight for (in contrast to [START_REF] Van De Geer | The lasso, correlated design, and improved oracle inequalities[END_REF]) arbitrary degrees of correlations. For no correlations, in particular, these results simplify to the classical results.

Highly Correlated Designs

First results for the highly correlated case are derived in [START_REF] Van De Geer | The lasso, correlated design, and improved oracle inequalities[END_REF]. Crucial in their study is the treatment of the stochastic term with metric entropy.

The bound on the prediction error for Lasso reads as follows:

Lemma 3.1. [vdGL12, Theorem 4.1 & Corollary 4.2] On the set

T α := sup β 2σ|ǫ T Xβ| Xβ 1-α 2 β α 1 ≤ λ we have for λ = (2 λn α-1 ) 2 1+α β 0 α-1 1+α 1 and 0 < α < 1 X( β -β 0 ) 2 2 ≤ 21 2 2 λn α-1 2 1+α β 0 2α 1+α 1 .
We show in the following that for high correlations the stochastic term T α has a high probability even for small α and thus favorable bounds are obtained. The parameter α can be thought of as a measure of correlations: α small corresponds to high correlations, α large corresponds to small correlations.

The stochastic term T α is estimated using metric entropy. We recall, that the covering numbers N(δ, F , d) measure the complexity of a set F with respect to a metric d and a radius δ. Precisely, N(δ, F , d) is the minimal number of balls of radius δ with respect to the metric d needed to cover F . The entropy numbers are then defined as H(δ, F , d) := log N(δ, F , d). In this framework, we say that the design is highly correlated if the covering numbers N(δ, sconv{X (1) , ..., X (p) }, • 2 ) (or the corresponding entropy numbers) increase only mildly with 1/δ, where {X (1) , ..., X (p) } are the columns of the design matrix and sconv denotes the symmetric convex hull. This is specified in the following lemma: Lemma 3.2. [vdGL12, Corollary 5.2] Let 0 < α < 1 be fixed. Then, assuming

log 1 + N( √ nδ, sconv{X (1) , ..., X (p) }, • 2 ) ≤ A δ 2α , 0 < δ ≤ 1, (4)
there exists a value C(α, A) depending on α and A only such that for all κ > 0 and for λ = σC(α, A) n 2-α log(2/κ), the following bound is valid:

P(T α ) ≥ 1 -κ.
We observe indeed that the smaller α is, the higher the correlations are. We also mention that Assumption (4) only applies to highly correlated designs. An example is given in [START_REF] Van De Geer | The lasso, correlated design, and improved oracle inequalities[END_REF]: the assumption is met if the eigenvalues of the Gram matrix

X T X n
decrease sufficiently fast.

Lemma 3.1 and Lemma 3.2 can now be combined to the following bound for the prediction error of the Lasso: Theorem 3.1. With the choice of λ as in Lemma 3.1 and under the assumptions of Lemma 3.2, it holds that

X( β -β 0 ) 2 2 ≤ 21 2 σC(α, A) n α log(2/κ) 2 1+α β 0 2α 1+α
1 with probability at least 1 -κ.

High correlations allow for small values of α and lead therefore to favorable bounds. For α sufficiently small, these bounds may even outmatch the classical fast rate bounds. For moderate or weak correlations, however, Assumption (4) is not met and therefore the above lemma does not apply.

Arbitrary Designs

In this section, we introduce bounds that apply to any degree of correlations. The correlations are in particular allowed to be moderate or small and the bounds simplify to the classical results for weakly correlated designs. We first introduce two measure for correlations that are then used to bound the stochastic term. These results are then combined with the classical slow rate bound to obtain a new bound for Lasso prediction. We finally give some simple examples.

Two Measures for Correlations

We introduce two numbers that measure the correlations in the design. For this, we first define the correlation function K : R + 0 → N as

K(x) := min{l ∈ N : ∃x (1) , . . . , x (l) ∈ √ nS n-1 , (5)
X (m) ∈ (1 + x) sconv{x (1) , . . . , x (l) } ∀1 ≤ m ≤ p},
where S n-1 denotes the unit sphere in R n . We observe that K(x) ≤ p for all x ∈ R + 0 and that K is a decreasing function of x. A measure for correlations should, as the metric entropy above, measure how close to one another the columns of the design matrix are. Indeed, for moderate x, K(x) ≈ p for uncorrelated designs, whereas K(x) may be considerably smaller for correlated designs. This information is concentrated in the correlation factors that we define as

K κ := inf x∈R + 0 (1 + x) log(2K(x)/κ) log(2p/κ) ,
κ ∈ (0, 1], and as

F := inf x∈R + 0 (1 + x) log(1 + K(x)) log(1 + p) .
Since K(0) ≤ p, it holds that F, K κ ∈ (0, 1]. We also note that similar quantities could can be defined for p = ∞ (removing the normalization). In any case, large F and K κ correspond to uncorrelated designs, whereas small F and K κ correspond to correlated designs.

Control of the Stochastic Term

We now show that small correlation factors allow for small tuning parameters and thus lead to favorable bounds for Lasso prediction. Crucial in our analysis is again the treatment of a stochastic term similar to the one above.

We prove the following bound in the appendix:

Theorem 3.2. With the definitions above, T as defined in Section 2 and for all κ > 0 and λ ≥ λ κ := K κ 2σ 2n log(2p/κ), it holds that

P (T ) ≥ 1 -κ.
Additionally, independently of the choice of λ,

E sup β 1 ≤M σ | ǫ T Xβ | ≤ F σ 8n log(1 + p) 3 M.
For small correlation factors K κ and F , the minimal tuning parameters λ κ and the expectation of the stochastic term are small. For K κ → 1, the minimal tuning parameters simplify to 2σ 2n log(2p/κ) (cf. [BRT09]). Similarly, for F → 1, the expectation of the stochastic term simplifies to σ 8n log(1+p) 3

M.

Together with the slow rate bound introduced in Section 2, this permits the following bound: Corollary 3.1. For λ ≥ λ κ it holds that X( β -β 0 ) 2 2 ≤ 2λ β 0 1 with probability at least 1 -κ.

Our contribution to this result concerns the tuning parameters: for the classical value λ = 2σ 2n log(2p/κ), the bound simplifies to the classical slow rate bounds. Correlations, however, allow for smaller λ and thus lead to more favorable bounds.

Let us have a look at some general aspects of the above bound. Most importantly, Corollary 3.1 applies to any degree of correlations. In contrast, Theorem 3.1 only applies for highly correlated designs, and the classical fast rate bounds only apply for weakly correlated designs. We also observe that the sparsity index s, which appears in the classical fast rate bounds, does not appear in Corollary 3.1. Hence, Corollary 3.1 can be useful even if the true regression vector β 0 is only approximately sparse. However, for very large β 0 1 , the bound is unfavorable. An example in this context can be found in [CP09, Section 2.2]. (However, we do not agree with their conclusions corresponding to this example. We stress, in contrast, that correlation are not problematic in general.) We finally refer to Remark 2.1 and to [START_REF] Massart | The Lasso as an ℓ 1 -ball model selection procedure[END_REF] for some additional considerations on the optimality of this kind of bounds.

Remark 3.1 (Weakly Correlated Designs). Corollary 3.1 holds for any degree of correlations because it follows directly from the classical slow rate bound and the properties of the refined tuning parameter. However, for weakly correlated designs, we can also make use of the refined tuning parameter to improve the classical fast rate bound by a factor K 2 κ . Usually, the gain is moderate, but in special cases such as sparse designs (see Example 3.2), the gain can be large. Finally, we note that optimal rates, that is, lower bounds for the prediction error X( β -β 0 ) 2 2 , are available for weakly correlated designs. The optimal rate is then s log(1 + p s ) and is deduced from Fano's Lemma (see [BTW07a, Theorem 5.1, Lemma A.1]). For similar results in matrix regression, we refer to [START_REF] Koltchinskii | Nuclear norm penalization and optimal rates for noisy low rank matrix completion[END_REF], where the assumptions on the design needed for such lower bounds are explicitly stated (see, in particular, their Assumption 2).

Examples In this final section, we illustrate some properties of the correlation numbers K κ and F for various settings. We consider, in particular, design matrices with different geometric properties and with different degrees of correlations.

֒→ Example 3.1 (Low Dimensional Design). Let dim span{X (1) , ..., X (p) } ≤ W . Then, X (j) ∈ √ W sconv{x 1 , ..., x W } for all 1 ≤ j ≤ p
for properly chosen x 1 , ..., x W ∈ √ nS n-1 (for example orthogonal vectors in a suitable subspace). Hence,

K κ ≤ (1 + √ W ) log(2W/κ) log(2p/κ) and F ≤ (1 + √ W ) log(1+W ) log(1+p) .
֒→ Example 3.2 (Sparse Design). Let the number of non-zero coefficients in X (j) be such that X (j) 0 ≤ d for all 1 ≤ j ≤ p, with d ≤ n. Then,

X (j) ∈ √ d sconv{x 1 , ..., x n } for all 1 ≤ j ≤ p for properly chosen x 1 , ..., x n ∈ √ nS n-1 (for example orthogonal vectors in R n ). Hence, K κ ≤ √ d log(2n/κ) log(2p/κ) and F ≤ √ d log(1+n) log(1+p) .
֒→ Example 3.3 (Weakly Correlated Design). We consider a weakly correlated design with p ≫ n. It turns out that the results in Section 3.2.1 lead to a tuning parameter λ ∼ n . This implies in particular that the results in Section 3.2.1, unlike the results in this section, do not simplify to the classical results for uncorrelated designs. We give a sketch of the proof in the Appendix.

֒→ Example 3.4 (Highly Correlated Design). We illustrate with an example that highly correlated designs indeed involve small correlation factors K κ . This implies, in accordance to the simulation results, that tuning parameters depending on the design rather than tuning parameters depending only on σ, n, and p should be applied. Let us describe how the design matrix X is chosen: Fix an arbitrary n dimensional vector X (1) of Euclidian norm √ n. Then, we add p -1 columns sampled according to X (1) + νN where ν is a positive (but small) constant and N is a standard normally distributed random vector. We finally normalize the resulting matrix.

For simplicity, we do not give explicit constants but present a coarse asymptotic result that highlights the main ideas. For this, we consider the number of variables p as a function of the number of observations n and impose the usual restrictions in high dimensional regression, that is, p ∼ n w for a w ∈ N. If the correlations are sufficiently large, more precisely ν ≤ 1 √ 8n , it holds that

K κ → 1 √ w as n → ∞
with probability tending to 1. This implies especially that the standard tuning parameters ∼ σ √ n log p can be considerably too large if ν is large. The proof of the above statement is given in the Appendix. ֒→ Example 3.5 (Equal Columns). Let the cardinality of the set |{X (j) : 1 ≤ j ≤ p}| = v. Then, K κ ≤ log(2v/κ) log(2p/κ) and F ≤ log(1+v) log(1+p) .

Experimental Study

We consider Algorithm 1 with different sets of parameters to make statements about the influence of the single parameters on Lasso prediction. In particular, we are interested in the influence of the correlations ρ.

Results

We collect in Table 1 the means of the optimal tuning parameters λ min and the means of the minimal prediction errors P E min for 1000 iterations and different parameter sets. Let us first highlight the two most important observations: first, correlations (ρ large) lead to small tuning parameters. Second, correlations do not necessarily lead to high prediction errors. In contrast, the prediction errors are mostly smaller for the correlated settings.

Let us now make some other observations. First, we find that the optimal tuning 2.90 ± 0.03 6.65 ± 0.06 20 40 4 3 0.99 2.42 ± 0.15 6.16 ± 0.17 0.9 5.33 ± 0.13 6.47 ± 0.14 0 12.33 ± 0.10 3.80 ± 0.03 parameters do not increase considerably when the number of observations n is increased. In contrast, the means of the minimal prediction errors decrease for the correlated case as expected, whereas this is not true for the uncorrelated case. Second, increasing the number of variables p leads to increasing optimal tuning parameters as expected (interestingly by factors close to log 400 log 40 , cf. Corrollary 3.1). The means of the minimal prediction errors do, surprisingly, not increase considerably. Third, as expected, increasing the sparsity s does not considerably influence the optimal tuning parameters but leads to increasing means of the minimal prediction errors. Forth, for σ = 3 both the optimal tuning parameter as well as the mean of the minimal prediction error increase approximately by a factor 3. The mean of the minimal prediction errors for σ = 3 and ρ = 0 is an exception and remains unclear. We finally mention that we obtained analogeous results for many other values of β 0 and sets of parameters.

Conclusions

The experiments illustrate the good performance of the Lasso estimator for prediction even for highly correlated designs. Crucial is the choice of the tuning parameters: we found that the optimal tuning parameters depend highly on the design. This implies in particular that choosing λ proportional to √ n log p independent of the design is not favorable.

Discussion

Our study suggests that correlations in the design matrix are not problematic for Lasso prediction. However, the tuning parameter has to be chosen suitable to the correlations. Both, the theoretical results and the simulations strongly indicate that the larger the correlations are, the smaller the optimal tuning parameter is. This implies in particular, that the tuning parameter should not be chosen only as a function of the number of observations, the number of parameters and the variance. The precise dependence of the optimal tuning parameter on the correlations is not known, but we expect that cross validation provides a suitable choice in many applications.

Combining these two arguments and using the triangular inequality, we get

X( β -β 0 ) 2 2 ≤ λ ( β -β 0 ) J 0 1 + λ βJ c 0 1 + λ β 0 1 -λ β 1 ≤ 2λ ( β -β 0 ) J 0 1 .
On the other hand, using the triangular inequality, we also get

X( β -β 0 ) 2 2 ≤ λ β -β 0 1 + λ β 0 1 -λ β 1 ≤ 2λ β 0 1 .
This completes the proof.

Proof of Theorem 3.2. We first show that for a fixed x ∈ R + 0 , the parameter space {β ∈ R p : β 1 ≤ M} can be replaced by the K(x) dimensional parameter space {β ∈ R K(x) : β 1 ≤ (1 + x)M}. Then, we bound the stochastic term in expectation and probability and eventually take the infimum over x ∈ R + 0 to derive the desired inequalities. As a start, we assume without loss of generality σ = 1, we fix x ∈ R + 0 and set K := K(x). Then, according to the definition of the correlation function (5), there exist vectors x (1) , ..., x (K) ∈ √ nS n-1 and numbers {κ j (m) : 1 ≤ j ≤ K} for all These two results imply sup

1 ≤ m ≤ p such that X (m) = K j=1 κ j (m)x (j) and K j=1 |κ j (m)| ≤ (1 + x). Thus, (
β 1 ≤M | ǫ T Xβ |≤ sup β 1 ≤(1+x)M | ǫ T X β |,
where β ∈ R K and X := (x (1) , ..., x (K) ). That is, we can replace the p dimensional parameter space by a K dimensional parameter space at the price of an additional factor 1 + x.

We now bound the stochastic term in expectation. First, we obtain by Cauchy-Schwarz's Inequality E sup

β 1 ≤(1+x)M | ǫ T X β | = E sup β 1 ≤(1+x)M | n i=1 K j=1 ǫ i X (j) i β j | ≤ E sup β 1 ≤(1+x)M β 1 max 1≤j≤K | ǫ T X (j) | = (1 + x)ME max 1≤j≤K | ǫ T X (j) | .
Next (cf. the proof of [vdGL11, Lemma 3]), we obtain for Ψ(x) := e x 2 -1

E max 1≤j≤K | ǫ T X (j) | ≤ Ψ -1 (K) max 1≤j≤K ǫ T X (j) Ψ ,
where • Ψ denotes the Orlicz norm with respect to the function Ψ (see [START_REF] Van De Geer | The bernstein-orlicz norm and deviation inequalities[END_REF] for a definition). Since ǫ T X (j) √ n is standard normally distributed, we obtain ǫ T X (j)

√ n Ψ = 8 3 (see for example [vdVW00, Page 100]). Moreover, one may check that Ψ -1 (y) = log(1 + y). Consequently,

E sup β 1 ≤(1+x)M | ǫ T X β | ≤ (1 + x)M 8n log(1 + K) 3
One can then derive the second assertion of the theorem by taking the infimum over x ∈ R + 0 . As a next step, we deduce similarly as above (compare also to [BRT09])

P sup β 1 ≤M 2 | ǫ T Xβ |≥ λM ≤P sup β 1 ≤1 2 | ǫ T X β |≥ λ 1 + x ≤K max 1≤j≤K P ǫ T X (j) √ n ≥ λ 2(1 + x) √ n =KP |γ| ≥ λ 2(1 + x) √ n ,
where γ is a standard normally distributed random variable. Setting λ := 2(1 +

x) 2n log(2K/κ), we obtain

P sup β 1 ≤M 2 | ǫ T Xβ |≥ λM ≤ 2K exp - λ 2 8(1 + x) 2 n = 2K exp (-log(2K/κ)) = κ.
The first assertion can finally be derived taking the infimum over x ∈ R + 0 and applying monotonous convergence.

Proposition 1. It holds for all 0 < δ ≤ 1 1 δ n ≤ N(δ, {x ∈ R n : x 2 ≤ 1}, • 2 ) ≤ 3 δ n .
Proof. If a collection of balls covers the set {x ∈ R n : x 2 ≤ 1}, the volume covered by these balls is at least as large as the volume of {x ∈ R n : x 2 ≤ 1}. Thus,

N(δ, {x ∈ R n : x 2 ≤ 1}, • 2 ) • δ n ≥ 1,
and the left inequality follows. The right inequality can be deduced similarly.

Proof Sketch for Example 3.3. For weakly correlated designs and for p very large, it holds that

log 1 + N( √ nδ, sconv{X (1) , ..., X (p) }, • 2 ) ≈ log (N(δ, {x ∈ R n : x 2 ≤ 1}, • 2 )) .
We can then apply Proposition 1 to deduce

log 1 + N( √ nδ, sconv{X (1) , ..., X (p) }, • 2 ) ≈ n log(1/δ).
Hence, Lemma 3.2 requires that A fulfills for all 0 < δ ≤ 1 (in fact, it is sufficient for

1 √ n
δ ≤ 1, see the proofs in [START_REF] Lederer | Bounds for rademacher processes via chaining[END_REF] for this refinement)

A δ 2α ≈ n log(1/δ) and thus

A α ≈ δ α n log(1/δ).

We can maximize the right hand side over 0 < δ ≤ 1 and plug the result into [vdGL12, Corollary 5.2] to deduce the result.

Proof for Example 3.4. The realizations of the random vectors X (1) + νN are clustered with high probability and can thus be described (in the sense of the definition of K κ ) by much less than p vectors. To see this, consider X (1) and a rotation R ∈ SO(n) such that RX (1) = ( √ n, 0, . . . , 0) T .

Moreover, let N be a standard normally distributed random vector in R n , n ≥ 2.

The measure corresponding to the random vector N is denoted by P. We then obtain with the triangle inequality and the condition on ν

P 1 - √ n RX (1) + νRN 2 ≥ 1 =P RX (1) + νRN 2 √ n ≤ 1 2 ≤P X (1) 2 -ν N 2 √ n ≤ 1 2 =P ν N 2 √ n ≥ 1 2 ≤P N 2 ≥ √ n 2ν ≤P N 2 ≥ √ 2n .
Now, we can bound the l 1 distance of the vector √ n(RX (1) +νRN )

RX (1) +νRN 2 generated according to our setting to the vector RX (1) : P RX (1) -√ n(RX

Finally, define for any 0 ≤ z ≤ 2 √ n the set C z := {( √ n, z, 0, . . . , 0) T , ( √ n, -z, 0, . . . , 0) T , ( √ n, 0, z, 0, . . . , 0) T , ( √ n, 0, -z, 0, . . . , 0) T , . . . } ⊂ R n These vectors will play the role of the vectors x (1) , ..., x (l) in Definition (5). It holds that card(C z ) = 2(n -1)

and

C z ⊂ √ n + z 2 S n-1 . (8) Moreover, {x ∈ R n : x -RX (1) 1 ≤ z} ∩ √ nS n-1 ⊂ sconv(C z ). (9)
Inequality (6) and Inclusion (9) imply that with probability at least 1 -2 exp -n 8 √ n(X (1) + νN)

RX (1) + νRN 2 ∈ sconv(R T C 2 √ n ).
Thus, using Equality (7) and Inclusion (8), we have K( √ 5 -1) ≤ 2(n -1) with probability at least 1 -2(p -1) exp -n 8 .

[BC11] A. Belloni and V. Chernozhukov. High dimensional sparse econometric models: An introduction. In P. Alquier, E. Gautier, and G. Stoltz, editors, Inverse Problems and High-Dimensional Estimation. Springer (Lecture Notes in Statistics), 2011.

[BJMO11] F. Bach, R. Jenatton, J. Mairal, and G. Obozinski. Convex optimization with sparsity-inducing norms. In S. Sra, S. Nowozin, S. J. Wright., editors, Optimization for Machine Learning, MIT Press, 2011.

[BRT09] P. Bickel, Y. Ritov, and A. Tsybakov. Simultaneous analysis of lasso and Dantzig selector. Ann. Statist., 37(4):1705-1732, 2009.

Figure 1 :

 1 Figure 1: We plot the mean values P E of the prediction errors A β -Aβ 0 2

Figure 2 :

 2 Figure 2: We plot the mean values P E of the prediction errors A β -Aβ 0 2

  j (m)| ≤ (1 + x) β 1 .

Table 1 :

 1 The means of the optimal tuning parameters λ min and the means of the minimal prediction errors P E min calculated according to Algorithm 1 with 1000 iterations and for different sets of parameters.

	n	p	s σ	ρ	λ min	P E min
	20 40 4 1 50 40 4 1 20 400 4 1 20 40 10 1	0.99 0.69 ± 0.03 1.77 ± 0.05 0.9 1.58 ± 0.03 2.37 ± 0.04 0 3.60 ± 0.03 3.17 ± 0.03 0.99 0.67 ± 0.03 1.29 ± 0.04 0.9 1.71 ± 0.03 1.85 ± 0.03 0 3.91 ± 0.03 3.48 ± 0.02 0.99 0.97 ± 0.03 1.75 ± 0.05 0.9 2.11 ± 0.04 2.58 ± 0.04 0 4.82 ± 0.03 3.34 ± 0.03 0.99 0.59 ± 0.03 6.50 ± 0.22 0.9 1.46 ± 0.03 7.97 ± 0.19 0

Acknowledgments We thank Sara van de Geer for the great support. Moreover, we thank Arnak Dalalyan, who read a draft of this paper carefully and gave valuable and insightful comments. Finally, we thank the reviewers for their helpful suggestions and comments.

* JCL acknowledges partial financial support as member of the German-Swiss Research Group FOR916 (Statistical Regularization and Qualitative Constraints) with grant number 20PA20E-134495/1.

Appendix

Proof of Remark 2.1. By the definition of the Lasso estimator, we have

Next, on the set T , we have 2σ|ǫ T X( β -β 0 )| ≤ λ β -β 0 1 . Additionally, by the definition of J 0 , β -β 0 1 = ( β -β 0 ) J 0 1 + βJ c 0 1 .


Konstantin K Starkov
email: k.starkov@tue.nl

Alexander Y Pogromsky

J E Rooda

Taylor & Francis

I T Consultant

Performance analysis for tandem manufacturing lines under variable structure production control method

Keywords: MANUFACTURING CONTROL SYSTEMS, PRODUCTION manufacturing systems, tracking systems, discrete time systems, boundary conditions, variable structure control

Performance analysis for tandem manufacturing lines under variable structure production control method.

Introduction

In the modern market keeping high competition in brands and varieties in type of products is the way for survival of manufacturing industries. Therefore production control methods with capabilities of quick responses to rapid changes in the demand and efficient distribution of the raw material throughout the network are of a big importance among leading manufacturers. Nowadays, the production control problem has been widely studied and a lot of valuable approaches including queuing theory, Petri nets, dynamic programming, linear programming, hybrid systems were proposed and some of them are implemented. Though up to this moment many methods have been developed, the factory performance remains unpredictable. [START_REF] Subramaniam | A WIP control policy for tandem lines[END_REF], [START_REF] Savkin | Optimal distributed real-time scheduling of flexible manufacturing networks modeled as hybrid dynamical systems[END_REF], [START_REF] Nilakantan | Enhancing supply chain performance with improved order-control policies[END_REF] and references there in), can be summarized as follows. The proposed tandem production model is considered in discrete time. The production speed of each machine is defined as deterministic with bounded perturbations. The future production demand is assumed to be unknown and with bounded fluctuations. As a result, strict, so called "worst case" bounds on a content of intermediate buffers and production tracking errors for one machine and for a unidirectional tandem production line are obtained as a starting point of our research. In order to monitor a current market situation, and the current inventory level of given network the proposed strategy can be easily combined with integrated computer-based systems (such as previously mentioned in this section).

In particular, this paper contains an explicit description of flow models (see, e.g. [START_REF] Dallery | Manufacturing flow line systems: a review of models and analytical results[END_REF], [START_REF] Alvarez-Vargas | A study of the continuous flow model of production lines with unreliable machines and finite buffers[END_REF], [START_REF] Pogromsky | Decentralized feedback control of manufacturing machines[END_REF] and references therein) for one machine and a line of manufacturing machines. Here a variable structure control (see, e.g. [START_REF] Utkin | Variable Structure Systems: Present and Future[END_REF], [START_REF] Khalil | Nonlinear Systems[END_REF]) is introduced as a control technique in order to give a solution to the demand tracking problem. Special attention is paid to the constraints presented in the network, such as capacity and buffer limitations. Each machine in the network can produce a restricted number of products in a fixed period of time, known as the capacity constraint. The content of the buffer between two machines is given by the difference between the total number of products produced by the upstream machine and the total number of products produced by the downstream machine. Considering that a manufacturing line has a unidirectional product flow implies that the buffer content can never be negative, e.g. the downstream machine cannot produce more than the upstream one.

The paper is organized as follows. First, in Section 2 the flow model of one manufacturing machine with variable structure control is presented, where the analysis of production error trajectories is developed. Then the flow model of a manufacturing line with variable structure control is analyzed in Section 3. Here necessary conditions are derived to guarantee the uniform ultimate boundedness of the production error trajectories of each machine. Performance and robustness issues of the closed-loop flow models are illustrated in numerical simulations in Section 4. Finally, Section 5 contains conclusions and implementation in a real environment.

One Machine

Flow Model

In discrete time a cumulative number of produced products in time k for a simple manufacturing machine can be described as a sum of its production rates at each time step till time k. Thus the flow model of one manufacturing machine in discrete time is defined as

y(k + 1) = y(k) + u(k) + f (k), (1)
where y(k) ∈ R is the cumulative output of the machine in time k, u(k) ∈ R is the control signal, and f (k) ∈ R is an unknown external disturbance. Under the assumption that there is always sufficient quantity of the raw material to feed the machine, the control aim is to track the non-decreasing cumulative production demand. We define the production Taylor & Francis and I.T. Consultant demand by using y d (k) ∈ R given by

y d (k) = y d0 + v d k + ϕ(k), (2)
where y d0 is a positive constant that represents the initial production demand, v d is a positive constant that defines the average desired demand rate, and ϕ(k) ∈ R is the bounded fluctuation that is imposed on the linear demand v d k. In order to give a solution to this tracking problem the following control input

u(k) = µsign + (ε(k)) (3) is considered.
Here µ is a positive constant that represents the processing speed of the machine, step function sign

+ (ε(k)) =(1, if ε(k) > 0|0, otherwise)
, and ε(k) ∈ R is the output production error with respect to the demand. This production error is given by ε

(k) = y d (k) -y(k)
, where ε(k + 1))ε(k) along the solutions of ε(k) is given by:

ε(k + 1) -ε(k) = v d + ∆ϕ(k) -µsign + (ε(k)) -f (k), (4)
with ∆ϕ(k) = ϕ(k + 1)ϕ(k). It is followed from (4) that in order to guarantee a proper demand trajectory tracking the product demand cannot be higher than the machine processing speed, which in this case is µ lots per time unit. Thus, let us consider that all machine perturbations W (k) = ∆ϕ(k)f (k) from ( 4) are bounded by

α 1 < W (k) < α 2 , ∀k ∈ N, (5)
where α 1 , α 2 are some constants that satisfy

α 2 < µ -v d , (6)
α 1 > -v d . (7)
By ( 6) and ( 7) we state that the machine can never produce products faster than its maximal speed and that considering the presence of perturbations bounded by (α 1 , α 2 ) the demand rate can only de positive, respectively. Thus, from (5), (6), and (7) the following condition (also known as capacity condition) holds

0 < v d + W (k) < µ. (8)
Basically, in order to follow the product demand, variable structure controller µsign + (ε(k)) is included in the flow model of one machine. The production error of a single machine is defined as the difference between the cumulative demand and the cumulative number of products produced up to this moment.

Production error analysis

In this section we analyse the tracking error trajectories behavior of the flow model (4). Theorem 2.1 : Assume that the discrete time system defined by (1) with control input (3) satisfies condition (8). Then all solutions of (4) are uniformly ultimately bounded by

lim sup k→∞ ε(k) ≤ v d + α 2 , (9) lim inf k→∞ ε(k) ≥ v d + α 1 -µ. ( 10
)
Proof : The proof of Theorem 2.1 is provided in the first part of the proof of Theorem 3.4. The flow model of a manufacturing line is presented in this section. Figure 1 presents a schematics of a line of N manufacturing machines with machines M j , buffers B j , and infinite product supply. Here the control strategy for one machine is modified with respect to the number of buffers and machines present in the line. A new limitation such as desired buffer content is considered in the model.

A Line of Machines

The flow model of the manufacturing line is defined as

y 1 (k + 1) = y 1 (k) + β 1 (k), (11)
y j (k + 1) = y j (k) + (β j (k))sign Buff (w j (k) -β j (k)), ( 12
) ∀j = 2, . . . , N, where y j (k) is the cumulative output of machine M j in time k, w j (k) = y j-1 (k) -y j (k) is the buffer content of buffer B j , β j (k) = u j (k) + f j (k), ∀j = 1, .
. . , N , f j is the external disturbance affecting machine M j (e.g. production speed variations, undesired delay), u j is the control input of machine M j and sign Buff (w

j (k)-β j (k)) = (1, if w j (k)-β j (k) ≥ 0|0, otherwise).
In order to give a solution to the demand tracking problem we consider the following control inputs:

u j (k) = µ j sign + (ε j+1 (k) + (w dj+1 -w j+1 (k))) (13) ∀j = 1, . . . , N -1, u N (k) = µ N sign + (y d (k) -y N (k)), ( 14
)
where µ j is the processing speed of machine j, w dj+1 is the desirable buffer level of buffer B j+1 and ε j+1 is the tracking error of machine M j+1 . The tracking error of each machine Taylor & Francis and I.T. Consultant is given by:

ε j (k) = ε j+1 (k) + (w dj+1 -w j+1 (k)) (15) ∀j = 1, . . . , N -2 ε N -1 (k) = ε N (k) + (w dN -w N (k)), ( 16
) ε N (k) = y d (k) -y N (k). ( 17
)
It follows from ( 17) that the error of machine M N is defined exactly as for the single machine case. The buffer restriction, as seen from ( 12), is the only difference in the flow model of machine M N with the flow model of (1). For ( 15), ( 16) new considerations are applied for the tracking error of each machine M j , where j = 1, ..., N -1. Thus, tracking error ε j (k) depends on the number of produced products y j (k) with respect to current demand y d (k) and desired buffer content w dj+1 of each downstream buffer. This means that every upstream machine needs to supply w dj+1 lots more than the downstream one.

Constant parameter w d is introduced in order to prevent downstream machines from lot starvation, e.g. in case of a sudden growth of the product demand. Basically, model ( 11), ( 12) describes the product flow through the line of N manufacturing machines (Fig. 2). The first machine described by ( 11) is considered to have always an access to the raw material and the raw material is always sufficient. The administration of this raw material to machine M 1 is decided by the control input ( 13). Here we consider that our control input is acting as an authorizing switch, which turns on M 1 if its tracking error ( 15) is positive and turns M 1 off if its tracking error is negative or zero. Tracking error ε 1 (k), see ( 15), consists of the difference between what is done (w 2 (k)), what have to be done (ε 2 (k)) and what have to be always in the buffer (w d2 ). It can be seen from ( 15), ( 16) and ( 13) that the same tracking error and control logics were applied for the rest of the machines till machine M N -1 . As for the last machine in the line, which is machine M N (Fig. 2), the control action is still based on the authorizing switch. However, the difference consists in the logic that triggers this switch. We expect that on the output of machine M N the cumulative product demand is followed by cumulative production of this machine. The control switch activates or deactivates the machine based directly on the production demand status (17). This control logic is the same as in the one machine case presented in the previous section. The substantial difference in models of the rest of machines from machine one can be appreciated through the general flow model ( 12). Here for each machine in the line we introduce an extra restriction on the buffer content of each upstream buffer. Function sign Buff (w j (k)β j (k)) is acting as an extra authorization together with the control input. Thus, any machine M j , with j = 2, . . . , N , is activated only if two authorizations are given. The first authorization comes from control input (u j (k)) of the machine, which is based on the current tracking error status of this machine (ε j (k)). The second authorization comes from the buffer content restriction which is granted if the buffer contains at least the minimal number of products required (β j (k)) in order for the machine M j to start its work. It is also important to take into account that the control actions are decentralized throughout the network. In other words the control action of each machine in the line depends only on the tracking error of its neighboring downstream machine (except for machine M N , which depends directly on cumulative demand input) and the current buffer content of its upstream buffer (Fig. 2). This gives our flow model an extra robustness with respect to the undesired events such as temporal machine setup or breakdown.

For further analysis, let us rewrite flow model ( 11), ( 12) in a closed-loop with ( 13), ( 14) as

∆ε 1 (k) = v d + ∆ϕ(k) -f 1 (k) -µ 1 sign + (ε 1 (k)), ( 18
) ∆ε j (k) = v d + ∆ϕ(k) -f j (k) -µ j sign + (ε j (k))sign Buff (w j (k) -β j (k)), (19)
where ∆ε

j (k) = ε j (k + 1) -ε j (k) and W j (k) = ∆ϕ(k) -f j (k), ∀j = 1, . . . , N .
Here we consider that system (18), ( 19) satisfies the following assumptions.

Assumption 3.1: (Boundedness of perturbations) There are constants α 1 , α 2 and α 3 such that W j (k) = ∆ϕ(k)f j (k) satisfies

α 1 < W j (k) < α 2 , ∀k ∈ N, j = 1, . . . , N, (20)
and f j (k) satisfies

f j (k) ≤ α 3 , ∀k ∈ N, j = 1, . . . , N. (21)
Assumption 3.2: (Capacity condition) Constants α 1 , α 2 satisfy the following inequalities

α 2 < µ j -v d , ∀j = 1, . . . , N, (22)
α 1 > -v d . (23)
Thus, from (20), ( 22), and ( 23) the following condition holds

0 < v d + W j (k) < µ j , ∀j = 1, . . . , N. (24)
It is important to notice that each M j machine in the line has a processing speed of µ j lots per time unit, which can differ from the rest of the machines, and the buffer condition is considered as Taylor & Francis and I.T. Consultant

w j (k) ≥ β j (k), ∀j = 2, . . . , N. (25)
Thus, from ( 15), ( 16) and ( 25) the following tracking error condition holds

ε j (k) ≥ β j (k) -w dj + ε j-1 (k), ∀j = 2, . . . , N,
where w dj satisfy the following assumption.

Assumption 3.3: (Desired buffer content condition) The constants w dj comply with the following inequality

w dj ≥ µ j + µ j-1 + α 3 + α 2 -α 1 ,
from where it follows that

w dj ≥ β j (k) + µ j-1 + α 2 -α 1 , k ∈ N, j = 2, . . . , N. (26)
If the condition ( 25) is not satisfied then

ε j-1 (k) (15,16,26) > µ j-1 + α 2 -α 1 + ε j (k), ∀j = 2, . . . , N. (27)

Production Error Analysis

In this section we present results respecting the production error trajectories behavior of the flow model ( 18), (19).

Theorem 3.4 : Assume that the discrete time system defined by ( 18), ( 19) satisfies Assumptions 3.1, 3.2, and 3.3 . Then all solutions of ( 18) and ( 19) are uniformly ultimately bounded by

lim sup k→∞ ε j (k) ≤ v d + α 2 , ( 28
) lim inf k→∞ ε j (k) ≥ v d + α 1 -µ j . ( 29
)
Proof : The proof of Theorem 3.4 is provided in Appendix A.

In consequence, for the buffer content w j (k) of each buffer B j defined by ( 15), ( 16), considering the obtained tracking error bounds (28), (29), it holds that lim sup

k→∞ w j (k) ≤ µ j-1 + α 2 -α 1 + w dj .
Now, in order to support the present development let us extend our analysis to simulation results.

One machine case

Simulation results for one manufacturing machine, driven by the variable structure regulator (3), are presented in this section. Here the processing speed was fixed to 5 lots per time unit.

Results from Figure 3 show that output y(k) follows product demand y d (k) with tracking error ε(k) bounded by 4 ≥ ε(k) ≥ 0 lots, which satisfies ( 9), (10). Here, initial conditions were set to y d0 = 20, y(0) = 0. Market fluctuations as well as external perturbations were set to the zero value. The tracking error from Figure 3 (b) is generated in this form due to the machine production speed limit of 5 lots per time unit with demand increment ∆y d (k) of 4 lots per time unit. This means that after the moment when the steady state is reached (20 time steps) demand signal y d (k) is 4 lots higher than output y(k). In this case the control responds with authorization for 5 lots, which are generated in the next time step, where the cumulative demand grows in another 4 lots. Thus, the current output is 3 lots lower than the current demand. Then for the next 4 time steps the machine keeps producing lots until the demand is reached again at time step 25 with y d (k) of 120 lots. Now during next time step the machine remains idle while the demand grows to 124 lots. Thus, the sequence from the 20th time step to the 25th is repeated. Finally we conclude that simulation results reflect the expected flow model behavior for all the product demands which are given in this section. These product demands were selected in order to test the model behavior inside the boundary of the given capacity condition (8).

Four machines case

M1 M2 M3 M4 B2 B3 B4 y1 y2 y3 y4
Figure 5. Schematics of a line of 4 manufacturing machines.

Simulation results for a line of 4 manufacturing machines (Fig. 5), driven by variable structure regulators (13), and ( 14) are presented in this section. For all the examples processing speed for each machine was set to µ j =< 8, 10, 7, 6 > (lots per time unit), with j = 1, .., 4, and the desired buffer content of each buffer was selected considering (26) as w dj =< 20, 18, 14 > (lots), with j = 2, .., 4.

The tracking error of each machine in the line is depicted in Figure 7. Here the initial conditions (y d0 , y 1 (0), y 2 (0), y 3 (0), y 4 (0)) were set to the zero value. After the first 24 time steps, as it is shown in Figures 6(a) and 7, the system reaches its steady state. Tracking errors are maintained inside [-2,5] lots for machine M 1 , [-4,1] lots for machine M 2 , [-1,5] lots for machine M 3 , and [0,5] lots for machine M 4 , which satisfy ( 28) and ( 29). From Figure 6(b) it can be observed that the inventory level of each buffer satisfies the upper bound restriction (30). Further, the output response and the tracking error of each machine are depicted in Figures 8 and9(a) for the same initial conditions as in the previous example. Here it is noticeable that after the first 7 time steps the output of machine M 4 reaches the current demand trajectory. Tracking errors are maintained inside [-6,4] lots for machine M 1 , [-8,4] lots for machine M 2 , [-5,4] lots for machine M 3 , and [-4,4] lots for machine M 4 , which satisfy (28), ( 29). The inventory level of each buffer is depicted in Figure 9(b). We indicate that the buffer content of each buffer satisfies the upper bound restriction (30), if given desired inventory level of each buffer satisfies (26).

Conclusion

The performances of a tandem production line under variable structure surplus control has been studied. Developed results prove the uniform ultimate boundedness for trajectories of each tracking error and provide the upper bound for each buffer content in a line of machines where each machine has a variable processing. Various simulation results are presented and discussed in order to illustrate and support analytical results. Simulation results reflect effectiveness and robustness of the flow models.

As can be seen from the model description, the control algorithm requires to have a knowledge of external and internal variables, which affect the production line functionality, as well as performance specifications of each machine. The external variable which is the current production demand yd(k) as well as internal variables which are buffer contents w j (k) of each buffer can be obtained from existing planning details of a given company. Performance specifications such as machine production speed can be deduced from the experimental data or obtained from the user specifications of a given machine. Thus having this knowledge the production method can be implemented in form of an algorithm as a manufacturing line control tool. Further, if the conditions from Theorem 3.4 are satisfied then the system's results on performance can be directly obtained from formulas ( 28), ( 29) and ( 30). This can also significantly reduce the computational time needed for an order release planning.

Furthermore, studies on variable structure control policy's application to re-entrant network, multiple part type production systems, and performance analysis with the presence of production delays and setup times as well as it comparison with another surplus based strategies will be pursued in our future research. IntJourProdResearchPa-per˙Dec2010*RevComments REFERENCES 15 which means that machine M 2 has sufficient material in its buffer B 2 in order to start working and machine M 1 always has an access to the infinite raw material supply. Thus these machines have an independent behavior and it will be sufficient to analyse the increment of only one of the functions V j,k in order to determine the behavior of ∆V 2M k . Let us assume that ε j (k) satisfies the following condition ) and in consequence from (A3) and (A4) it follows that η j,k = 1. Then, ∆V j,k along the solutions of ε j (k) is given by

F

ε j (k) > 0, ( A6
∆V j,k = max    -ε j (k) -W j (k) + α 1 , ε j (k) + W j (k) -α 2 -µ j , 0    Vj,k+1 + min    ε j (k) + µ j -v d -α 1 , -ε j (k) + v d + α 2 , 0    -Vj,k , ∀j = 1, 2.
From where with help of Assumptions 1 and 2 it can be easily deduced that

∆V j,k =      0 if ε j (k) ≤ v d + α 2 , -ε j (k) + v d + α 2 < 0 if v d + α 2 < ε j (k) ≤ µ j + α 2 -W j (k), -µ j + v d + W j (k) < 0 if ε j (k) > µ j + α 2 -W j (k). (A7)
Now, suppose that for ε j (k) the following condition holds ) and in consequence from (A3) and (A4) it yields that η j,k = 0. Then ∆V j,k along the solutions of ε j (k) is given by

ε j (k) ≤ 0, ( A8
∆V j,k = max    -ε j (k) -W j (k) + α 1 -µ j , ε j (k) + W j (k) -α 2 , 0    Vj,k+1 + min    ε j (k) + µ j -v d -α 1 , -ε j (k) + v d + α 2 , 0    -Vj,k , ∀j = 1, 2.
Here with help of Assumptions 1 and 2 it can be easily deduced that

∆V j,k =      0 if v d + α 1 -µ j ≤ ε j (k) ≤ 0, -v d -W j (k) < 0 if ε j (k) < -µ j + α 1 -W j (k), ε j (k) + µ j -v d -α 1 < 0 if -µ j + α 1 + W j (k) ≤ ε j (k) < v d + α 1 -µ j . (A9)
Summarizing, for conditions (A5), (A6), and (A8) from (A7) and (A9) it holds that if V j,k > 0 its increment ∆V j,k < 0. From the definition of min it yields that for ∆V 2M k given by (A2) the following inequality is satisfied

∆V 2M k ≤ V j,k+1 -V j,k ≤ 0, ( A10

where i = arg max j=1,2 {V j,k+1 }. Note that ∆V j,k = V j,k+1 -V j,k = 0 only if either first condition of (A7) or first condition of (A9) is satisfied and for all ε

j (k) ∈ [v d + α 1 - µ j , v d + α 2 ] it follows that ∆V 2M k < 0. Thus, in this case it holds that for V 2M k > 0 given by (A1) its increment ∆V 2M k < 0.
Case 2 (Insufficient Buffer Content) Let us assume that w 2 (k) satisfies the following inequality

w 2 (k) < β 2 (k), ∀k ∈ N, (A11) and ε 2 (k) satisfies ε 2 (k) ≤ 0, ∀k ∈ N. (A12)
Then from ( 27) it holds that

ε 1 (k) > µ 1 + α 2 -α 1 + ε 2 (k). ( A13
)
Here similarly to Case 1 the behavior of these two machines can be considered independently. Thus, for ε 1 (k) satisfying (A13) it holds that ∆V 1,k is given by (A7) or (A9) if ε 2 (k) < -µ 1α 2 + α 1 or ∆V 1,k is given by (A7) if -µ 1α 2 + α 1 ≤ ε 2 (k) ≤ 0. For ε 2 (k) satisfying (A12) the increment ∆V 2,k is given by (A9). In consequence for ∆V 2M k given by (A2) the inequality (A10) in this case is also satisfied. Now, let us assume that ε 2 (k) satisfies ε 2 (k) > 0, ∀k ∈ N.

(A14)

Then from (27) it holds that ε 1 (k) is given by (A13). In this case M 2 has a positive tracking error, but its buffer B 2 has insufficient raw material content (A11) in order to start working (η 2,k = 0). Machine M 1 has a positive error as well, but due to its infinite raw material supply connection it can immediately initiate its production process (η 1,k = 1). Thus, for (A14) and (A13) let us rewrite ∆V 2M k from (A2) as

∆V 2M k = max            -ε 1 (k) -W 1 (k) + α 1 , ε 1 (k) + W 1 (k) -α 2 -µ 1 , -ε 2 (k) -W 2 (k) + α 1 -µ 2 , ε 2 (k) + W 2 (k) -α 2 , 0            V 2M k+1 + min            ε 1 (k) + µ 1 -v d -α 1 , -ε 1 (k) + v d + α 2 , ε 2 (k) + µ 2 -v d -α 1 , -ε 2 (k) + v d + α 2 , 0            -V 2M k .(A15)
It follows from ( 20), (A14) and (A13) that ∆V 2M k from (A15) can be reduced to

∆V 2M k = max    ε 1 (k) + W 1 (k) -α 2 -µ 1 , ε 2 (k) + W 2 (k) -α 2 , 0    V 2M k+1 + min    -ε 1 (k) + v d + α 2 , -ε 2 (k) + v d + α 2 , 0    -V 2M k . (
V N M k = max {V 1 (ε 1 ), ..., V N (ε N )} .
Here the similar reasoning is followed as for the proof for 2 machines.

The analysis is subdivided into the same 2 cases. Case 1 (Sufficient Buffer Content) and the first part of Case 2 (w j (k) < β j (k) and ε j (k) ≤ 0, ∀j = 2, . . . , N ) are solved identically to the proof for the line of 2 machines. For the second part of Case 2 the proof mostly relies on the condition ( 27) and the assumption that machine M 1 has always an access to the infinite raw material supply. Due to the extensive technical details we omit the complete analysis for a line of N machines and restrict ourselves by only giving this general idea of the procedure.

Figure 1 .

 1 Figure 1. Schematics of a line of N manufacturing machines.

Figure 2 .

 2 Figure 2. Flow model diagram for a line of N manufacturing machines.

Figure 3 .Figure 4 .

 34 Figure 3. Demand vs Output (a) and Tracking Error (b), with v d = 4 and y d0 = 20.

Figure 4

 4 Figure 4(b) shows the tracking error and the output response of the machine to the nonlinear demand growth. Here product demand y d (k) = 35+ 2k + sin(5k). The resulting tracking error is bounded by 3 > ε(k) > -4 lots, which satisfies (9) and (10). It can be

Figure 6 .

 6 Figure 6. Outputs y j (k) vs. Demand yd(k) (a) and Buffer Content w j (k) (b), with v d = 5.

Figure 7 .Figure 8 .

 78 Figure 7. Tracking Error ε j (k), with v d = 5.

Figure 9 .

 9 Figure 9. Outputs y j (k) vs. Demand yd(k) (a) and Buffer Content w j (k) (b), with v d = 3 and ∆ϕ(k) = sin(50k).

Finally

  , presented simulation results on selected examples reflect the desired flow model behavior. All technical conditions proposed in this section corresponds to analytical results described in Section 3.

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.ukInternational Journal of Production Research

Acknowledgement

The research leading to these results has received partial funding from the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement no. INFSO-ICT-223844.

Appendix A. Proof of Theorem 3.4

Let us prove that Theorem 3.4 holds for a line of 2 manufacturing machines (j = 1, 2) defined by ( 18) and ( 19). With this goal, let us introduce the following Lyapunov function

where

Here for the sake of brevity

Thus, ∆V 2M k along the solutions of ε 1 (k) and ε 2 (k) is given by

where

Here for the sake of brevity we introduce η j,k as

In order to perform a more detailed analysis on ∆V 2M k , let us divide this proof into 2 cases.

Case 1 (Sufficient Buffer Content) Suppose that w 2 (k) satisfies the following inequality

International Journal of Production Research IntJourProdResearchPa-per˙Dec2010*RevComments REFERENCES 17 Now, let us prove that for ε 1 (k) given by (A13) inequality

is satisfied. Indeed, from condition (A13) it yields that

Thus, inequality (A17) is satisfied. Also, from (A18) it holds that

Now, considering (A17) and (A19) we can rewrite V 2M k+1 given by a first term of (A16) as

Let us prove that for ε 1 (k) given by (A13) inequality

is satisfied. Here from condition (A13) it yields that

Thus, inequality (A21) is satisfied. From inequalities (A13), ( 22) it follows that

From (A21), (A23) we can rewrite V 2M k given by a second term of (A16) as

Having V 2M k+1 given by (A20) and V 2M k given by (A24), we can finally reduce ∆V 2M k from (A16) to

Thus, for this Case it holds that for V 2M k > 0 given by (A1) its increment ∆V 2M k < 0. Summarizing for 2 cases, we have shown that for V 2M k > 0 given by (A1) its increment ∆V

Thus, lim sup k→∞ V 2M k = 0 which completes our proof. In this proof we have analyzed the increment of the proposed Lyapunov function by means of 2 cases. Now for a line of N manufacturing machines (j = 1, ..., N ) defined by


David Peidro

Josefa Mula
email: fmula@cigip.upv.es

Alemany Diaz

Francisco-Cruz Lario
email: fclario@cigip.upv.es

David ; Peidro

David Peidro

M M E Alemany

Dpto David Peidro

Fuzzy multi

Keywords: Master planning, supply chain, ceramic sector, fuzzy multi-objective model. Master planning, supply chain, ceramic sector, fuzzy multi-objective model

published or not. The documents may come

Fuzzy multi-objective optimisation for master planning in a ceramic supply chain 1 Introduction

A supply chain (SC) may be considered an integrated process in which a series of several organisations (e.g., suppliers, manufacturers, distributors, shops), work together to acquire raw materials and to convert these materials into finished goods to be distributed to shops [START_REF] Beamon | Supply chain design and analysis: Models and methods[END_REF]. Integrated SC management optimises the SC as a whole by simultaneously considering planning for all the links involved. Thus to obtain optimum solutions to minimise costs, several authors have studied SC production planning with mathematical programming models (for instance, von Lanzenauer and Pil-Glombik (2002); [START_REF] Kreipl | Planning and scheduling in supply chains: An overview of issues in practice[END_REF]; [START_REF] Chen | Multi-objective optimization of multi-echelon supply chain networks with uncertain product demand and prices[END_REF]; [START_REF] Park | An integrated approach for production and distribution planning in supply chain management[END_REF]; Peidro et al. (2009b); Mula et al. (2010b)). For an extensive review, we refer readers to Peidro et al. (2009a) and Mula et al. (2010a).

Recently, [START_REF] Alemany | Mathematical programming model for centralized master planning in ceramic tile supply chains[END_REF] proposed a centralised mixed integer linear programming model for the master planning problem for the replenishment-production-distribution of SCs in the ceramic sector for the purpose of maximising the total gross margin in multi-supplier, multi-plant, multi-type, multi-level distribution centres with a multi-item and multi-period logistic environment. The contemplated model was inspired by an actual case and includes the most relevant characteristics of the ceramic sector. It also considers the several distribution levels integrated by nodes of various types. * * This research has been carried out in the framework of a project funded by the Science and Technology Ministry of the Spanish Government, entitled "Project of reinforcement of the competitiveness of the Spanish managerial fabric through the logistics as a strategic factor in a global environment". Ref. PSE-370000-2008-8. [START_REF] Alemany | Mathematical programming model for centralized master planning in ceramic tile supply chains[END_REF], this paper proposes a fuzzy multi-objective linear programming model (FMOLP) for the master planning problem for replenishment, production, and distribution in a ceramic tile SC. The objectives in conflict to simultaneously optimise are to: (i) maximise the fuzzy gross margin (incomes minus costs), (ii) minimise the fuzzy idle time of the production resources, and (iii) minimise the fuzzy backorder quantities. The aspiration levels of these three objectives can be considered fuzzy in nature because of the incompleteness and/or unavailability of the data required over the planning horizon, which could be obtained by being subjectively based on the planner's experience. Then, an interactive solution methodology to solve the FMOLP problem for the purpose of finding a preferred compromise solution is applied, which is based on the work of Torabi and Hassini (2008). For other solution methodologies to solve fuzzy multi-objective approaches, readers are referred to Zimmerman (1978), [START_REF] Lai | Possibilistic linear programming for managing interest rate risk[END_REF], [START_REF] Li | Computing efficient solutions to fuzzy multiple objective linear programming problems[END_REF] and [START_REF] Selim | Collaborative production-distribution planning in supply chain: a fuzzy goal programming approach[END_REF].

In relation to the previous work by [START_REF] Alemany | Mathematical programming model for centralized master planning in ceramic tile supply chains[END_REF], this fuzzy model is able to obtain better performance in terms of minimising the idle time of the production resources without lowering gross margin levels or increasing backorder levels, which is desirable for the entire SC. Other multiobjective applications for SC planning may be found in [START_REF] Chen | Multi-objective optimization of multi-echelon supply chain networks with uncertain product demand and prices[END_REF], [START_REF] Chern | A heuristic algorithm for master planning that satisfies multiple objectives[END_REF] and [START_REF] Selim | Collaborative production-distribution planning in supply chain: a fuzzy goal programming approach[END_REF].

The rest of this article is arranged as follows. Section 2 formulates a FMOLP model for the master planning problem in a ceramic tile SC. Section 3 describes its solution methodology. Section 4 validates and evaluates the proposed model by using an example based on a real-world problem. Section 5 offers conclusions and further research.

Fuzzy multi-objective linear programming for ceramic supply chain planning

Assumptions and nomenclature

This section formulates the FMOLP model for master planning in a ceramic tile SC. This model considers the following assumptions which are justified in detail in [START_REF] Alemany | Mathematical programming model for centralized master planning in ceramic tile supply chains[END_REF]:

• The flow of parts, components, raw materials (RMs) and finished goods (FGs) that might circulate between the nodes is known beforehand. • There are several production plants supplied with various RMs provided by different suppliers with a limited supply capacity. • Each production plant has one or several parallel production lines, which process different FGs with a limited capacity. • There are FGs that are manufactured only in production plants; others may be partly subcontracted, while some may be totally subcontracted to external suppliers. • For commercial reasons and for saving setups, FGs are grouped into product families, which are assigned to the production lines. • The minimum run lengths for product families are specified.

• The item setups also exist among the products belonging to the same product family.

• Minimum lot sizes for those products manufactured on a specific line are considered.

• The distribution of FGs from production plants to end customers is carried out in various stages by different types of distribution centres, such as central warehouses, logistic centres and shops. • Neither manufactured nor subcontracted FGs can be stored in manufacturing plants.

They are sent to the central warehouses, which have a limited storage capacity. • The outgoing FGs from the central warehouses cover the demand of both end customers and supply logistics centres. • Logistics centres supply FGs to the shops previously assigned to them.

• Shops attend to end customers' demands. The nomenclature defines the indices, sets of indices, parameters and decision variables (Table 1).

Objective functions

There are three objectives in conflict to simultaneously optimise.

Objective 1 (z 1 )

To maximise the total gross margin over the time periods computed by subtracting costs from the revenues produced by sales in central warehouses and shops. Total costs include transport and purchase costs of RMs; production costs of FGs; setup costs of product families and FGs in the production plants; transport costs of FGs; storage costs; and subcontracting costs of FGs. To minimise the backlogged demand over the whole planning horizon, expressed in the product units (m 2 ) in both warehouses and shops.

∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ - - - - - - - - -       + ≅ t q Wq(q) w Iw(w) i iqwt iqw t a Qa(a) q Iq(q) i iaqt iaq t i Bi(i) b Ab(b) a ibat ib t a
∑ ∑ ∑ ∑ ∑ ∑ ∈ ∈ + ≅ t w Iw(w) i iwt t a Ia(a) i iat DIFTK DIFA z Min 2 (2) Objective 3 (z 3 )
To minimise idle time, expressed in time units (hours), since production lots may be variable (although they must be above a specified minimum) and preparation times are fixed; thus the model should minimise the resulting idle time in terms of the preparation times employed when launching those lots that better adjust to demand which, on the other hand, will generate less storage costs.

( ) ∑∑ ∑ ∑ ∑ ∈ ∈ ∈         - + + ≅ t p p Lp l lpt l Il i ilp ilpt ilpt ilp l Fl f flpt flp caf tfab MP ZI tsetupi ZF tsetupf z Min ) ( ) ( ) ( 3 * * * (3)
For each objective function, the decision maker (DM) has imprecise objectives. Symbol " ≅ " is the fuzzified version of "=" and refers to the fuzzification of the aspiration levels. Accordingly, Eqs.

(1), ( 2) and (3) are fuzzy, and it is necessary for the DM to simultaneously optimise these conflicting objectives in the imprecise aspiration levels framework.

Constraints

In this section, we briefly review the constraints originally proposed by [START_REF] Alemany | Mathematical programming model for centralized master planning in ceramic tile supply chains[END_REF], which are necessary for our proposal.

Constraint (4) is the inventory balance equation for RMs.

∑ ∑ ∑ ∈ ∈ ∈ - - + = ) ( ) ( ) ( 1 ) * ( c Ic i p Lp l ilpt ic c Rc r crpt cpt cpt MP v CTP INC INC c,p,t ∀ (4)
Constraint ( 5) establishes safety stocks for RMs.

cp cpt ssc INC ≥ c,p,t ∀ (5)
Constraint ( 6) defines the supply capacity available for RMs suppliers.

crt p crpt ca CTP ≤ ∑ Rc(c),t c,r ∈ ∀ (6)
Constraint ( 7) establishes the capacity available for production lines.

( ) lpt l Il i ilpt ilp ilpt ilp l Fl f flpt flp caf MP tfab ZI tsetupi ZF tsetupf ≤ + + ∑ ∑ ∈ ∈ ) ( ) ( * * * Lp(p),t p, l ∈ ∀ (7)
Constraint ( 8) is related to the product families to be produced in each line.

∑ ∈ = ) ( f If i ilpt flpt MP MPF Fl(l),t Lp(p),f l , p ∈ ∈ ∀ (8)
Constraint (9) establishes the minimum lot sizes for FGs production.

t ), l ( Il i ), p ( Lp l , p X * lmi MP ilpt ilp ilpt ∈ ∈ ∀ ≥ (9)
Constraints ( 10) and ( 11) allocate products and product families to each line. Parameters M1 and M2 are large enough integer numbers.

t l Il i p Lp l p X M MP ilpt ilpt ), ( ), ( , * 1 ∈ ∈ ∀ ≤ (10) t l Fl f p Lp l p Y M MPF flpt flpt ), ( ), ( , * 2 ∈ ∈ ∀ ≤ (11)
Constraints ( 12) to (15) guarantee the control of the setup of FGs and product families.

t l Il i p Lp l p X X ZI ilpt ilpt ilpt ), ( ), ( , 1 ∈ ∈ ∀ - ≥ - (12) t ), p ( Lp l , p 1 X ZI i ilpt i ilpt ∈ ∀ - ≥ ∑ ∑ (13) t l
∈ ∀ = - ∑ ∑ ∈ ∈ (17)
Constraints ( 18) to ( 21) are related to subcontracting decisions. These constraints also ensure that the amount of subcontracted FGs is transported to the warehouses.

S * sc min CSC ibt ib ) b ( Ab a ibat ≥ ∑ ∈ Bi(i),t PFSP,b i ∈ ∈ ∀ (18) S * sc min CSC ibt ib ) b ( Ab a ibat ≥ ∑ ∈ Bi(i),t PFST,b i ∈ ∈ ∀ (19) S * casc CSC ibt ibt ) b ( Ab a ibat ≤ ∑ ∈ Bi(i),t PFSP,b i ∈ ∈ ∀ (20) S * casc CSC ibt ibt ) b ( Ab a ibat ≤ ∑ ∈ Bi(i),
∑ ∈ ) ( a,t ∀ (23)
Constraints ( 24) and ( 25) are inventory balance equations for the FGs in warehouses.

∑ ∑

∈ ∈ - - - + = ) ( ) ( 1 a Qa q iaqt iat a Pa p ipat iat iat CTCL VEA CTA INA INA PFNS, a, t i ∈ ∀ ( 24
) CTCL VEA CSC CTA INA INA ) a ( Qa q iaqt iat ) i ( Bi b ) a ( Ba b ibat ) a ( Pa p ipat 1 iat iat ∑ ∑ ∑ ∈ ∈ ∧ ∈ ∈ - - - + + = PFSP, a, t i ∈ ∀ (25)
Constraint ( 26) is similar to ( 24) and ( 25), but also ensures that subcontracted FGs only come from FG suppliers.

9 CTCL VEA CSC INA INA ) a ( Qa q iaqt iat ) i ( Bi b ) a ( Ba b ibat 1 iat iat ∑ ∑ ∈ ∈ ∧ ∈ - - - + = PFST, a, t i ∈ ∀ (26)
The backorder quantities in warehouses are calculated using Constraint (27).

iat iat iat iat da DIFA DIFA VEA = - + -1 Ia(a),t a, i ∈ ∀ (27)
Constraint ( 28) limits the backorder quantities in warehouses.

iat iat da DIFA * 1 α ≤ Ia(a),t a, i ∈ ∀ (28)
Constraints ( 29) and ( 30) are the inflows and outflows of FGs through each logistic centre and shop, respectively. 31) determines the backorder quantities in shops.

∑ ∑ ∈ ∈ = ) ( ) ( q Wq w iqwt q Aq a iaqt CTTK CTCL Iq(q), t q,i ∈ ∀ (29) iwt iqwt VETK CTTK = Iw(w),t Qw(w),i w,q ∈ ∈ ∀ (30) Constraint (
Iw(w),t w,i dtk DIFTK DIFTK VETK iwt 1 iwt iwt iwt ∈ ∀ = - + - (31)
Constraint ( 32) limits the backorder quantities in shops.

Iw(w),t w,i dtk * 2 DIFTK iwt iwt ∈ ∀ ≤ α (32)
The model also contemplates non-negativity constraints and the definition of binary variables (31).

MPF flpt , MP

) T t B, b R, r W, w Q, q A, a P, p L, l C, c I, i F, f ∈ ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈ ∀
Finally, some decision variables can be defined as integers, but they could change depending on the real-world problem to which the model is applied.

Solution methodology

In this section, the solution methodology proposed by Peidro et al. ( 2010) is adopted as the basis of this work. This previous work defines an approach to transform a fuzzy multi-objective linear programming model (FMOLP) into an equivalent auxiliary crisp mathematical programming model for a supply chain transport planning problem in the automobile industry. This approach adopts linear membership functions to represent all the fuzzy objective functions together with the fuzzy programming solution method of Torabi and Hassini (TH) (2008).

According to Torabi and Hassini [3], a multi-objective model could be transformed into a singleobjective model as follows:

Max ) ( ) ( ) ( x x k z k k ∑ - + = µ θ γ γλ λ 1 0 s.t. n k x k z ,..., ) ( 1 0 = ≤ µ λ F o r P e e r R e v i e w O n l y 10 ) (x F x ∈ ] , [ , 1 0 0 ∈ γ λ ( 34
)
where

. , 0 1 > = ∑ k k k θ θ
The compensation coefficient determines whether the solution obtained will be balanced; in other words, if it gives the same importance to all the objectives; or if it is unbalanced; that is, prioritising those objectives whose θ k. weight is greater. Any parameter γ values close to one will provide more balanced solutions; otherwise, values close to zero will generate solutions with better results for those objectives whose weight is greater. It is the DM's task to determine which type of solutions he/she wishes to obtain when applying the present solution methodology.

The interactive solution procedure, adopted from [START_REF] Peidro | An interactive fuzzy multi-objective approach for operational transport planning in an automobile supply chain[END_REF], is summarised as follows:

Step 1. Formulate the original FMOLP model according to Eqs. ( 1) to (33).

Step 2. Specify the corresponding linear membership functions for all the fuzzy objectives (upper and lower limits).

Step 3. Determine the corresponding relative importance of the objective functions (θ k ) and the compensation coefficient (γ).

Step 4. Transform the original FMOLP problem into an equivalent single-objective MILP form using the TH fuzzy programming method.

Step 5. Solve the proposed auxiliary crisp single-objective model by the MIP solver and obtain the initial solution.

Step 6. If the DM is satisfied with this current efficient compromise solution, stop. Otherwise, go back to Step 2 and provide another efficient solution by changing the controllable parameters value (γ, and θ k ).

Application to a ceramic tile supply chain

This section uses the example provided by [START_REF] Alemany | Mathematical programming model for centralized master planning in ceramic tile supply chains[END_REF] to validate and evaluate the results of our proposal. It is an SC that is representative of the ceramic tile sector. There are 3 production plants which produce 4 FGs grouped into 3 product families whose rates, minimum run lengths and fixed costs are provided. Each plant has 2 production lines. All the product families may be manufactured on the production lines in the 3 plants. Six-week periods are considered in the planning horizon. In addition, the following information is provided: bill of materials, transportation costs, setup costs, initial inventory, available production and storage capacities, raw material costs, safety stocks, inventory costs, setup times, production costs, sale prices, subcontracting costs, backorder costs, production run times, minimum lot sizes and demand. Details of this data used can be found in [START_REF] Alemany | Mathematical programming model for centralized master planning in ceramic tile supply chains[END_REF].

Implementation and resolution

The proposed model has been developed in the MPL language, V4.2. The resolution has been carried out with optimisation solver CPLEX 9. Finally, the input data and the model solution values have been processed with the Microsoft Access database (2007). The experiment was run on a PC with a 2.40 GHz processor and 4 GB of RAM, and a limited calculation time of 120 seconds was defined. In the above-described solution methodology application, the DM provides the relative importance of the objectives as: θ 1 > θ 2 > θ 3 , θ = (0.5, 0.35, 0.15). In this case, the total gross margin is more important for the DM over the time periods followed by minimising backlogged demand and, finally, by minimising idle time. Thus, an unbalanced compromise solution with the highest degree of satisfaction for z1 and z2 is of particular interest.

Evaluation of the results

This section compares the results obtained by the mixed integer linear programming model proposed in [START_REF] Alemany | Mathematical programming model for centralized master planning in ceramic tile supply chains[END_REF] with those obtained by the heuristic procedure and the FMOLP solution methodology proposed in this work.

[ u l z z 1 1 , ] l z 1 = 200,000 u z 1 = 370,000 [ u l z z 2 2 , ] l z 2 = 0 u z 2 = 500 [ u l z z 3 3 , ] Not applicable l z 3 = 0 u z 3 = 620
Table 2 presents the comparisons of the results obtained by the two approaches for all the objectives identified. Furthermore, details are provided of the total solution time employed by each approach, the minimum degree of satisfaction of the objectives (λ 0 ), the equivalent crisp model goal value according to the TH approach (λ(x)), the degree of satisfaction of the goal functions, as well as the upper and lower limits specified by the DM to define the linear membership functions that are representative of the fuzzy objectives. The λ 0 , λ(x) µ z1 , µ z2 , µ z3 , values and the upper and lower limits are only significant for the proposed approach, and are not applicable to the deterministic model.

As Table 2 shows, the results obtained by the approach proposed for the case of γ=0.1 are the most appropriate. As mentioned previously, it is worth remembering that when the DM defines the weights of the objectives, he/ she should have opted for a more unbalanced solution that prioritises those objectives with greater weights. In the case of γ=0.1, the two proposed approaches offer 12 similar results for objectives z1 and z2, and substantial improvement (40.25%) in idle time terms. In this particular case, the minimum degree of satisfaction is lower (λ 0 = 0.4052) than for the other compensation coefficient values. This is mainly due to the proposed model prioritising objectives z1 and z2 more than z3 (whose weight is θ 3 = 0.15); z3 sets the minimum degree of satisfaction value.

For those compensation coefficient values above 0.1 (γ > 0.1), the proposed approach performs in a similar fashion. For these γ values, the proposed model tends to obtain more balanced solutions, without paying attention to the weights defined by the DM. This is the reason why the result for objective z3 considerably improves (an improvement of up to 70% if compared to the deterministic model), and obtains higher minimum degree of satisfaction values. In contrast, the objective z1 values slightly worsen if compared to the deterministic model (by up to 5.2%). The λ(x) value lowers as the compensation coefficient value increases. According to Eq. ( 11), this is due to the fact that when γ increases, the weight of the minimum degree of satisfaction becomes higher than the result of λ(x), and is below the weight of the degrees of satisfactions of all the objectives as a whole.

The approach proposed in this work performs than the single-objective deterministic approach proposed in [START_REF] Alemany | Mathematical programming model for centralized master planning in ceramic tile supply chains[END_REF]. The defined solution methodology offers better results than the deterministic model, similar results in terms of maximising the total gross margin, equal results as regards minimising backlogged demand (which is zero in all cases), and better results when minisiming idle time. When it comes to solving the model, it also offers flexibility for the DM as it enables both unbalanced and balanced solutions. It is the DM who must decide which combination of parameters best satisfies his/her needs.

Conclusions

This paper continues the work of [START_REF] Alemany | Mathematical programming model for centralized master planning in ceramic tile supply chains[END_REF] and proposes an FMOLP model to solve the master planning problem for the replenishment-production-distribution of SCs in the ceramic sector for the purpose of maximising the total gross margin, minimising backorder quantities and minimising the idle production time in multi-supplier, multi-plant, multi-type, multi-level distribution centres with a multi-item and multi-period logistic environment. For the purpose of solving the FMOLP model, we followed an interactive solution methodology.

The interactive solution methodology yields an efficient compromise solution and presents the DM's overall satisfaction with the determined goal values in a multi-objective master planning problem for a centralised replenishment, production and distribution ceramic tile supply chain. This approach provides solutions that are consistent with the DM's preferences (i.e., the consistency between weight vector θ k and the satisfaction vector) because it is able to find different efficient solutions for a specific problem with a given weight vector θ k by changing the γ value. We show that a compromise solution is obtained which, in relation to the previous work [START_REF] Alemany | Mathematical programming model for centralized master planning in ceramic tile supply chains[END_REF], improves the utilisation of the production resources by minimising the idle time without considerably lowering the gross margin and/or increasing backorder quantities. Moreover, our fuzzy model solution is robust as regards the fuzziness in the problem without causing an explosive growth of computational effort. At this point, the use of metaheuristics to improve the efficiency of the solution methodology will be a forthcoming work. Finally, this work assumes that the linear membership functions for the related objective functions are reasonably provided. In real-world situations, however, the DM should generate suitable membership functions based on subjective judgement and/or historical resources. Future studies on non-linear membership functions to solve the multiobjective master planning problem in a ceramic tile SC may apply. Manufacturing Technology, 36, 401-418. Torabi, S., Hassini, E., 2008, An interactive possibilistic programming approach for multiple objective supply chain master planning. Fuzzy Sets and Systems, 159, 193-214. von Lanzenauer, Ch., Pilz-Glombik, H.K., 2002, Coordinating supply

•

  Limited backorders are permitted in both central warehouses and shops.

  of satisfaction of the kth objective function and the minimum degree of satisfaction of the objectives, respectively. Moreover, θ k and γ indicate the relative importance of the kth objective function and the compensation coefficient, respectively. The θ k parameters are determined by the DM based on her/his preferences so that

  ., Ozkarahan, I., 2008, A supply chain distribution network design model: An interactive fuzzy goal programming-based solution approach. The International Journal of Advanced

  chain decisions: an optimization model. OR Spectrum 24 (1), 59-78. Zimmermann, H., 1978, Fuzzy programming and linear programming with several objective functions. Fuzzy Sets and Systems, 1, 45-46. Page 14 of 14 http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.uk International Journal of Production Research

Table 1 .

 1 Nomenclature.

	Indices
	c	RMs, items, and components (c=1…C)	q	Logistic centres (q=1…Q)
	i	FGs (i=1…I)	w	Shops (w=1…W)
	f	Product families (f=1…F)	r	Suppliers of RMs, items, and components
	l	Production lines (l=1…L)		(r=1…R)
	p	Production plants (p=1…P)	b	Suppliers of FGs (b=1…B)
	a	Warehouses (a=1…A)	t	Time periods (t=1…T)
	Sets of Indices
	Il(l)	Set of FGs that can be manufactured on	Lp(p)	Set of manufacturing lines that belong to
		manufacturing line l		production plant p
	Fl(l)	Set of product families that can be manufactured	Pa(a)	Set of production plants that can send FGs to
		on manufacturing line l		warehouse a
	If(f)	Set of FGs that belong to product family f	Aq(q)	Set of warehouses that can supply logistics centre
				q
	Ip(p)	Set of FGs that can be produced in production	Rc(c)	Set of suppliers that can supply RM c
		plant p
	Ia(a)	Set of FGs that can be stored in warehouse a	Rp(p)	Set of suppliers of RMs that can supply
				production plant p
	Ic(c)	Set of FGs of that RM c form part	Cr(r)	Set of RMs that can be supplied by supplier r
	PFNS Set of FGs that cannot be subcontracted	Qa(a)	Set of logistics centres that can be supplied by
				warehouse a
	PFSP Set of FGs that can be subcontracted either	Wq(q) Set of shops that can be supplied by logistics
		partially or completely		centres q
	PFST Set of FGs that are compulsorily subcontracted	Qw(w) Set of logistics centres capable of supplying shop
		completely		w
	Iq(q)	Set of FGs that can be sent to logistics centre q	Bi(i)	Set of suppliers of FGs i to which the FG may be
				subcontracted
	Iw(w) Set of FGs that can be sent to shop w	Ba(a)	Set of suppliers of FGs that can supply warehouse
				a
	Lf(f)	Set of manufacturing lines that may produce	Ab(b)	Set of warehouses that can be supplied by the
		product family f		supplier b of FGs
	Model Parameters
	ca crt	Capacity (units) of supplying RM c of	M1,M2		Very large integers
		supplier r in period t
	costtp crp	Cost of purchase and transport of one unit	capal a		Storage capacity (m 2 ) in warehouse a
		of RM c from supplier r to production plant
		p
	caf lpt	Production capacity available (time) of	costtcl iaq	Cost of transporting one m 2 of FG i from
		production line l at plant p during time			warehouse a to logistics centre q
		period t
	cm i	Loss ratio of FG i. It represents the percentage of faulty m 2 obtained due to the	costina ia	Cost of making an inventory of one m 2 of FG i in the warehouse over a time period
		intrinsic characteristics of the production
		process in the ceramics sector.
	cq i	First quality coefficient of FG i. It represents the percentage of m 2 that can be	costdifa ia	Cost of backordering one m 2 of demand of FG i in warehouse a over a time period
		sold as first quality.
	costta ipa	Cost of transporting one m 2 of FG i from	pa ia		Sales value of one m 2 of FG i in warehouse
		production plant p to warehouse a			a
	costp ilp	Cost of producing one m

2 of FG i on production line l of production plant p da iat External demand (m 2 ) of FG i at the warehouse a in period t costsetupf flp Setup costs for product family f on production line l of production plant p ssa ia Safety stock (m 2 ) of FG i at warehouse a cp Safety stock of RM c in production plant p 2 α Maximum backorder quantity permitted in a period in shops expressed as a percentage of the demand of that period casc ibt Supply capacity (m 2 ) of FG i of supplier b over time period t Decision Variables CTP crpt Amount of RM c to be purchased and transported from supplier r to production plant p over time period t INA iat Inventory (m 2 ) of FG i in warehouse a over time period t INC cpt Inventory of the RM c at plant p at the end of time period t CSC ibat Amount (m 2 ) of FG i subcontracted to supplier b for warehouse a over time period t MPF flpt Amount (m 2 ) of product family f manufactured on production line l of production plant p over time period t S ibt Binary variable with a value of 1 if FG i is subcontracted to supplier b over time period t

  ilpt , CTP crpt , CTA ipat , INA iat , INC cpt , CTCL iaqt , CTTK iqwt , VEA iat , DIFA iat , VETK iwt , DIFTK iwt , CSC ibat ≥ 0 and,

	X	ilpt	,	Y	flpt	,	ZF	flpt	,	ZI	ilpt	,	S	ibt	∈	{ } 1 , 0	(33

Table 2 .

 2 Comparison of the results.

	Item	Deterministic	Proposed	Proposed	Proposed	Proposed	Proposed
			method	method	method	method	method
			(γ=0.1)	(γ=0.3)	(γ=0.5)	(γ=0.7)	(γ=0.9)
		z1=368,232.4 €	z1=364,246.6 €	z1=349,243.1 €	z1=349,172.9 €	z1=349,117.2 €	z1=349,079.2 €
	Objective	z2=0	z2=0	z2=0	z2=0	z2=0	z2=0
	values	z3=617.23 h.	z3=368.78 h.	z3=185.30 h.	z3=186.30 h.	z3=184.96 h.	z3=185.30 h.
	Total
	solution	40.76 sec.	120 sec.	120 sec.	120 sec.	120 sec.	120 sec.
	time
	λ0 value		λ0 = 0.4052	λ0 = 0.7011	λ0 = 0.6995	λ0 = 0.7016	λ0 = 0.7011
	λ(x)value		λ(x) = 0.8233	λ(x) = 0.7752	λ(x) = 0.7528	λ(x) = 0.7330	λ(x) = 0.7116
	µz1		0.9662	0.8780	0.8775	0.8772	0.8769
	µz1		1	1	1	1	1
	µz3		0.4052	0.7011	0.6995	0.7016	0.7011

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.ukInternational Journal of Production Research

This research has been carried out in the framework of a project funded by the Science and Technology Ministry of the Spanish Government, entitled "Project of reinforcement of the competitiveness of the Spanish managerial fabric through the logistics as a strategic factor in a global environment". Ref. PSE-370000-2008-8.