Andrzej T Kowal
email: andrzej.t.kowal@pwr.wroc.pl

geometry and anharmonic vibrational spectra of PbX 2 ••NO and PbX 2 ••ON (X=F, Cl, Br, I) complexes

Keywords: lead(II) halide, anharmonic, effective core potential, VSCF, nitrogen monoxide lead(II) halide, effective core potential, anharmonic, nitrogen monoxide

    Computational study of the equilibrium geometry and anharmonic vibrational spectra of PbX2**NO and PbX2**ON (X=F, Cl, Br, I) complexes

INTRODUCTION

Blue shift of the intra-ligand stretching mode of σ-donors like CO, NO and N 2, observed in the infrared spectra of fourth row divalent metal halides MX 2 (M=Ca, Cr, Mn, Ni, Cu, Zn) upon MX 2 ••L complex formation [1,2] prompted increased interest in such complexes. Subsequently, similar effect has been detected in the low temperature Ar matrix infrared spectra of the species formed by SnX 2 and PbX 2 halides (X=F, Cl, Br, I) with carbon monoxide, dinitrogen and nitrogen monoxide [3]. Analysis of their infrared spectra led to the conclusion that the structure of SnX 2 ••L and PbX 2 ••L complexes consists of a bent MX 2 moiety with end-on coordinated ligand molecule and the blue shift of intra-ligand stretching mode can be explained in terms of electron density transfer from anti-bonding σ * orbitals of the ligand to empty d orbitals of the metal [3]. Ab initio study of PbX 2 ••L compounds (L=CO, N 2 ; X=F, Cl, Br, I) at MP2 level with relativistic effective core potentials (RECP) and basis sets of SBKJC+(d) type afforded description of equilibrium geometry of these species and correctly reproduced, at anharmonic level, shift of ν(N≡N) and ν(C=O) stretching modes to higher energy upon complex formation [4]. Electronic structure of nitrogen monoxide [5] and its complexes with transition metals [6,7] have been extensively studied both experimentally and by computational chemistry methods. Interaction of NO with silver [8], gold [9], and platinum [10] small metal clusters as well as with transition metal clusters embedded in zeolites [11] has also been a subject of considerable interest with relation to heterogeneous catalysis. Variety of binding modes of the NO molecule to metal atoms, resulting in linear Ncoordinated, linear O-coordinated, bent N-or O-coordinated, and side-on η 2 coordinated structures have been either observed experimentally [6][7][8][9][10][11] or investigated by quantum chemistry methods [6,12]. Amongst the computational chemistry methods employed to the elucidation of electronic structure and properties of open-shell systems like nitrogen monoxide complexes with metals and metal ions, perturbational treatment of electron been preferred to unrestricted Möller-Plesset (MP) perturbation theory because of better convergence properties and lesser sensitivity to spin contamination [13][14][15][16][17]. Characteristics of Z-averaged perturbation theory (ZAPT) [13,14], based on ROHF reference wave function, appear particularly interesting as the computational demand of the method is comparable to that of the closed-shell MP method [14] and it is entirely devoid of spin contamination at second order (ZAPT2) [14]. Moreover, since the analytic gradient expressions for second order ZAPT theory have been derived [15] and implemented in revised form within distributed data interface [16], thus enabling parallel execution of the code on computer clusters, the method seems even more attractive for treatment of large open-shell molecular systems and detailed PES scans on such molecules. Binding mode of the NO molecule and geometry of the complexes have been commonly inferred from the infrared spectra of species trapped in low temperature rare gas matrices [3,6] aided with density functional theory (DFT) calculations of the harmonic approximation vibrational spectra [6]. Properties of metal cluster bound NO have been studied by DFT method [8,10,12] or first-principles based simulation packages like VASP [9,11]. In most cases, vibrational spectra computed in the harmonic approximation have been used to establish NO binding mode and provide support to experimental transition assignments, except for the most recent VASP study on nitric oxide adsorption on Pd clusters in mordenite [11], where one-dimensional treatment of NO anharmonicity has been applied. The importance of anharmonic treatment of nitric oxide complexes becomes evident from the analysis of the latter data [11], where the red shift of ν(N=O) stretching mode of Pd cluster bound NO is contrasted with the blue shift of this mode in mordenite supported Pd clusters.

Vibrational self-consistent field (VSCF) theory [18][19][20][21][22], an analog of self-consistent field method in electronic structure theory, provides a route to computation of anharmonic approximation vibrational spectrum by taking into account intrinsic anharmonicity related to vibrational motion along single normal mode in the mean field of the remaining modes, as well as that originating in inter-mode coupling. PES points required by the direct VSCF method are computed on a rectangular grid derived from normal mode displacements in rectilinear [23,24] or curvilinear (internal) [25] coordinates along single modes and mode pairs (and possibly mode triples) by electronic structure method of choice. It has been shown that VSCF method based on PES grid derived from normal mode displacements in curvilinear coordinates provides better description of low energy large amplitude vibrations [25] as compared to that employing grid deduced from rectilinear coordinates [23,24]. Correlation corrected extensions of VSCF method, termed CC-VSCF, based on either virtual configuration interaction (VCI) including singly and doubly excited states [23,26] or secondorder non-degenerate [23,26] and degenerate [27,28] perturbation theory account for intermode correlation effects. CC-VSCF method has been successfully applied to weakly bound species like rare gas (Rg) hydrogen fluoride HRgF complexes [29], HXeI [30], HXeOH [31], as well as to PbX 2 ••L complexes (X= F, Cl, Br, I; L= CO, N 2 ) [4]. Noticeably, VSCF method provided correct explanation of the anomalous isotope effect in the spectrum of HXeOH [31], and a blue shift of ν(N≡N) and ν(C=O) stretching modes upon PbX bonding isomers (X=F, Cl, Br, I). Anharmonic wavenumber shifts of ν(N=O), ν s (Pb-X), and ν as (Pb-X) fundamentals upon PbX 2 ••NO complex formation are compared to those observed in the low temperature Ar matrix spectra [3] and correspondingly to those of PbX

METHODS

All computations of equilibrium geometry, harmonic approximation vibrational spectra, CC-VSCF fundamentals and overtone transitions were done within electronic structure suite of programs GAMESS [32]. Equilibrium geometry parameters of PbX 2 ••NO and PbX 2 ••ON (X=F, Cl, Br, I) complexes were evaluated under C s symmetry constraints using parallel implementation [16] of second-order Z-averaged perturbation theory (ZAPT2) [13][14][15] and Stevens-Basch-Krauss-Jasien-Cundari (SBKJC) scalar-relativistic effective core potentials (RECP) and valence basis sets on all atoms [33][34][35][36]. Replacement of all-electron basis sets by effective core potentials does not lead to significant loss of accuracy even in case of 2-nd and 3-rd row elements [37]. SBKJC basis set was augmented with single d-type polarization and single L-type diffuse functions, denoted as SBKJC+(d). Harmonic approximation vibrational spectra of the complexes were computed from Hessians evaluated by numerical differentiation of analytic ZAPT2 gradients. Two-mode representation of the VSCF system potential [23,26,27], including single mode (diagonal) V i diag (Q i ) and mode-mode coupling V ij coup (Q i , Q j ) terms was assumed, based on adequate description of PbX 2 ••L anharmonicity [4] by such truncated expression:

∑ ∑∑ - N j N i N i > j j i coup ij j diag j N ) Q , (Q V + ) (Q V = ) Q , V(Q F o r P e e r R e v i e w O n l y 0.5 ]
range (where ω i is the frequency of i-th normal mode). Direct CC-VSCF method required calculation of 9360 PES points for each species. Intensities of fundamental transitions were evaluated using ZAPT2 dipole moment surfaces and VSCF wavefunctions of the corresponding states,

2 0 3 3 8π | ) (Q Ψ | ) (Q µ | ) (Q Ψ | ω hc N = I i (m) i i i ) ( i i a i
r where ψ i (0) and ψ i (m) denote VSCF wavefunctions of the ground and the m-th excited vibrational states of i-th normal mode, ω i is the CC-VSCF vibrational frequency of this mode, and the remaining symbols have their regular meaning [26]. (267.7 pm) [46] ions, computed at MP2 level in scalar relativistic ECP basis [45,46] Table 2) compares favorably against that computed by the same method in triple-zeta valence basis set, TZ2P (1895.8 cm -1 ) [14]. VSCF computed fundamental of this mode (1869.1 cm -1 ), including intrinsic (diagonal) anharmonicity, appears reasonably close to the transition wavenumber observed in the Ar matrix spectrum of NO (1875.0 cm -1 ) [3], showing deviation typical of VSCF method for stretching modes [23][24][25][26][27]. Since there is no contribution to anharmonic correction from mean field effect or inter-mode correlation in the case of molecule with single vibrational mode, the aforementioned wavenumber of ν(N=O)

RESULTS AND DISCUSSION

fundamental will be used to evaluate anharmonic shift of ν(N=O) mode upon complex formation. VSCF computed anharmonic spectra of PbF 2 ••NO and PbF 2 ••ON complexes (Table 2) provide fundamentals corrected for intrinsic (diagonal) anharmonicity, mean field and mode-mode coupling (VSCF), and mode correlation at degenerate second order perturbational (DPT2-VSCF) level. For the reason stated above, ν(N=O) fundamental of the coordinated NO, corrected for intrinsic anharmonicity, will be used to assess wavenumber shift due to complex formation. Fundamentals of PbF 2 , computed at MP2/SBKJC+(d) level (Table 2), show negligible deviation from low temperature Ne matrix spectrum [42] whereas those of ν s (Pb-F) and ν a (Pb-F) modes in PbF 2 ••NO complex overestimate transition wavenumbers measured in low temperature Ar matrix [3] by ~10 cm -1 , partly due to matrix polarizability effects [START_REF] Hargittai | [END_REF]. It seems interesting to note that the harmonic approximation provides qualitatively incorrect description of ν(N=O) mode shift on PbF 2 ••NO complex formation (-11.7 cm -1 ) whereas the wavenumber difference of the respective fundamentals (+15.8 cm -1 ) nears that deduced from experiment (+16.4 cm -1 ) [3]. Although wavenumber shifts of ν s (Pb-F) and ν a (Pb-F) fundamentals on going from PbF 2 to PbF 2 ••NO are lower than those estimated in the harmonic approximation, they are still larger than shifts assessed experimentally [3]. Wavenumber of high intensity ν(Pb-N) fundamental at 232.8 cm -1 is 10 cm -1 lower than its harmonic counterpart, with majority of anharmonic correction coming from intrinsic anharmonicity and mean field / mode coupling effects. Same pattern of anharmonic corrections can be observed for all PbF represented as a two-dimensional mapping of V i,j coup (Q i , Q j ) mode coupling potential (Figure 3) on NxN PES grid,

              ∑∑ ∑ N = p N = q jq ip N = p ip j i otherwise | ) n , V(n | j = i if | ) V(n | N = ) Q , P(Q 1 1 1
where n ip and n jq are grid points and N is the number of grid points (N=16), reveals that intermode coupling contribution to VSCF anharmonic correction seems substantial for ν s (Pb-F) and ν a (Pb-F), ν(Pb-N) and ρ w (PbF 2 ), ρ w (PbF 2 ) and ρ t (PbF 2 ) mode pairs. Since no experimental evidence of PbX 2 ••ON complexes has been reported to date, their anharmonic spectra can be compared against those of PbX 2 ••NO isomer or structurally similar lead(II) compounds. Fundamental transitions of PbF 2 ••ON isomer (Table 2) can be viewed as characteristics of O-coordinated NO molecule, as evidenced by 33.4 cm -1 red shift of ν(N=O) stretching mode upon complex formation and less pronounced shifts of ν s (Pb-F) and ν a (Pb-F) fundamentals on going from PbF 2 to PbF 2 ••ON. Wavenumber of ν(Pb-O) fundamental, computed at 183.2 cm -1 , shows significant anharmonic correction due to intrinsic anharmonicity and mean field / mode coupling and appears remarkably close to the lowest energy ν(Pb-O) stretching mode of Pb 6 O(OH) 6

4+ ion (179 cm -1 ) [46], evaluated in the harmonic approximation at MP2/ECP level. Anharmonic corrections of PbF 4). Therefore, a reversal of experimental assignment placing ν s (Pb-Br) mode higher in energy than ν a (Pb-Br) would be consistent with VSCF calculation and VSCF computed intensity pattern of these two fundamentals (Table 4). Fundamentals of PbBr 2

Anharmonic spectrum of PbI 2 ••ON complex

Fundamental transitions of PbI 2 computed by DPT2-VSCF method (Table 5) are weakly anharmonic and their wavenumbers do not differ appreciably from Ar and Kr matrix spectra extrapolated to zero matrix polarizability [START_REF] Hargittai | [END_REF]. Since the VSCF procedure did not converge in the case of PbI 2) and the relative strength of the potential is coded in shades of gray (0, white and 255, black).

Figure 4. Two-dimensional mapping of V i,j coup (Q i , Q j ) mode coupling potential of PbF [START_REF] Hargittai | [END_REF] b -r e Ref. [START_REF] Hargittai | [END_REF] c -Ref. [41] d -Ref. [40] e -r e (NO)=115.1 pm Ref. [ a -Ar matrix data Ref. [3] b -Ne matrix data Ref. [42] cνstretching, δ -deformation, ρ r -rocking, ρ wwagging, ρ t -twisting, τ -torsional, ipin-phase mode, op -out-of-phase mode a -Ar and Kr matrix data extrapolated to zero polarizability Ref. [START_REF] Hargittai | [END_REF] b -Ar matrix data Ref. [3] cν -stretching, δ -deformation, ρ r -rocking, ρ wwagging, ρ t -twisting, τ -torsional, ipin-phase mode, op -out-of-phase mode 2) and the relative strength of the potential is coded in shades of gray (0, white and 255, black).

  restricted open-shell Hartree-Fock (ROHF) reference wave function has

3 . 2 .

 32 , are very close to those of PbX 2 ••ON series. One can also note that lengthening of r(Pb-L) distance and increase in (Pb-L-L) angle in both PbX 2 ••NO and PbX 2 ••ON isomers follows the order of decreasing electronegativity of halogen atoms and, in consequence, decreasing positive charge on Pb atom. Anharmonic spectra of PbF 2 ••NO and PbF 2 ••ON complexes Harmonic wavenumber of ν(N=O) mode computed at ZAPT2/SBKJC+(d) level (1890.1 cm -1 ;

  3 cm -1 , respectively, while formation of PbBr 2 ••ON complex shifts both modes to lower energy by -3.1 cm -1 . A decrease in the energy of ν(Pb-N) and ν(Pb-O) fundamentals in PbBr 2 ••NO and PbBr 2 ••ON isomers with relation to those of the fluoride and chloride analogs follows the order of decreasing positive charge on lead(II) ion and, in consequence, decreasing Pb-N and Pb-O bond strengths. Interestingly, sensitivity of the vibrational modes involving PbBr 2 moiety as a whole, ρ w (PbBr 2 ) and ρ t (PbBr 2 ), to the NO coordination mode can be viewed as a demonstration of difference in Pb-N and Pb-O bond strengths.

Figure 1 .

 1 Figure captions

Figure 2 .

 2 Figure 2. Equilibrium geometry of PbX 2 ••ON complexes computed at ZAPT2/SBKJC+(d)

Figure 3 .

 3 Figure 3. Two-dimensional mapping of V i,j coup (Q i , Q j ) mode coupling potential of PbF 2 ••NO

Figure 1 .Figure 2 .

 12 Figure 1. Equilibrium geometry of PbX 2 ••NO complexes computed at ZAPT2/SBKJC+(d) level (Molden [44] drawing).

Figure 3 .

 3 Figure 3. Two-dimensional mapping of V i,j coup (Q i , Q j ) mode coupling potential of PbF 2 ••NO at ZAPT2/SBKJC+(d) level. Modes are labeled in the order of decreasing harmonic wave number (first column of Table2) and the relative strength of the potential is coded in shades of gray (0, white and 255, black).

Figure 4 .

 4 Figure 4. Two-dimensional mapping of V i,j coup (Q i , Q j ) mode coupling potential of PbF 2 ••ON at ZAPT2/SBKJC+(d) level. Modes are labeled in the order of decreasing harmonic wave number (first column of Table2) and the relative strength of the potential is coded in shades of gray (0, white and 255, black).

Page 29 of 29 URL

 29 : http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics

  2

  ••L complex formation [4]. Having this in mind, one can anticipate that anharmonic level description of open-shell PbX 2 ••NO and PbX 2

  ••ON species (X=F, Cl, Br, I) vibrational spectra by correlation corrected VSCF method seems feasible.

	Present work reports on ZAPT2/SBKJC+(d) level computation of equilibrium geometry
	parameters and CC-VSCF anharmonic vibrational spectra of PbX 2 ••NO and PbX 2 ••ON

  2

  the difference is almost negligible, even when r e < r g relationship is taken into account[START_REF] Hargittai | [END_REF]. A noticeable shortening of r(N=O) distance upon PbX 2 ••NO complex formation (0.2 -0.3 pm) occurs in concord with halide dependent lengthening of r(Pb-X) distance, amounting to 1.4pm (PbF 2 ••NO) -2.3 pm (PbI 2 ••NO) (Table1). Opposite change or no change in r(N=O) bond length is observed on formation of PbX 2 ••ON species in accord with less pronounced increase in r(Pb-X) distance (0.7 -1.0 pm) on going from fluoride to iodide. Equilibrium geometry of PbX 2 ••NO species (Figure1), characterized by r(Pb-N) bond distance within 266.6 -271.7 pm, and (Pb-N-O) angle of 109.2 -120.7 deg, differs considerably from that of PbX 2 ••N 2[4]

	3.1. Equilibrium geometry It seems worth noting that equilibrium bond length of free NO computed at ZAPT2/SBKJC+(d) level (Table 1) is longer than experimental bond distance of 115.1 pm [38] by 4.1 pm, and longer by 3.4 pm with relation to that calculated by ZAPT2 method using triple-zeta quality basis set TZ2P [14]. Although overestimation of r(N=O) distance appears significant, it seems of secondary importance in computation aimed at anharmonic level o r P e e r R e v i e w O n l description of vibrational transitions in PbX F y

2 ••NO and PbX 2 ••ON complexes, including blue shift of ν(N=O) stretching mode upon complex formation. On the other hand, equilibrium bond distances and angles of lead halides evaluated under C 2v symmetry by MP2 method in the same basis set (Table

1

) overestimate experimental gas phase geometry parameters derived from electron diffraction

[START_REF] Hargittai | [END_REF]

to a less significant extent with the exception of PbF 2, where r(Pb-F) distance is longer than experimental r g value by 1.4 pm (Table

1

). Nonetheless, complex by much shorter Pb-N bond length and less obtuse Pb-L-L angle (where L-L represents ligand molecule). Moreover, computed (Pb-N-O) angle in PbX 2 ••NO series nears that characteristic of N-coordinated NO molecule with sp 2 hybridized nitrogen atom in transition metal complexes

[12]

. Isomer complex molecule, PbX 2 ••ON (Figure

2

) shows slightly longer Pb-O bond length (Table

1

), varying between 267.8 and 275.8 pm, and more obtuse (Pb-O-N) angle, as compared to PbX 2

••NO species. Despite remote similarity, it can be noted that the longest Pb-O bond distances of Pb 4 (OH) 4 4+ (252 pm) [45] and Pb 6 O(OH) 6 4+

  2

  contribution to anharmonicity resulting from mode correlation appears significant only in the case of ν a (Pb-F) mode, leading to 0.8 cm -1 difference between VSCF and DPT2-VSCF wavenumbers. Analysis of the mode-mode coupling pattern of PbF 2 ••NO vibrational modes,

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

••NO fundamentals (Table 2), whereas

  2

  coupled, which leads in consequence to sizable VSCF anharmonic corrections of the respective fundamentals. Similarity of mode coupling patterns of PbF 2 ••NO and PbF 2 ••ON implies that interaction between mode pairs does not depend substantially on NO coordination mode, but rather involves fundamentals originating in analogous normal mode displacements, like ρ w (PbF 2 ) and ρ t (PbF 2 ) mode pair of both isomers. Same characteristics of the two-dimensional mapping of mode coupling potential pertain to the remaining members of PbX 2 ••NO and PbX 2 ••ON series, so that their coupling patterns are not reported here for the sake of brevity. VSCF computation predicts the most intense overtone and combination transitions of PbF 2••NO and PbF 2 ••ON complexes (Table 2a) to originate in ν Anharmonic spectra of PbCl 2 ••NO and PbCl 2 ••ON complexes Fundamentals of PbCl 2 , computed at MP2/SBKJC+(d) level, are reported in Table3as a reference enabling evaluation of wavenumber shift of ν s (Pb-Cl) and ν a (Pb-Cl) modes resulting from complex formation. Their wavenumbers overestimate experimental Ar matrix spectrum[3] by 12 -29.7 cm -1 due to matrix polarizability effects as well as inaccuracy of the VSCF method itself. Vibrational transitions of PbCl 2 ••NO and PbCl 2 ••ON isomers (Table3) are weakly anharmonic with the exception of ν(Pb-N), ν(Pb-O), ρ w (PbCl 2 ), ρ t (PbCl 2 ), and for the most part in diagonal anharmonicity and mean field / coupling effects with minor contribution from mode-mode correlation. Wavenumber shift of ν(N=O) fundamental on PbCl 2 ••NO complex formation is smaller than that calculated for PbF 2 ••NO species and amounts to 14.2 cm -1 , in accord with reduced σ donation from NO molecule to less positively charged lead ion[3]. Similar effect can be observed in the VSCF spectrum of PbCl 2 ••ON complex, where O-coordination of NO molecule generates oppositely signed shift of ν(N=O) fundamental by -30.7 cm -1 as compared to -33.4 cm -1 in PbF 2 ••ON (Table6). For the same reason the change in energy of ν s (Pb-Cl) and ν a (Pb-Cl) fundamentals on going from PbCl 2 to PbCl

	••ON modes lack significant contribution from mode correlation with the exception of ρ t (PbF 2 ) twisting mode, where DPT2-VSCF fundamental appears 2.8 cm -1 above its VSCF analog. Mapping of PbF F o r P e e r R e v i e w O n l y F o r P e e r R e v i e w O n l y 3.4. Anharmonic spectra of PbBr 2 ••NO and PbBr 2 ••ON complexes The assignment of ν s (Pb-Br) and ν a (Pb-Br) streching modes in the gas phase spectrum of o r P e e r R e v i e w O n l PbBr F y

2 ••ON coupling potentials (Figure

4

), evaluated at ZAPT2/SBKJC+(d) level, indicates that ν s (Pb-F) and ν a (Pb-F), ν(Pb-O) and ρ w (PbF 2 ), ρ w (PbF 2 ) and ρ t (PbF 2 ) mode pairs are strongly s (Pb-F), ν a (Pb-F), ν(Pb-N), ν(Pb-O), ρ w (PbF 2 ), and τ(N=O) modes, which display high intensity fundamentals and are mutually coupled to an appreciable extent (Figures

2 and 3

). Anharmonic corrections of overtones and combination transitions (Table

2a

) almost exclusively stem from diagonal anharmonicity and mean field / coupling effects whereas contribution of inter-mode correlation appears nearly negligible.

3.3.

τ(N=O) modes, whose DPT2-VSCF fundamentals are 4 -10 cm -1 lower in energy than their harmonic equivalents. Anharmonic corrections of PbCl 2 ••NO and PbCl 2 ••ON modes originate 2 ••NO, amounting to -8.9 cm -1 and -8.6 cm -1 , respectively, is much smaller than that computed for the fluoride analog. The same effect is also responsible for the decreased wavenumber shifts of ν s (Pb-Cl) and ν a (Pb-Cl) fundamentals upon PbCl 2 ••ON complex formation. Discernible diminution of the wavenumbers of ν(Pb-N) and ν(Pb-O) fundamentals in PbCl 2 ••NO and PbCl 2 ••ON complexes with relation to fluoride analogs (Table 3) can be attributed to the weakening of Pb-N and Pb-O bonds on fluoride to chloride substitution caused by diminished σ donation from NO molecule, mentioned above. Anharmonic corrections of τ(N=O) torsional mode in PbCl 2 ••NO and PbCl 2 ••ON isomers stem for the most

part from mean field / mode coupling effects and appear coordination mode dependent, leading to 18.6 cm -1 difference between harmonic and VSCF wavenumber of this mode in the case of O-coordinated NO molecule. 2 , placing antisymmetric mode at higher energy than symmetric one

[43]

contradicts the order observed in other PbX 2 halides

[3]

and also that of VSCF computed fundamentals

[4]

(Table

  2

Table 5 .

 5 The red shift of ν(N=O) fundamental on PbI 2 ••ON complex formation (-30.4 cm -1 ) exceeds that calculated in the harmonic approximation (-25.5 cm -1 ), revealing the importance of anharmonic correction in the case of ν(N=O) mode. The wavenumber of ν(Pb-O) fundamental (150.7 cm -1 ) differs from its harmonic counterpart by 8.6 cm -1 due to nearly equal diagonal and mean field / mode coupling contributions to anharmonic correction. The energy of that mode is the lowest in PbX 2 ••ON series, in agreement with the trend observed for other members of the series. VSCF computed shifts of ν s (Pb-I) and ν a (Pb-I) fundamentals upon complex formation, -1.4 and -2.1 cm -1 , respectively, are significantly lower than those evaluated for the other PbX 2 ••ON species, due to dramatically weakened σ-donation from the ligand. In line with the characteristics of PbF 2 ••ON and PbCl 2 ••ON anharmonic spectra, the

	F o r
	P
	e
	e r
	R
	e
	4. CONCLUSIONS Equilibrium geometry parameters of the open-shell PbX 2 ••NO and PbX 2 ••ON species, v i e
	w evaluated at ZAPT2/SBKJC+(d) level, depict both isomers as NO molecule end-on
	coordinated to PbX 2 halide with NO molecular axis nearly parallel to the PbX 2 plane (Figures O 1 and 2). Computation of the anharmonic approximation vibrational spectra of PbX 2 ••NO n l (X=F, Cl, Br) and PbX 2 ••ON (X=F, Cl, Br, I) complexes by DPT2-VSCF method at the same y level of theory afforded satisfactory description of the experimental Ar matrix spectrum of
	PbF 2 ••NO species [3] and provided sensibly consistent predictions of the spectra of the
	remaining compounds. Harmonic approximation based wavenumber shift of the ν(N=O)

••NO complex, only anharmonic spectrum of PbI

2

••ON species is reported in lowest energy τ(N=O) torsional mode fundamental shows exceptionally large anharmonic correction, stemming mostly from intrinsic anharmonicity with smaller contributions arising from VSCF and mode correlation effects. mode upon PbX 2 ••NO complex formation (Table 6) appears qualitatively incorrect as being oppositely signed to the experimental one, whereas that computed from VSCF fundamentals

  2

Table 1 .

 1 Structural parameters of PbX 2 halides and PbX 2 ••L (X=F, Cl, Br, I; L=NO)

	complexes computed at MP2/SBKJC+(d) and ZAPT2/SBKJC+(d) level, respectively.
	Compound	r(Pb-X) [pm]	<(X-Pb-X) [ o ]	r(Pb-L) [pm]	<(Pb-L) [ o ]	r(L-L) e [pm]
		Calcd.	Exptl. Calcd.	Exptl.	Calcd.	Calcd.	Calcd.
	NO							119.2
	PbF 2	205.0	203.6 a	98.0	97.8 c , 98.5 a
	PbF 2 ••NO	206.4		97.3		266.6	109.2	119.0
	PbF 2 ••ON	205.7		98.1		267.8	117.1	119.3
	PbCl 2 PbCl 2 ••NO PbCl 2 ••ON PbBr 2 PbBr 2 ••NO	244.1 F 245.7 o 244.8 258.6 r 244.5 a 257.9 b 260.6	98.3 96±3 d , 100.8 98.6 98.8 99.5 98.8 d 99.9	269.1 270.7 270.7	115.5 124.5 117.8	119.0 119.2 118.9
	PbBr 2 ••ON PbI 2 PbI 2 ••NO PbI 2 ••ON	259.5 279.9 282.1 280.9	100.1 280.7 a 100.6 P 101.1 e r e 101.1	99.7 d	273.5 271.7 275.8	125.5 120.7 127.8	119.2 118.9 119.2
				R
				e
					v i e
					w
						O n l
							y

••ON at ZAPT2/SBKJC+(d) level. Modes are labeled in the order of decreasing harmonic wave number (first column of Table

2

) and the relative strength of the potential is coded in shades of gray (0, white and 255, black). a -r g Ref.

Table 2 .

 2 Vibrational transitions (cm -1 ) of PbF 2 , and PbF 2 ••NO, PbF 2 ••ON complexes calculated at MP2/SBKJC+(d), ZAPT2/SBKJC+(d) level by direct VSCF method.

	38]

Table 2a .

 2a Most intense (I > 1.0 km/mole) overtone and combination transitions (cm -1 ) of PbF 2 , and PbF 2 ••NO, PbF 2 ••ON complexes calculated at MP2/SBKJC+(d), ZAPT2/SBKJC+(d) level by direct VSCF method.

		Harmonic Diagonal	VSCF	DPT2-VSCF Intensity	Assignment a
						[km/mole]
	PbF 2	1078.2	1078.4	1071.8	1069.9	1.2	ν s (Pb-F) + ν a (Pb-F)
	PbF 2 ••NO	1048.2 603.3 581.2 487.2 371.8 341.4 F o r 1048.5 599.9 581.5 471.9 365.8 335.4	1042.6 593.3 570.5 458.1 353.5 325.2	1040.8 592.9 569.7 458.1 353.5 325.2	1.2 2.9 4.7 1.6 2.8 1.1	ν s (Pb-F) + ν a (Pb-F) ν s (Pb-F) + τ(N-O) ν a (Pb-F) + τ(N-O) 2 ν(Pb-Ν) ν(Pb-Ν) + ρ w (PbF 2 ) ip ν(Pb-Ν) + ρ w (PbF 2 ) op
	PbF 2 ••ON	311.8 196.4 136.3 1062.5 618.4 595.8 392.6	305.1 193.6 131.8 P 1062.8 e 611.5 592.5 374.9 e 291.3 181.9 117.8 1057.3 602.4 579.0 357.5 r	291.3 181.9 117.7 1055.4 602.0 578.3 357.5	1.2 4.6 3.3 1.2 2.8 5.3 2.5	ν(Pb-Ν) + τ(N-O) ρ w (PbF 2 ) ip + τ(N-O) 2 τ(N-O) ν s (Pb-F) + ν a (Pb-F) ν a (Pb-F) + τ(N-O) ν a (Pb-F) + τ(N-O) 2 ν(Pb-Ο)
	301.4 275.5 180.9 158.3 155.0 151.7 150.8 a -ν -stretching, δ -deformation, ρ w -wagging, ρ t -twisting, τ -torsional, ip -in-phase 294.2 276.7 276.6 2.3 ν(Pb-Ο) + ρ w (PbF 2 ) ip 268.9 257.0 257.0 1.6 ν(Pb-O) + ρ w (PbF 2 ) op 174.2 156.9 156.9 4.1 R ρ w (PbF 2 ) ip + τ(N-O) 155.5 150.7 150.7 1.3 2 ρ w (PbF 2 ) op 148.8 136.8 136.8 6.0 e ρ w (PbF 2 ) op + τ(N-O) 137.0 116.8 116.8 3.3 2 τ(N-O) 147.0 126.4 126.4 1.2 τ(N-O) + ρ t (PbF 2 ) v i e mode, op -out-of-phase mode w
						O n l
							y

Table 3 .

 3 Vibrational transitions (cm -1 ) of PbCl 2 and PbCl 2 ••NO, PbCl 2 ••ON complexes calculated at MP2/SBKJC+(d), ZAPT2/SBKJC+(d) level by direct VSCF method.

		Mode Harmonic Diagonal VSCF DPT2-VSCF Intensity	Exptl a	Assignment c
							[km/mole]
	PbCl 2	1	349.1	348.2	347.9	347.9	50.5	321.6	ν s (Pb-Cl)
		2	331.1	332.1	329.6	329.6	106.6	299.9	ν a (Pb-Cl)
		3	111.5	111.5	111.0	111.0	3.9	99 b	δ(Cl-Pb-Cl)
	NO	1	1890.1	1869.1			12.4	1875.0	ν(N-O)
	PbCl 2 ••NO	1	1886.0	1883.3 1881.9	1881.9	4.1		ν(N-O)
		2	340.0	339.2	339.0	339.0	47.6		ν s (Pb-Cl)
		3	322.3	323.3	321.0	321.0	99.0		ν a (Pb-Cl)
		4	228.3	223.9	219.5	219.5	106.6		ν(Pb-N)
		5	122.1	121.8	118.2	118.3	0.9		ρ w (PbCl 2 ) ip
		6	107.3	107.2	106.1	106.1	10.5		δ(Cl-Pb-Cl)
		7	79.4	80.1	76.2	76.8	7.2		ρ t (PbCl 2 )
		8	76.6	76.2	74.2	74.1	21.2		ρ w (PbCl 2 ) op
		9	36.7	40.5	34.9	34.1	128.3		τ(N-O)
	PbCl 2 ••ON	1	1864.3	1838.4 1838.7	1838.7	25.6		ν(N-O)
		2	344.9	344.1	344.0	344.0	46.3		ν s (Pb-Cl)
		3	327.3	328.3	326.1	326.1	102.5		ν a (Pb-Cl)
		4	180.6	174.9	170.4	170.4	79.9		ν(Pb-O)
		5	111.6	111.7	111.0	111.0	3.4		δ(Cl-Pb-Cl)
		6	93.9	92.4	88.6	88.8	10.9		ρ w (PbCl 2 ) ip
		7	67.0	68.0	63.4	63.5	8.5		ρ t (PbCl 2 )
		8	63.2	62.8	61.8	61.5	27.1		ρ w (PbCl 2 ) op
		9	53.9	41.8	35.3	35.2	138.6		τ(N-O)

a -Ar matrix data Ref.

[3]

b -gas phase data Ref.

[43]

cν -stretching, δ -deformation, ρ r -rocking, ρ wwagging, ρ t -twisting, τ -torsional, ipin-phase mode, op -out-of-phase mode

Table 4 .

 4 Vibrational transitions (cm -1 ) of PbBr 2 PbBr 2 ••NO, PbBr 2 ••ON complexes calculated at MP2/SBKJC+(d), ZAPT2/SBKJC+(d) level by direct VSCF method.

		Mode Harmonic Diagonal	VSCF DPT2-VSCF Intensity	Exptl.	Assignment c
							[km/mole]
	PbBr 2	1	232.2	231.8	231.6	231.6	22.6	247 a	ν s (Pb-Br)
		2	222.8	223.3	222.1	222.1	60.0	251 a	ν a (Pb-Br)
		3	71.7	71.7	71.5	71.5	1.1	59 a	δ(Br-Pb-Br)
	NO	1	1890.1	1869.1			12.4	1875.0 b	ν(N-O)
	PbBr 2 ••NO	1	1888.6	1883.7	1881.4	1881.3	4.7		ν(N-O)
		2	231.2	230.7	228.6	228.6	55.4		ν s (Pb-Br)+(Pb-N)
		3	216.4	216.8	215.8	215.8	55.7		ν a (Pb-Br)
		4	213.9	213.4	209.0	209.0	65.0		ν(Pb-N)+ν s (Pb-Br)
		5	114.6	113.4	109.8	109.8	4.8		ρ w (PbBr 2 ) ip
		6	72.7	72.7	71.8	71.8	6.6		δ(Br-Pb-Br)
		7	62.0	62.7	59.8	61.1	5.0		ρ t (PbBr 2 )
		8	58.1	58.0	56.9	56.8	8.5		ρ w (PbBr 2 ) op
		9	12.9	31.9	26.5	24.9	119.0		τ(N-O)
	PbBr 2 ••ON	1	1863.9	1841.1	1840.0	1840.0	25.7		ν(N-O)
		2	229.4	228.9	228.5	228.5	21.3		ν s (Pb-Br)
		3	219.8	220.3	219.0	219.0	57.6		ν a (Pb-Br)
		4	170.2	166.9	165.1	165.2	75.6		ν(Pb-O)
		5	90.3	91.0	89.2	89.2	15.4		ρ w (PbBr 2 ) ip
		6	71.1	71.1	70.9	70.9	4.1		δ(Br-Pb-Br)
		7	55.2	57.7	52.8	53.1	77.9		ρ t (PbBr 2 )
		8	48.5	48.3	47.3	47.3	13.9		ρ w (PbBr 2 ) op
		9	47.8	47.3	47.2	47.1	65.3		τ(N-O)

a -gas phase data Ref.

[43]

b -Ar matrix data Ref.

[3]

cν -stretching, δ -deformation, ρ r -rocking, ρ wwagging, ρ t -twisting, τ -torsional, ipin-phase mode, op -out-of-phase mode

Table 5 .

 5 Vibrational transitions (cm -1 ) of PbI 2 and PbI 2 ••ON complex calculated at MP2/SBKJC+(d), ZAPT2/SBKJC+(d) level by direct VSCF method.

		Mode Harmonic Diagonal	VSCF DPT2-VSCF Intensity	Exptl.	Assignment c
							[km/mole]
	PbI 2	1 2	176.7 171.8	176.5 172.1	176.3 171.3	176.3 171.3	13.2 46.0	168 a 163 a	ν s (Pb-I) ν a (Pb-I)
		3	54.3	54.3	54.1	54.1	0.4	43 a	δ(I-Pb-I)
	NO	1	1890.1	1869.1			12.4	1875.0 b	ν(N-O)
	PbI 2 ••ON	1	1864.6	1838.7	1838.9	1838.9	22.4		ν(N-O)
		2	175.5	175.3	174.8	174.9	17.1		ν s (Pb-I)
		3	169.4	169.7	169.2	169.2	45.0		ν a (Pb-I)
		4	159.4	155.7	150.8	150.7	57.0		ν(Pb-O)
		5	85.9	83.4	79.3	79.4	11.8		ρ w (PbI 2 ) ip
		6	54.8	54.8	54.6	54.6	2.9		δ(I-Pb-I)
		7	48.6	48.7	43.8	44.8	51.9		ρ t (PbI 2 )
		8	42.3	42.3	41.6	41.4	7.9		ρ w (PbI 2 ) op
		9	40.5	30.7	25.1	23.9	72.5		τ(N-O)

Table 6 .

 6 Effect of complex formation on wavenumbers of ν(N-O) and ν(Pb-X) stretching modes computed at ZAPT2/SBKJC+(d) level in harmonic (Harm) and DPT2-VSCF (Anh)

	Page 25 of 29	Molecular Physics

URL: http://mc.manuscriptcentral.com/tandf/tmph a -Ar matrix spectra Ref.

ACKNOWLEDGEMENTS

A generous grant of the computer time from the Wrocław Center for Networking and Supercomputing (WCSS) is gratefully acknowledged.


Alemany Diaz

Jorge Boj Viudez

Josefa Mula Bru

C L Esteban

M M E Alemany

J J Boj

J Mula
email: fmula@cigip.upv.es

Francisco C Lario
email: fclario@cigip.upv.es.

Eva María Del Mar

Dpto Alemany

Mathematical programming model for centralised master planning in ceramic tile supply chains

Keywords: MIXED INTEGER LINEAR PROGRAMMING, SUPPLY CHAIN MANAGEMENT, MASTER PRODUCTION SCHEDULE, Master planning, supply chain, ceramic sector, mixed-integer linear programming model Master planning, supply chain, ceramic sector, mixed-integer linear programming model Master planning, supply chain, ceramic sector, mixed-integer linear programming model

published or not. The documents may come   L'archive ouverte pluridisciplinaire

4 described a multi-product, bi-stage MILP model that considered production, distribution, and marketing aspects and also assumed varying time periods for production and distribution. [START_REF] Jayaraman | Planning and coordination of production and distribution facilities for multiples commodities[END_REF] described a multi-product, multi-stage and multi-period MILP model. It included strategic decisions (warehouse and plant localisation) and tactical decisions (product mix, shipments, etc.) and used a heuristic technique based on Lagrangean relaxation to solve the problem. Lavoie and Abdulnour (2003) described an MILP model that considered four first-tier suppliers and a factory to contribute to coordination at a tactical level. This model was multi-periodbased and minimised the total production, inventory and distribution costs. It also used a simulation model to analyse costs at the distribution level by considering vehicle routing, direct shipping, and an intermediate consolidation centre. [START_REF] Lin | A supply network planning problem in a multi-stage and multi-site environment[END_REF] developed a multi-product, multi-stage and multi-period integer linear programming (ILP) model that considered different time periods, i.e., days and months. Wu and Golbasi (2004) proposed an MILP model to solve a multi-item, singlestage, multi-facility lot-sizing problem. It employed Lagrangean decomposition to achieve the model resolution. [START_REF] Park | An integrated approach for production and distribution planning in supply chain management[END_REF] presented an MILP model to jointly solve the production and distribution planning problem in a multi-plant, multi-retailer, multi-item and multi-period logistic environment whose objective was to maximise the total net profit. Although production capacity was assumed to be known, it was necessary to determine the number of vehicles required for distribution from plants to shops. [START_REF] Spitter | Linear programming models with planned lead times for supply chain operations planning[END_REF] proposed an MILP for capacity-constrained supply chain operations planning that used planned lead times with multi-period capacity consumptions where items may be allocated to multiple resources. [START_REF] Schulz | Supply chain optimization of large-scale continuous processes[END_REF] presented a multi-product, multi-stage and multi-period mixed-integer nonlinear programming (MINLP) model applied to a petrochemical SC to consider different types of delivery. [START_REF] Ekşioğlu | Cross-facility management of production and transportation planning problem[END_REF] formulated a uni-product model for transport and production planning as a graph associated with MILP formulation. [START_REF] Ekşioğlu | A Lagrangean heuristic for integrated production and transportation planning problems in a dynamic, multi-item, two-layer supply chain[END_REF] extended this model to a multi-model product and used Lagrangean decomposition to solve it. [START_REF] Lim | A simulation approach for productiondistribution planning with consideration given to replenishment policies[END_REF] presented a hybrid model that comprised not only an MILP model for the decision-making process in the production and inventory capacities in each SC node, but also a simulation model for production and distribution planning. [START_REF] Oh | Global multiproduct production-distribution planning with duty drawbacks[END_REF] presented an LP model for the production and distribution planning of a global chemical company. [START_REF] Bilgen | A mixed-integer linear programming model for bulk grain blending and shipping[END_REF] developed an MILP model for bulk grain blending and shipping. [START_REF] Meijboom | Tactical coordination in a multi-location and multi-stage operations structure: A model and a pharmaceutical company case[END_REF] developed an MILP model for production planning in an SC that was applied to a pharmaceutical company.

With regard to mathematical programming models for ceramic SC planning, [START_REF] Liberatore | A hierarchical production planning system[END_REF] proposed a hierarchical production planning model with two decision levels based on mathematical programming models for planning in a ceramic SC. This model focused on production and distribution stages. It did not consider lot sizing, percentage of first qualities and defects, and part supplying limitations by suppliers. Furthermore, and unlike this article, it considered a single level for the distribution stage. Other works on production planning in the ceramic SC were [START_REF] Özdamar | Simultaneous lot sizing and loading of product families on parallel facilities of different classes[END_REF], [START_REF] Özdamar | Hybrid Heuristics for the capacitated lot sizing and loading problem with setup times and overtime decisions[END_REF], and [START_REF] Özdamar | A hierarchical planning system for energy intensive production environments[END_REF]. These works addressed the problem known as the Capacitated Lot Sizing and Loading Problem (CLSLP) and mainly focused on the production stage without considering the replenishment of raw materials and components (RMs), or the subsequent distribution of finished goods (FGs). After reviewing the literature, we were unable to find studies on mathematical models for the master planning of replenishment, production, and distribution for a ceramic tile SC, which is the main objective of this article.

Problem description

SCs belonging to the ceramic sector present a series of peculiarities, which are described here. Figure 1 provides an overview of the problem resources in terms of replenishment, production, and distribution (physical view). It is assumed that the possibilities of flow between the nodes of the various stages (arcs), as well as the parts, components, raw materials (RMs), and finished goods or articles (FGs) that might circulate through them, have been considered beforehand. The existence Each production plant has one or several production lines (processors in parallel) with a limited capacity. Different FGs can be processed by each production line. There are FGs with high added value that are manufactured only in production plants; others may be partly subcontracted, while some may be totally subcontracted to external suppliers (normally products with a low added value). FGs are grouped into product families for production and commercial reasons. A product family is defined as a group of FGs of identical use (flooring or coverings), format (size), grout (white or red), and whose preparation on production lines is similar. This is done to minimise setup times and costs. Changeovers from one product family to the next incur setup costs owing to the time spent in changing, for instance, moulds [START_REF] Özdamar | Simultaneous lot sizing and loading of product families on parallel facilities of different classes[END_REF]. Lines may not be standardised, in which case, each product family can be processed according to specific facilities with the appropriate technical features. Therefore, not all production lines are capable of processing all the product families, although the product families that may be processed on each line is known. Given the important setup times between product families on production lines, production within a minimum number of consecutive time periods should be carried out whenever a production line is set up for a specific product family (minimum run length). Item setups among the products belonging to the same product family also exist. Because of technological factors involved in the production process itself, when a certain product is manufactured on a specific line, it should be produced in an equal or greater amount than the minimum lot size. This is partly because a certain percentage of defects occurs during the production process, and only a percentage of the manufactured items may be sold as first quality FG.

In the majority of the production planning models developed at the tactical level, the capacities at each stage are aggregated and setup changes are not explicitly considered. However, if at this level the setup times involve an important consumption capacity and have been completely ignored, this 6 may lead to an overestimation of the real capacity availability which, in turn, may lead to unfeasible events during the subsequent disaggregation of tactical plans. Considerable savings may be also be achieved through optimum lot-sizing decisions. However, accounting for setup times at the tactical level would mean simultaneously including decisions about the allocation and lot sizing of production. This problem is known as the Capacitated Lot-Sizing and Loading Problem (CLSLP) [START_REF] Özdamar | Hybrid Heuristics for the capacitated lot sizing and loading problem with setup times and overtime decisions[END_REF]. Given the lengthy setup times involved in the manufacturing of ceramic floorings and coverings, these setup times need to be considered at the tactical level. This work aims to solve this problem within the CLSLP framework.

The distribution of several FGs (multi-item) from production plants to end customers is carried out in various stages (multi-level) by different types of distribution centres (multi-type), such as central warehouses, logistic centres and shops. Neither manufactured nor subcontracted FGs can be stored in manufacturing plants. So they are sent to the first distribution level, which is made up of a number of central warehouses with a limited storage capacity. Outgoing FGs from central warehouses are designed to not only cover the demand of certain end customers (for instance, independent distributors that do not belong to the firm, construction firms, etc.), but to also supply logistics centres. Logistics centres, unlike warehouses, do not have the required storage capacity and only supply FGs to shops that have been previously assigned to them. Finally, shops, which do not have storage capacity, attend to end customers' demands. Although this type of SC attempts to achieve a maximum customer service level, backorders are permitted in both central warehouses and shops. However, backorders quantities are limited to a certain demand percentage to ensure the accomplishment of an objective customer service level defined by the SC. This is a usual situation in the ceramic tile sector, given its limited production flexibility owing to setup costs and times.

Model formulation

A centralised MILP model has been developed to solve the planning problem of replenishing RMs in production plants, producing and subcontracting FGs, and distributing FGs. This model assumes that the decisions made about FGs and product families, which may be manufactured on each line, and those decisions made about establishing distribution routes for the different FGs throughout the network, have been previously made and are known for certain. The objective is to maximise the total net profit over the time periods of the planning horizon by: (i) optimising the total SC network starting with replenishment and going through production on to distribution; and (ii) optimising the production in the production nodes. To accomplish the latter, production planning in plants will have to simultaneously deal with not only the allocation of FGs to production lines with a limited capacity, but also with the determination of lot sizing.

The following MILP formulation, MP-RDSINC, provides a deterministic model for both multisupplier, multi-plant, multi-type, multi-level distribution centres and multi-item and multi-period logistic problems when product demands are known with certainty. The indices, sets of indices, model parameters, and decision variables are described in Tables 1234, respectively. With the proposed definition of the sets of indices, the model is able to capture any flow of materials among the nodes belonging to the various stages of the physical overview (Figure 1).

∑∑ ∑ ∑∑ ∑ ∑ ∑∑ ∑ ∑ ∑∑ ∑ ∑ ∑∑ ∑ ∑∑ ∑ ∑∑ ∑ ∑ ∑∑ ∑ ∑ ∑∑ ∑ ∑ ∑∑ ∑ ∑ ∑∑ ∑ ∑ ∑∑ ∑ ∑ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ - - - - - - - - - - -         + t w Iw(w) i iwt iw 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
∑ ∑ ∑ ∈ ∈ ∈ - - + = ) ( ) ( ) ( 1 ) * ( c Ic i p Lp l ilpt ic c Rc r crpt cpt cpt MP v CTP INC INC c,p,t ∀ (2) cp cpt ssc INC ≥ c,p,t ∀ (3) crt p crpt ca CTP ≤ ∑ Rc(c),t c,r ∈ ∀ (4) ( ) 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 ∈ ∀ (5) ∑ ∈ = ) ( f If i ilpt flpt MP MPF Fl(l),t Lp(p),f l , p ∈ ∈ ∀ (6) t ), l ( Il i ), p ( Lp l , p X * lmi MP ilpt ilp ilpt ∈ ∈ ∀ ≥ (7) t l Il i p Lp l p X M MP ilpt ilpt ), ( ), ( , * 1 ∈ ∈ ∀ ≤ (8) t l Fl f p Lp l p Y M MPF flpt flpt ), ( ), ( , * 2 ∈ ∈ ∀ ≤ (9) t l Il i p Lp l p X X ZI ilpt ilpt ilpt ), ( ), ( , 1 ∈ ∈ ∀ - ≥ - (10) t ), p ( Lp l , p 1 X ZI i ilpt i ilpt ∈ ∀ - ≥ ∑ ∑ (11) t l Fl f p Lp l p Y Y ZF flpt flpt flpt ), ( ), ( , 1 ∈ ∈ ∀ - ≥ - (12) F o r P e e r R e v i e w O n l y 8 t ), p ( Lp l , p 1 Y ZF f flpt f flpt ∈ ∀ - ≥ ∑ ∑ (13) 1 ,.., 1 ' ), ( ), ( , 1 1 ' ' + - = ∈ ∈ ∀ ≤ ∑ - + = flp tmf t t t flpt tmf T t l Fl f p Lp l p ZF flp (14) t ), p ( Ip i , p CTA MP * cq * ) cm 1 ( ) p ( Ap a ipat ) p ( Lp l ilpt i i ∈ ∀ = - ∑ ∑ ∈ ∈ (15) S * sc min CSC ibt ib ) b ( Ab a ibat ≥ ∑ ∈ Bi(i),t PFSP,b i ∈ ∈ ∀ (16) S * sc min CSC ibt ib ) b ( Ab a ibat ≥ ∑ ∈ Bi(i),t PFST,b i ∈ ∈ ∀ (17) S * casc CSC ibt ibt ) b ( Ab a ibat ≤ ∑ ∈ Bi(i),t PFSP,b i ∈ ∈ ∀ (18) S * casc CSC ibt ibt ) b ( Ab a ibat ≤ ∑ ∈ Bi(i),t PFST,b i ∈ ∈ ∀ (19) ia iat ssa INA ≥ Ia(a),t a,i ∈ ∀ (20) a a Ia i iat capal INA ≤ ∑ ∈ ) ( a,t ∀ (21) ∑ ∑ ∈ ∈ - - - + = ) ( ) ( 1 a Qa q iaqt iat a Pa p ipat iat iat CTCL VEA CTA INA INA PFNS, a, t i ∈ ∀ (22) 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 ∈ ∀ (23) CTCL VEA CSC INA INA ) a ( Qa q iaqt iat ) i ( Bi b ) a ( Ba b ibat 1 iat iat ∑ ∑ ∈ ∈ ∧ ∈ - - - + = PFST, a, t i ∈ ∀ (24) iat iat iat iat da DIFA DIFA VEA = - + -1 Ia(a),t a, i ∈ ∀ (25) iat iat da DIFA * 1 α ≤ Ia(a),t a, i ∈ ∀ (26) ∑ ∑ ∈ ∈ = ) ( ) ( q Wq w iqwt q Aq a iaqt CTTK CTCL Iq(q), t q,i ∈ ∀ (27) iwt iqwt VETK CTTK = Iw(w),t Qw(w),i w,q ∈ ∈ ∀ (28) F o r P e e r R e v i e w O n l y 9 Iw(w),t w,i dtk DIFTK DIFTK VETK iwt 1 iwt iwt iwt ∈ ∀ = - + - (29) Iw(w),t w,i dtk * 2 DIFTK iwt iwt ∈ ∀ ≤ α (30) MPF flpt , MP 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, { } 1 , 0 S , ZI , ZF , Y , X ibt ilpt flpt flpt ilpt ∈ (31) T t B, b R, r W, w Q, q A, a P, p L, l C, c I, i F, f ∈ ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈ ∀
The objective function (1) expresses the total net profit over the time periods that have been computed by subtracting total costs from total revenues. Total revenue is included in the income produced by sales in central warehouses and shops. Total costs include transport and purchases costs of RMs, production costs of FGs and setup costs of product families and FGs in production plants, the cost of transport from production plants to warehouses, storage costs, backorder costs in warehouses, costs of subcontracting the FG in the warehouses, transport costs from warehouses to logistics centres, and from logistics centres to shops, and, finally, the costs of backorders in shops.

Constraint (2) is the inventory balance equation of the RM required to manufacture various FGs, of which it forms part. Constraint (3) establishes that the RM inventory does not fall below a minimum value assigned by the safety stock that has been established for each RM, so that a lack of supply does not disrupt production. Constraint (4) ensures that the amount of each RM supplied to all the production plants in each time period does not exceed the capacity of the supply available for each supplier of each RM.

Constraint (5) ensures that the capacity required for the setup of product families and articles in each time period, and the processing of the lots assigned to each line, do not exceed the capacity available on each line at each time period. Constraint (6) ensures that the amount to be produced per product family on each production line in each time period equals the sum of the amounts of the articles belonging to each product family to be manufactured on each particular production line. Constraint (7) guarantees that should a certain amount of an article be manufactured on a production line, it is equal to or above the minimum lot size established for the article on that line. Constraints (8) and ( 9) ensure that only a certain amount of an article or product family could be manufactured on one line, provided that the production of this product or product family, respectively, has been assigned to this line. Parameters M1 and M2 are large enough integer numbers.

Constraints (10)-( 14) guarantee the control of the setup of articles and families. Constraint (10) relates the setup variable (ZI ilpt ) to the production scheduling variables in the current and previous periods (X ilpt and X ilpt-1 , respectively). A changeover to product i on line j occurs in period t, and only if X ilpt =1 and X ilpt-1 =0. If the line was inactive in period t-1, and only a certain amount of i was produced in time period t, the model would calculate a setup for this product. If products were produced on the same production line in time periods t and t-1, Constraint (11) enables saving a single changeover by consecutively sequencing production of one of them, i.e., it is possible to save a changeover for one of the products if its production on the line was done at the end of time period t-1 and continued to the start of time period t. Constraints ( 12) and ( 13) are similar to Constraints (10) and (11), respectively, but refer to product families. Constraint (14) ensures that a production line is setup for a product family at least during the minimum run length of that product family on that particular line (tmf flp ).

As this model focuses only on first quality FGs planning, Constraint (15) ensures that part of the non-faulty production of a given FG, which is first quality, is transported to the central warehouses that can be supplied by this plant. 16) and ( 17) guarantee that the minimum amounts to be subcontracted to external suppliers of partly and completely subcontracted products, respectively, are met should a decision be made to subcontract them during a specific time period. Constraints ( 18) and ( 19) indicate that the subcontracting binary variable takes the value of one when a decision is made to subcontract a certain amount of FGs that can be subcontracted partly or completely, respectively. These constraints also ensure that the amount of FGs subcontracted and transported to warehouses, at each time period and for each FG, that can be subcontracted partly and completely, do not exceed the amount that the corresponding FG supplier has available. As it is possible to meet a certain independent demand through warehouses, Constraint (20) maintains a certain amount of safety stock for each FG to meet any demand variability. Constraint ( 21) is the limitation in the warehouses' capacity that is assumed to be shared by all the FGs. Constraints ( 22)-( 24) are inventory balance equations corresponding to FGs that cannot be subcontracted, those that are usually subcontracted partly and when applicable, and those that can only be subcontracted. As Constraint ( 22) indicates, FGs that cannot be subcontracted and that arrive at the warehouses only come from production plants. Constraint (24) shows how the entry of subcontracted FGs only comes from FG suppliers. Finally, Constraint (23) contemplates the entry of those FGs that can be subcontracted partly in the warehouses, from both production plants and FG suppliers.

As backorders are permitted in both central warehouses and shops, sales may not coincide with the demand for a given time period. Backorder quantities are calculated using Constraint (25). Constraint (26) limits the backorder quantities in each period in terms of a percentage of the demand of each time period. Constraints ( 27) and ( 28) are the inflows and outflows of FGs through each logistic centre and shop, respectively. Constraints ( 29) and ( 30) are similar to Constraints ( 25) and ( 26), but refer to shops.

The model also contemplates non-negativity constraints and the definition of binary variables (31). Finally, some decision variables can be defined as integers, but could change depending on the real-world problem where the model is applied.

The model is easily adapted to cases wherein backorders are not allowed (with a value of zero for both 1 α in Constraint (25) and 2 α in Constraint (30), or by omitting the corresponding backorder decision variables), and when inventory in both the production plants and shops may be held. By defining the different sets, the model is capable of not only representing situations in which all the products may be distributed via all the routes, but also those situations in which this product flow is restricted: for instance, cases in which only certain products may be supplied by certain nodes, and/or the limitations with which the nodes supply others. All this can be achieved without having to introduce additional constraints. Although the minimum run lengths for families (tmf flp ) and minimum lot sizes for manufacturing articles (lmi ilp ) are considered in industry practice, another possible way to extend the model would be to consider minimum lot sizes or run lengths for both.

A detailed computational example

Input data description

We describe in detail an example to validate the described model. The SC physical configuration, cost data, and SC operation characteristics are inspired by a real SC in the Spanish ceramic sector. The model data are based on historical SC information (demand data) and on the mean real values provided by this SC (times and costs). While the physical configuration has been slightly modified for reasons of confidentiality, it respects the structure described in the article. Additionally, we decided to consider an example problem of a size that represents the main relevant sector characteristics, but not excessively large so that it could be described in detail here. Nonetheless, a real industrial problem would be on a larger scale than that considered in this work and would require more data that are beyond the validation purpose of our proposal. The model's planning horizon is assumed to be six weeks (which is usual in the ceramic sector), but excludes the frozen horizon length (usually two/three weeks). In the example considered, 4 FGs, labelled i 1 through to i 4 and 3 FG families, f 1 , f 2 , and f 3 , were included in such a way that i 1 and i 2 belong to f 1 , i 3 belongs to f 3 , and i 4 belongs to f 4 . The FGs, i 1 and i 2 , have the same format (30*30) and are innovative products chosen from the SC catalogue with high added value. They belong to the PFNS set and cannot therefore be subcontracted. The FGs, i 3 and i 4 , have the same format (20*20). The FG, i 3 , belongs to PFSP, so it may be partly subcontracted. i 4 belongs to PFST, meaning that it has to be totally subcontracted. Other relevant FGs data for the model can be consulted in Tables 5678. The FG demand in warehouses (da iat ) and shops (dtk iwt ) was obtained from the SC historical data, and the sale prices in warehouses (pa ia ) and shops (pw iw ) were provided by the SC. Characteristics of central warehouse customers and shop customers differ greatly. Typical examples of central warehouse customers are independent distributors not belonging to the firm and construction firms; meanwhile, those who intend to do house renovations are typical examples of shop customers.

Three production plants exist, p 1 , p 2 , and p 3 , with two production lines per plant, so: l 1 and l 2 belong to p 1 , l 3 and l 4 belong to p 2 , and l 5 and l 6 and belong to p 3 . All the product families may be manufactured on the production lines at the various plants. The production capacity assumed to be available on each line for each time period (caf lpt ) is 50 hours. It is usual in the ceramic sector that production lines work on a two-shift basis of 8 hours per shift (80 hours/week). Because our example only deals with the demand of four FGs, which is obviously only one part of the real demand for the time horizon considered, the production line capacity has been reduced to be more realistic.

The SC considers a mean total production cost of 9 €/m 2 where around 50-60% corresponds to productive costs (direct labour and energy), excluding setup costs, and where around 50-40% corresponds to replenishment costs (RMs and corresponding transport costs). This information has been used to obtain the FG productive costs (costp ilp ), as seen in Table 6. A mean value of an 8hour setup time (tsetupf flp ) is assumed for all the product families, while a 4-hour setup time is assumed for all the FGs within their product family (tsetupi ilp ). In the ceramic tile sector, setups are usually carried out by maintenance personnel. The corresponding setup costs (costsetupf flp , costsetupi ilp ) mainly include the costs of maintenance personnel and the material used for the setup. Table 6 provides the data of the minimum lot size (lmi ilp ) and the unit processing time (tfab ilp ). The minimum product family run length (tmf flp ) is three periods for all the families and lines, and the FG unitary inventory holding cost per period (costina ia ), were set to €0.16 for i1, €0.14 for i2, €0.1 for i3, and €0.09 for i4. All these information were obtained from the ceramic SC. Some of the most common RM of ceramic FGs are glazes and clays. We only consider the RM clay (c1) with a supply cost per kilogram (costtp crp ) of €0.24. The supplier of this RM (r1) is able to supply it to any production plant with a total supply capacity (ca crt ) of 180000 kilograms. A safety stock (ssc cp ) of 2000 kilograms of c1 is taken into account for each production plant. Each m 2 of FG, i1, i2, and i3 will require 12.2, 12.2, and 11.1 kilograms of c1 (v ic ), respectively. There is also an FG supplier (b1) capable of supplying all the products that can be subcontracted (i3 and i4) to all the warehouses up to a limit of 850 m 2 (casc ibt ). There are two central warehouses (a1, a2) with a storage capacity (capal a ) of 200000 m 2 , three logistics centres (q1, q2, and q3) and six shops (w1, w2, w3, w4, w5, and w6). A safety stock (ssa ia ) of 2000 m 2 for i 1 and i 2 and 1600 m 2 for i 3 and i 4 is contemplated for each warehouse assigned. All the production plants may supply all the products to all the warehouses, and all the warehouses may supply any logistics centre. Nonetheless, each 12 logistics centre (q1, q2, and q3) only supplies two shops, i.e., logistics centre q1 supplies shops w1 and w2, logistics centre q2 supplies shops w3 and w4, and logistics centre q3 supplies shops w5 and w6 with any product. Unitary transportation costs for each FG have been calculated based on their corresponding weight and the distance for transport.

The percentage of faults and defects of all the FGs is 10% (cm i ). Likewise, 98% of the production can be sold as first quality (cq i ). Both percentages are mean values provided by the SC. A cost of €100 per m 2 of FG backorder and the time period in warehouses and shops is also contemplated. All initial inventories have been set to their corresponding safety stock level. Details on the data used are presented in Tables 56789101112.

[Insert Table 5 The proposed model has been done in the MPL language, V4.2. The resolution has been carried out with optimisation solver CPLEX 6.6. Finally, the input data and the model solution values were 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.

Evaluation of results

Before evaluating the model, Tables 1314151617show the most representative values of the decision variables that lead to the optimum model solution and used to evaluate the model.

[Insert Table 13 It is worth highlighting that the value of the decision variable sales in the optimum solution of the previous example coincides with the demand data because backorders do not exist in either the central warehouses or shops at any time period.

The assessment method used consisted in analysing the MP-RDSINC model operation according to a set of indicators: computational efficiency, level of service, levels of inventory, and 13 the total net profits in terms of planning. Evidently, any other parameter that is characteristic of the setting in which the method has been applied could be included as an indicator. In particular, such parameters are of interest for the example, which is the objective of the application. All the four indicators selected are explained in detail as follows:

The computational efficiency parameter measures the computational effort required to solve models. The indicators are: the number of iterations needed by the solver and used to reach the optimum solution, the number of model variables, the number of integers in the model, the number of constraints in the model, the number of non-zero elements in the constraints matrix that the model contains, the density of the constraints matrix that the model contains, and the CPU time required to obtain the model solution.

The level of service parameter, expressed as a percentage, is 100% when no delay in demand occurs. If a delay in demand occurs, the level of service could be determined for each planning horizon time period as follows:

Level of Service (%) = 100 1 ×       - eDemand Cummulativ Backorders (30)
To evaluate the inventory parameter levels, the total mean amount of the products inventory (FGs and RMs) for all the planning horizon time periods is calculated.

Alternatively, the inventory maintenance costs are considered in the total net profits parameter. These are the sum of all the incomes less the costs generated in each planning horizon time period considered, i.e., the costs incurred by the production plans provided by the model: replenishment costs, inventory maintenance costs, costs involved in backorders, set up costs involved in the production of a product family and FG, fluctuating production costs, and costs of deliveries among the SC nodes.

Table 18 summarises the computational effort of the model. The Iterations column indicates the number of iterations needed to reach an optimum solution. The number of Variables also appears, these being the Integers, Constraints, and Non-zero elements of the constraints matrix that the model contains. In addition, the Density of the matrix appears. The Solution Time column indicates the time taken to reach a solution in seconds.

[Insert Table 18 about here]

A 100% level of service was achieved with data in the proposed model. Table 19 shows the mean inventory levels of FGs and RM.

[Insert Table 19 about here] Finally in Table 20, the value of the net profit, incomes, and total costs are shown. Given the simplifications assumed in the example, these results are an approximation of the reality. The real SC scenario includes a wide variety of production mix and other additional costs (for instance, the replenishment costs of other RMs). In our example, the model generated a total net profit of € 368,232.45. The different components of the MP-RDSINC objective function appear in Table 20: the incomes are those produced by sales in the central warehouses and shops, replenishment costs include transport costs, and the purchases of RMs. Production costs are those incurred in manufacturing FGs in plants. Inventory costs are the storage costs throughout the SC. Setup costs include those related to the setup of product families and FGs in production plants. Outsourcing costs represent the costs of subcontracting FGs in the warehouses. Transport costs include FGs distribution costs throughout the SC from production plants to warehouses, from warehouses to logistics centres, and from logistics centres to shops. Finally, backorder costs reflect the SC penalty in case of deferred demands in warehouses and shops. The total net profit is calculated as the difference between incomes and total costs.

Conclusions

This work presents the MILP model that has been developed to solve the master planning problem for the replenishment-production-distribution of SCs in the ceramic sector for the purpose of maximising the total net profit in multi-supplier, multi-plant, multi-type, multi-level distribution centres with a multi-item and multi-period logistic environment. The model contemplated was inspired by an actual case and includes the most relevant characteristics of the ceramic sector, e.g., those related to setups, minimum lot sizes for articles, and minimum product family run lengths, faults, and first quality products. This model also considers the several distribution levels integrated by nodes of various types. This model was validated by a real-world case example, but one on a smaller scale for the purpose of providing details of all the input data and of the solution obtained.

Many firms in the ceramic sector only work with planning overviews based on spreadsheets. However, given the increasing complexity of product catalogues and current market pressure to reduce supply times, more rigorous methods are needed to optimise resources. The model presented here offers advantages, such as better coordination and communication among the various departments and SC nodes to substantially lower replenishment, production, and distribution costs, and to improve customer service. Furthermore, given the dramatic increase of end products, the possibilities for assigning and establishing lots on production lines multiply. Therefore, the expected reduction of production plan costs on production lines stands out as the model establishes the product families and articles and the corresponding amounts produced on each line in an attempt to save changeovers as far as possible, this being an important objective in the ceramic sector. This model may also assist managers in evaluating the financial impact of the modification of certain data, e.g., the capacities of plants and warehouses, minimum manufacturing times of product families, and so on.

It is crucial to highlight that demand uncertainty is a major challenge for the ceramic SC in particular and for SCs in general. In this context, we are currently conducting research on managing uncertainty in parameters such as demand, based on our previous works [START_REF] Mula | Models for production planning under uncertainty: A review[END_REF][START_REF] Mula | Capacity and material requirement planning modelling by comparing deterministic and fuzzy models[END_REF]. However, some customers in the ceramic sector only purchase first quality FGs, while others would accept inferior quality. This involves using the nomenclature of first, second, and third quality. Although the production capacity in our proposed model has been developed to take all the qualities into account, the distribution of only first quality FGs has also been planned. One way of extending the model would be to include all the qualities of the products in the proposed model, which would impact the distribution of these products to customers, and to consider not only cost, but also other service parameters. Another alternative way of extending the model would be to consider uncertainty in relation to product quality. We may again consider a centralised view of the problem, which this study has done, because as the analysis was carried out, work was being done with an SC of a firm that formed part of a group of firms working in a centralised manner. However, it would also be interesting to consider a collaborative view of the model and to compare both views. Finally, another line of research would be to study the collaborative decision-making processes between RM suppliers and the ceramic tile SC presented in this article.

F

  Ministry of the Spanish Government, entitled 'Hierarchical methodology in the context of uncertainty in the collaborative planning of a supply-distribution chain/network. Application to the ceramic sector'. Ref. DPI2004-06916-C02-0.

Figure 1 :

 1 Figure 1: Multi-supplier, multi-plant, multi-type-multilevel-distribution centres, multi-item and multi-period logistic model.

Figure 1 :

 1 Figure 1: Multi-supplier, multi-plant, multi-type and multi-level distribution centres, multi-item and multi-period logistic model.

Table 1 .

 1 Indices.

	o
	r
	P
	e e r
	R e v i e w
	O n
	l y

Table 2 .

 2 Sets of Indices

	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

Table 3 .

 3 Model parameters

	Model Parameters
	ca crt	supplier r in period t Capacity (units) of supplying RM c of	M1,M2	Very large integers
	costtp crp	Cost of purchase and transport of one unit of RM c from supplier r to production plant	capal a	Storage capacity (m 2 ) in warehouse a
		p
	caf lpt	Production capacity available (time) of production line l at plant p during time	costtcl iaq	Cost of transporting one m 2 of FG i from 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 during 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 in a time period
		sold as first quality.
	costta ipa	Cost of transporting one m 2 of FG i from production plant p to warehouse a	pa ia	Sales value of one m 2 of FG i in warehouse 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
	costsetup ilp	Setup costs for FG i on production line l of production plant p	1 α	Maximum backorder quantity permitted in a period in warehouses expressed as a
					percentage of the demand of that period
	tfab ilp	Time to process one m 2 of FG i on production line l of production plant p	costsc ib	Cost of subcontracting one m 2 of FG i to FG supplier b
	tsetupf flp	Setup time for product family f on production line l of production plant p	minsc ib	Minimum amount (m 2 ) of FG i to be subcontracted to FG supplier b
	tsetupi ilp	Setup time for article i on production line l of production plant p	costttk iqw	Cost to transport one m 2 of FG i from logistics centre q to shop w
	lmi ilp	Minimum lot size (m 2 ) of FG i on production line l of production plant p	costdiftk iw	Cost to backorder one m 2 of the demand of FG i in a time period at shop w
	tmf flp	Minimum run length (expressed as multiples of the time period used) of	pw iw	Sales price of one m 2 of FG i in shop w
		product family f on production line l of
		production plant p
	v ic	Units of RM c needed to produce one m 2 of FG i	dtk iwt	External demand (m 2 ) of FG i in shop w during the time period t
	ssc 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 in time period t

Table 4 .

 4 Decision Variables

	Decision Variables
	CTP crpt	Amount of RM c to be purchased and transported from supplier r to production plant p in time	INA iat	Inventory (m 2 ) of FG i in warehouse a in time period t
		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 in time period t
	MPF flpt	Amount (m 2 ) of product family f manufactured on production line l of production plant p in time	S ibt	Binary variable with a value of 1 if FG i is subcontracted to supplier b in time period t
		period t
	MP ilpt	Amount production line l of production plant p in time (m 2 ) of FG i manufactured on	VEA iat	Amount (m 2 ) of FG i sold in warehouse a during time period t
		period t
	X ilpt	Binary variable with a value of 1 if FG i is manufactured on production line l of production	DIFA iat	Backorder quantity (m 2 ) of FG i in warehouse a during time period t
		plant p in time period t, and with a value of 0
		otherwise
	Y flpt	Binary variable with a value of 1 if product family f is manufactured on production line l of	CTCL iaqt Amount (m 2 ) of FG i transported from warehouse a to logistics centre q in time
		production plant p in time period t, and with a		period t
		value 0 otherwise
	ZI ZF flpt	Binary variable with a value of 1 if a setup takes place of product family f on production line l of	VETK iwt	Amount (m 2 ) of FG i sold in shop w during time period t
		production plant p in time period t, and with a
		value of 0 otherwise
	CTA ipat	Amount (m 2 ) of FG i to be transported from production plant p to warehouse a in time period	DIFTK iwt Backorder quantity (m 2 ) of FG i in shop w during time period t
		t

ilpt Binary variable with a value of 1 if a setup takes place of product i on production line l of production plant p in time period t, and with a value of 0 otherwise CTTK iqwt Amount (m 2 ) of FG i transported from logistics centre q to shop w in time period t

Table 5 .

 5 Data of FG in warehouses

	daiat

Table 6 .

 6 Data of FG in plants

	Page 21 of 28		International Journal of Production Research
	Family	FG	Line	Plant	costpilp	costsetupilp	tfabilp	lmiilp	costsetupfflp
	f1	i1	l1	p1	5,76	459,43	0,021	1500	930,46
			l2	p1	5,76	459,43	0,019	1500	930,46
			l3	p2	5,76	459,43	0,019	1500	930,46
			l4	p2	5,76	459,43	0,021	1500	930,46
			l5	p3	5,76	459,43	0,019	1500	930,46
			l6	p3	5,76	459,43	0,019	1500	930,46
		i2	l1	p1	5,28	351,16	0,015	1500	930,46
			l2	p1	5,28	351,16	0,017	1500	930,46
			l3	p2	5,28	351,16	0,017	1500	930,46
			l4	p2	5,28	351,16	0,017	1500	930,46
			l5	p3	5,28	351,16	0,015	1500	930,46
			l6	p3	5,28	351,16	0,015	1500	930,46
	f2	i3	l1	p1	4,8	322,29	0,014	1500	798,32
			l2	p1	4,8	322,29	0,012	1500	798,32
			l3	p2	4,8	322,29	0,014	1500	798,32
			l4	p2	4,8	322,29	0,014	1500	798,32
			l5	p3	4,8	322,29	0,012	1500	798,32
			l6	p3	4,8	322,29	0,014	1500	798,32

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.uk

Table 7 .

 7 Data of FG in shops

	dtkiwt

Table 8 .

 8 Data of FG in suppliers

	FG	Supplier	costscib	minscib
	i3	b1	10	100
	i4	b1	7,1	100

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.

uk International Journal of Production Research

	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y
	23

Table 9 .

 9 FG transport costs (1)

	FG	Plant	Warehouse costtaipa
	i1	p1	a1	0,054
			a2	0,05
		p2	a1	0,05
			a2	0,054
		p3	a1	0,088
			a2	0,08
	i2	p1	a1	0,054
			a2	0,05
		p2	a1	0,05
			a2	0,054
		p3	a1	0,088
			a2	0,08
	i3	p1	a1	0,054
			a2	0,05
		p2	a1	0,05
			a2	0,054
		p3	a1	0,088
			a2	0,08

Table 10 .

 10 FG transport costs (2)

	FG	Warehouse Logistics Centre	costtcliaq	FG	Warehouse Logistics Centre	costtcliaq
	i1	a1	q1	0,305	i3	a1	q1	0,77
			q2	1,06			q2	2,7
			q3	1,24			q3	3,15
		a2	q1	0,305		a2	q1	0,77
			q2	1,06			q2	2,7
			q3	1,24			q3	3,15
	i2	a1	q1	0,305	i4	a1	q1	0,77
			q2	1,06			q2	2,7
			q3	1,24			q3	3,15
		a2	q1	0,305		a2	q1	0,77
			q2	1,06			q2	2,7
			q3	1,24			q3	3,15

Table 11 .

 11 FG transport costs (3)

	FG	Logistics Centre	Shop	costttkiqw	FG	Logistics Centre	Shop	costttkiqw
	i1	q1	w1	0,081	i3	q1	w1	0,206
			w2	0,09			w2	0,229
		q2	w3	0,072		q2	w3	0,183
			w4	0,081			w4	0,206
		q3	w5	0,068		q3	w5	0,172
			w6	0,09			w6	0,229
	i2	q1	w1	0,081	i4	q1	w1	0,206
			w2	0,09			w2	0,229
		q2	w3	0,072		q2	w3	0,183
			w4	0,081			w4	0,206
		q3	w5	0,068		q3	w5	0,172
			w6	0,09			w6	0,229

Table 12 .

 12 RM transport costs

	RM	Supplier	Plant	costtpcrp
	c1	r1	p1	0,24
			p2	0,24
			p3	0,24

Table 13 .

 13 Amount (m 2 ) of product i manufactured on production line l of production plant p in time period t

	Mpilpt

Table 14 .

 14 Amount (m 2 ) of subcontracted FG i subcontracted to supplier b for warehouse a in time period t

	CSCibat

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.

uk International Journal of Production Research

	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y
	25

Table 15 .

 15 Amount of the RM c to be purchased and transported from supplier r to production plant p in time period t

	CTPcrpt

Table 16 .

 16 Inventory (m 2 ) of FG i in warehouse a in time period t

	INAiat

Table 20 .

 20 Total net profit

		MP-RDSINC
	Incomes	1059935,65
	Replenishment costs	208156,59
	Production costs	381360,54
	Inventory costs	11259,79
	Outsourcing costs	36024,3
	Set-up costs	6463,63
	Transport costs	48438,35
	Backorder costs	0
	Total costs	691703,2
	Profits	368.232,45 €

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

project funded by the Science and Technology Ministry of the Spanish Government, entitled "Hierarchical methodology in the context of uncertainty in the collaborative planning of a supply-distribution chain/network. Application to the ceramic sector". Ref. DPI2004-06916-C02-0. project funded by the Science and Technology