Stochastic model for heterogeneous rail freight car fleet management based on the model predictive control

CHAPTER 6

Stochastic model for heterogeneous rail freight car fleet management based on the model predictive control Contents 6.1 Discrete time MPC framework for rail freight car fleet sizing and allocation problem 6.2 Design variables 6.3 System performance measure and constraints 6.3.1 System dynamics 6.3.2 Nonnegative control constraints 6.3.3 Nonnegative state constraints 6.3.4 Capacity constraints 6.4 State vector estimation 6.5 Forecasting the state of rail freight cars by ARIMA-Kalman method 6.5.1 Components of weighted matrices A, B, and L 6.5.2 The components of matrix Γ(P) 6.6 MPC controller 6.6.1 Optimization problem 6.6.2 Detailed description of the MPC approach 6.6.3 Numerical experiments

174 175 181 182 183 183 183 184 185 188 196 198 198 199 204

From the literature review, it can be seen that the previous contributions were mainly concentrated on a homogeneous rail freight car fleet and/or empty flows in a static deterministic environment. Namely, the approaches were based on a multiperiod integer linear programming or mixed integer linear programming formulation that maximize the profit or minimize the cost of the rail freight car fleet system and these approaches are valid only if the demands, forecasts, disturbances, and parameters are constant overtime. This is never the case, so it is necessary to find an optimal policy in a dynamic environment. For that purpose, in this chapter, a model predictive control (MPC) framework was developed for obtaining an optimal set of control actions that involve state, control, and capacity constraints. The demand Optimization Models for Rail Car Fleet Management https://doi.org/10.1016/B978-0-12-815154-9.00006-X

© 2020 Elsevier Inc. All rights reserved.

173

174

Optimization models for rail car fleet management

as a stochastic variable and heterogeneous rail freight car fleet with a partial substitutability among various freight car types are considered.

6.1 Discrete time MPC framework for rail freight car fleet sizing and allocation problem MPC represents a multivariable control method that is usually used to optimally control complex systems while explicitly considering constraints. The system dynamics is represented by a discrete time predictive model where the next state is a function of the current state, current demands (disturbance), and the current control vector with constraints. In each time instance, the optimal control problem is solved online based on the current state (at time n ¼ n0) and the predictive demands over a N step finite horizon (at time n ¼ n0, …, n0 + N  1). The result of the optimization is a sequence of control vectors overtime n ¼ n0, …, n0 + N  1 but only the first control vector (at time n ¼ n0) is applied to the system. In the next state (at time n ¼ n0 + 1), the optimal control problem is solved again for the time horizon n ¼ n0, …, n0 + N then only the control vector at time n ¼ n0 + 1 is applied and so forth. In this chapter, a general MPC framework for the simultaneous rail freight car fleet size and allocation problem is considered (Milenkovi
c et al., 2015a, b). The proposed methodology is graphically illustrated in Fig. 6.1. A linear dynamic model of the system considers loaded and empty freight car flows of different types. Partial substitutability between freight car types

Fig. 6.1 Rail freight car fleet size and allocation approach.

Stochastic model for heterogeneous rail freight car fleet management

175

has also been covered by the model. State vector estimation is based on the forecasts of the number of freight cars of each type in each station and the dynamic model of the system. Forecasts are obtained by the autoregressive integrated moving average (ARIMA)-Kalman forecasting method. The Kalman filter method is used for deriving the optimal system state over the prediction horizon. For the evaluation of rail freight car inventories by stations, the stochastic multiperiod economic order quantity (EOQ) model and single period random newsboy model are applied. Based on the optimal state estimate, the objective functional, and constraints, the MPC controller is designed and used for solving the resulting rail freight car fleet sizing and allocation problem. Presented approach demonstrates the applicability of the receding horizon control for determining the size of rail freight car fleet based on establishing a balance between exact cost parameters, the cost of unmet demand on one side, and the sum of car ownership and utilization costs on the other. All these costs are given as input parameters for each freight car type. The mathematical model enables through the vector difference state equation, an adequate presentation of the dynamic nature of freight car demand. The proposed optimal MPC controller will be applied for deriving the optimal rail freight car fleet investment decisions.

6.2 Design variables As in the case of the fuzzy random modeling approach, let’s consider a network composed of a set of stations N. The planning horizon P has been divided into discrete time intervals 0, 1, …, m, …, n, …, P. This model takes into account the possibility of substitution of the desired freight car type by an appropriate alternative. Namely, in some cases, a transport of specific commodities can be performed by more than one type of freight cars. In that case, the demand may be specified for an aggregate car type. The aggregate car type is only used in distribution planning, there are no cars of aggregate car type. The demand for an aggregate car type is fulfilled only by original car types. The aggregation of car types is prespecified, that is, which original car types can be used to satisfy the demand for an aggregate car type ( Joborn, 2001). The total number of freight car types is denoted as Tσ , composed of the original types T(T  Tσ ) and aggregate car types Ta(Ta  Tσ ), respectively. For each, original and aggregate car type, t (t 2 Tσ ) requests for rail transportation service exist between stations

176

Optimization models for rail car fleet management

i and j in period n, Dijt(n). Dijt(n) represents the demand for freight car of type t from station i to station j in period n, and it is assumed to have a stochastic character. The demand generates loaded Fijt(n) as well as empty freight car flows Eijt(n) on the network. Control actions which represent the set of all loaded and empty freight car departures are under the influence of freight car arrival times. As in the case of the fuzzy random model, rather than to model the traveling times directly, the problem is formulated in terms of freight car arrivals. More precisely, given that Fijt(m) freight cars of type t (t 2 T) were dispatched from station i in the period m, in this approach the question is how many of these cars actually arrive at station j in the period n. Therefore, θijt(m, n)[αijt(m, n)] are defined as the proportions of loaded [empty] freight cars of type t dispatched from i to j in the period m which actually arrive in the period n. These proportions are random variables since the information about freight car arrival is not known in advance. In order to provide an adequate response regarding new freight car requests in a very uncertain environment, it is necessary to maintain a pool of freight cars in the stations. In this approach, variable Sit(n) represents the number of freight cars of type t (t 2 T ) present in the station at the end of period n. The system state represents the total number of a certain type t (t 2 T ) of freight cars offered at each station for each day over the planning horizon. In case if insufficient freight cars of type t are available at the station i in period n to meet all demand, some demand will be either backordered or lost from the system. Variable Vijt(n) represents the number of units of demand, or freight car loads of type t (t 2 Tσ ) that remains unmet at the end of the period n. The notation for the model formulation is given in Table 6.1. The basic formulations between the defined quantities can be stated as follows: Sit ðn + 1Þ ¼ Sit ðnÞ +

N X X  X

Fjitt0 ðmÞ  θjit ðm, n + 1Þ + Ejitt0 ðmÞ  αjit ðm, n + 1Þ



j¼1 t0 2T t m


N  XX

 Fijtt0 ðn + 1Þ + Eijtt0 ðn + 1Þ , t ¼ 1,…, T

t 0 2T t j¼1

(6.1) Vijt0 ðn + 1Þ ¼ Vijt0 ðnÞ + Dijt0 ðn + 1Þ 

X

Fijtt0 ðn + 1Þ, t 0 ¼ 1,…, T σ

(6.2)

t2T t0

Fijtt0 ðnÞ, Eijtt0 ðnÞ, Vijt0 ðnÞ, Sit ðnÞ  0

(6.3)

Stochastic model for heterogeneous rail freight car fleet management

177

Table 6.1 Notation for model formulation Notation

Definition

Sets

N P Po Tσ T Ta Tt Tt

Set of stations in a railway network Set of periods over a planning horizon Set of observations in a time series of the number of rail freight cars by stations Set of all original and aggregate car types Set of all original car types, T  Tσ Set of all aggregate car types, Ta  Tσ Set of car types, original and aggregate, that can be substituted by the type t. If type t is an original type, then the type itself is included in T t Set of all car types that can substitute for the type t. If type t is an original type, then the type itself is included in T t . Thus, T t are all car types that can be used to fulfil a demand for cars specified for type t

Parameters

Dit Dijt(n) ωijt(n) λijt(n) hit pijt eijt lijt rijtt0 Qt θijt(m, n)

Total demand for freight cars of type t at station i in the first P  1 periods of the planning horizon Random demand for a rail freight car of type t between station i and station j over period n, i, j 2 N; t 2 Tσ ; n 2 P, i 6¼ j Random component of the demand for cars of type t between stations i and j during the period n, i, j 2 N; t 2 Tσ ; n ¼ 0, …, P  1, i 6¼ j Exact component of the demand for cars of type t between stations i and j during the period n i, j 2 N; t 2 Tσ ; i 6¼ j Unit holding cost for a rail freight car of type t at station i per period, i 2 N; t 2 T Unit shortage cost for a rail freight car of type t between stations i and j per period i, j 2 N; t 2 T; i 6¼ j Unit cost of the empty movement of rail freight cars of type t from station i to station j, i, j 2 N; t 2 T; i 6¼ j Unit cost of the loaded movement of rail freight cars of type t from station i to station j, i, j 2 N; t 2 T; i 6¼ j Cost of substituting a car of type t0 by a car of type t, i, j 2 N ;t 2 T ;t 0 2 T t ;n 2 P;i 6¼ j Unit ownership cost for a car of type t traveling between stations, per period Portion of loaded cars of type t arrived during the period n, and dispatched during the period m from station i to station j, n > m, i, j 2 N; t 2 T; n ¼ 0, …, P  1, i 6¼ j Continued

178

Optimization models for rail car fleet management

Table 6.1 Notation for model formulation—cont’d Notation

Definition

αijt(m, n)

Portion of empty cars of type t arrived during the period n, and dispatched during the period from station i to station j, n > m, i, j 2 N; t 2 T; n ¼ 0, …, P  1, i 6¼ j Number of loaded cars of type t dispatched from station i during the period m and arrived at station j during the period n, i, j 2 N ;t 2 T ;t 0 2 T t ;n > m,i 6¼ j Number of empty cars dispatched from station i during the period m and arrived at station j during the period n, i, j 2 N ;t 2 T ;t 0 2 T t ;n > m,i 6¼ j Forecasting error of the number of cars of type t at station i during the period n, i 2 N; t 2 T; n ¼ 0, …, P Ratio of the number of empty cars of type t to the total number of freight cars Random variable describing travelling time for freight cars of type t from station i to station j i, j 2 N; t 2 T; i 6¼ j

xijtt0 m(n)

yijtt0 m(n) Δit(n) ϑt τijt

Objective functional & decision variables

J ϕ1(X(P)) ϕ2[X(n), U(n), n]

Sit(n) Vijt(n) Eijtt0 (n)

Fijtt’(n) γ it

System performance measure, which represents the optimality criterion for dispatching empty and loaded cars. In this paper it is also denoted as the cost functional Function of state vector in the last period of the planning horizon derived from the unit holding cost, unit shortage cost and car ownership cost during travel Function of state and control vectors during the period n; it is derived from the unit holding cost, unit shortage cost, car ownership cost during travel and unit cost of moving empty and loaded cars from station to station Number of empty and loaded cars of type t at station i at the end of the period n, i 2 N; t 2 T; n ¼ 0, …, P Unmet demand for freight car loads of type t between stations i and j over the period n, i, j 2 N; t 2 Tσ ; n ¼ 1, …, P, i 6¼ j Number of empty car flows of type t used to fulfil demand for car type t0 between station i and station j over the period n, i, j 2 N ;t 2 T ;t 0 2 T t ;n ¼ 1, …,P, i 6¼ j Number of loaded car flows of type t used to fulfil the demand for car type t0 between station i and station j over the period n, i, j 2 N ;t 2 T ;t0 2 T t ;n ¼ 1, …,P, i 6¼ j; The average freight car of type t inventory status and shortage at the station i over the planning horizon per period considered, i 2 N; t 2 T;

Stochastic model for heterogeneous rail freight car fleet management

179

Table 6.1 Notation for model formulation—cont’d Notation

Definition

ηit

The average number of (loaded and empty) cars dispatched from the station i over the planning horizon per period considered i 2 N; t 2 T;

Vectors and matrices

X(n) U(n) d(n) Z(n) ν(n) Λ(n) G Φ M(0)

R1(n) R2(n) H A(n) B L Γ(P)

State vector with components Sit(n), Vit(n), Dijt(n), xijtt0 m(n), yijtt0 m(n) Control vector with components Eijtt0 (n), Fijtt0 (n) Vector of disturbances Vector of the forecasted number of freight cars by station during the period n ¼ 0, …, P Vector of error in the forecasted number of freight cars by station during the period n ¼ 0, …, P State transition matrix represented in terms of unit matrices, zero matrices and matrices Λ1, Λ2(n), Λ3(n), and Λ4 System control matrix System disturbance matrix Symmetric, positive-semidefinite matrix of covariances of the state vector at the initial moment and of the preceding state estimation for the same period Symmetric, positive-semidefinite matrix of covariances of the random component of demand Symmetric, positive-semidefinite matrix of covariances of white noise in the forecast Matrix in the linear term of the state vector in the expression for the forecasted number of freight cars by stations Positive-definite matrix with components π sit and P P qt Ps¼m saijts ðn + s  mÞ and qt Ps¼m sbijts ðn + s  mÞ, t 2 T Positive-definite matrix with components π uit, t 2 T Matrix with components π hit, and rijtt0 t 2 T, t0 2 Ta Positive-semidefinite matrix in the quadratic term of the state vector in the last period of the planning horizon

Relation (6.1) represents conservation of flow constraints for each type of freight cars that include the effects of uncertain travel times for car movements through random θ and α terms. The number of freight cars of type t in the next period is equivalent to the number from the last period plus loaded and empty freight car flows which will arrive during the n + 1 period minus the dispatched loaded and empty freight car flows of type t used to fulfil the demand for type t as well as the types t0 which can be substituted by cars of type t. Relation (6.2) exists for all original and aggregate car types.

180

Optimization models for rail car fleet management

The demand for t0 can be fulfilled by any car types t 2 T t0 that can substitute for the type t0 . Unmet demand in the period n + 1 is equal to the unmet demand from the previous period plus the new demand minus the cars that are used to fulfil demand in the time period n + 1. The last relation reflects the natural assumption of nonnegativity of freight car flows, states of cars in every station, and unmet demands. The demand Dijt(n) (t ¼ 1, …, Tσ ) is modeled as a sum of random (stochastic) and deterministic components. This sequence is assumed to be of Gauss-Markovian character. Considering that the Gaussian form is retained in the case of linear transformation, it can be described by a state vector of a linear dynamic system excited by a purely Gaussian random sequence (uncontrollable system input), with initial Gaussian state vector. A Gaussian pure random process represents the boundary value of a Gauss-Markovian process with a very large covariance and very short correlation time. As the Fourier transformation of the correlation function with respect to time is a constant, its spectrum is white (Bryson and Ho, 1975). That is why a pure Gaussian random process is often referred to as white noise. Let us introduce an auxiliary variable dijt(n) for all car types for which the following holds: dijt ðnÞ ¼ Dijt ðn + 1Þ

(6.4)

dijt ðn + 1Þ ¼ λijt ðnÞλijt ðn  1Þ⋯λijt ð0Þdijt ð0Þ + ωijt ðnÞ, n ¼ 0,1,…, P  1 (6.5) λijt ðnÞ ¼

μijt ðn + 2Þ , n ¼ 0,1, …,P  2 μijt ðn + 1Þ   E ωijt ðnÞ ¼ 0

(6.6) (6.7)

Relations (6.4)–(6.6) enable converting any random process to an equivalent Markov random process by a proper enlargement of the state vector for variable dijt(n). Relation (6.7) represents the expectation of the Gaussian purely random process, which is the random component of demand in this case. In order to determine random variables θijt(m, n) and αijt(m, n), we applied Liou’s numerical procedure for evaluating a transient system response (Liou, 1966). More precisely, a certain period of time must elapse between the dispatch of a freight car from one station and arrival at another and therefore, the system’s response to a specified control action is delayed in time. As the dynamics of car arrivals is a stochastic function of the delayed system response to a specified demand, it is necessary to know its probability density distribution. A number of distributions may be used for modeling uncertainty in rail freight car traveling times (Spiegelman et al., 2010). For the

Stochastic model for heterogeneous rail freight car fleet management

181

purpose of this study, travel times of freight cars are assumed to follow a negative binomial distribution ( Jordan and Turnquist, 1983; Bojovic, 2002). In that context, this distribution can be interpreted for describing the delay time that has to elapse before a car dispatched from station i appears at station j. Let us now reformulate θijt(m, n) and αijt(m, n) as follows: ajit, n + 1m ðnÞ ¼ θjit ðm, n + 1Þ, m < n , n ¼ 0,…, P,t 2 T bjitt0 , n + 1m ðnÞ ¼ αjit ðm, n + 1Þ, m < n , n ¼ 0, …,P,t 2 T

(6.8) (6.9)

By associating another index, the variables that denote the arrivals of loaded and empty cars from preceding time intervals can also be described: xjitt0 , n + 1m ðnÞ ¼ Fjitt0 ðmÞ, m < n , n ¼ 0,…, P,t 2 T, t 0 2 T t yjitt0 , n + 1m ðnÞ ¼ Ejitt0 ðmÞ, m < n , n ¼ 0,…, P,t 2 T, t 0 2 T t

(6.10) (6.11)

A pair of additional variables that are defined as follows: fijtt0 ðnÞ ¼ Fijtt0 ðn + 1Þ, eijtt0 ðnÞ ¼ Eijtt0 ðn + 1Þ, 8i , j, t,t 0 , n

(6.12)

has also been introduced.

6.3 System performance measure and constraints The performance measure in case of stochastic rail freight car fleet sizing and allocation model has the same structure as in the “fuzzy” case, except that it contains only crisp values and it contains cost terms related to different types of freight cars, as well as to substitution possibilities. Therefore, the crisp form of the objective functional looks as follows: J ¼ ϕ1 ½X ðP Þ +

P1 X

ϕ2 ½X ðnÞ, U ðnÞ, n

(6.13)

n¼0

The first term is a cost of the rail freight car fleet system in the last period over the planning horizon which is the function (ϕ1) of the state vector as a result of control actions during the previous periods. It is derived based on the same cost parameters but defined for each type of freight cars independently: the unit cost of holding a freight car of type t for one period in a station i, hit, the unit shortage cost or penalty cost per period for one freight car of unmet demand of type t from station i to station j, pijt, and the daily cost of ownership or leasing a freight car of type t while traveling on route, Qt. The cost component Qt is determined per car per period (Beaujon and Turnquist, 1991). The holding cost for a non-traveling freight car of type t is at least the ownership cost Qt, but in general it would include additional cost associated with storage and management of freight cars in a station.

182

Optimization models for rail car fleet management

In this paper, hit is defined as the unit cost of holding a freight car of type t for one period at a location, where hit Qt. The second term represents a cost of rail freight car fleet system in P  1 periods over the planning horizon which is the function (ϕ2) of the state vector X(n) and control vector U(n) and it is derived from the unit holding cost, hit, the unit shortage cost, pijt, the cost of substituting a car of type t0 by a car of type t, rijtt’ (t 2 T ,t 0 2 T t ) and the unit ownership or leasing cost during travel, Qt, and the unit costs of moving empty, ejit, and loaded, ljit, cars from one station to another at the remaining discrete instants. MPC-based formulation contains a set of constraints that restrict the feasibility region.

6.3.1 System dynamics The model is described by a time discrete control system whose states vary according to the following difference equation: X ðn + 1Þ ¼ Λ½n, aðnÞ, bðnÞX ðnÞ + GU ðnÞ + ΦdðnÞ, n ¼ 0,1,…, P  1 (6.14) Therefore, the state in the next period is a function of the state, control, and vector of disturbances d(n). Control vector U(n) contains control actions of empty and loaded freight cars dispatched from station i in the period n for all original and aggregate car types.   fijtt0 ðnÞ U ðnÞ ¼ X , i, j 2 N, t 2 T , t 0 2 T t (6.15) eijtt0 ðnÞ 2NT ðN 1Þ + 2N ðN 1Þ T t0 t’2T a

Matrix Λ[n, a(n), b(n)] is a state transition matrix whose dimensions are  X  ½TN ð2N  1Þ + 2N ðN  1Þ T a + ðP  1ÞðT + T tÞ  a t2T  X  T tÞ   ½TN ð2N  1Þ + 2N ðN  1Þ T a + ðP  1ÞðT + t2T a

Matrix "

G

matrix with dimensions ! !# X a  ½2NT NT ð2N  1Þ + 2N ðN  1Þ T + ðP  1Þ T + Tt t2T a X ðN  1Þ + 2N ðN  1Þ T t Φ is the disturbance matrix defined as t2T a " a

unit

represents

matrix

of

an

input

dimensions

NT ð2N  1Þ + 2N ðN  1Þ

183

Stochastic model for heterogeneous rail freight car fleet management

T a + ðP  1Þ T +

X

!!# Tt

 ½NT ð2N  1Þ + 2N ðN  1ÞðT a +

t2T a

ðP  1Þ T +

X

!!# d(n) is the vector of disturbances which represent

Tt

t2T a

a random component of demand for cars during the period n with components ωijt(n). The state vector X(n) is of dimensions ½NT ð2N  1Þ + 2N ðN  1Þ !!# X T a + ðP  1Þ T + and describes the current state in the rail Tt t2T a

freight car system through the number of loaded and empty freight cars in the period n, the empty and loaded flows arriving in all stations from the preceding time periods and unmet demand over the same period for the original and aggregate car types.

6.3.2 Nonnegative control constraints For each period n ¼ 0, 1, …, P  1 and each car type t 2 T, the number of loaded and empty freight car flows is a nonnegative number: U ðnÞ  O2 where O2 represents a 2NT ðN  1Þ + 2N ðN  1Þ

X

(6.16) T t zero matrix.

t2T a

6.3.3 Nonnegative state constraints For each period n ¼ 0, 1, …, P  1 and car type t 2 T, the current state in the rail freight car system must be a nonnegative number: X ðnÞ  O3

(6.17)

"

where O3 represents a

NT ð2N  1Þ + 2N ðN  1Þ T a + ðP  1Þ

T+

X

!!# Tt

t2T a

zero matrix.

6.3.4 Capacity constraints We also impose the constraints in which the number of rail freight cars in each station will be limited by the capacities of stations: DX ðnÞ  K

(6.18)

where D is a matrix ½N  ðNT ð2N  1Þ + 2N ðN  1Þ  a P of dimensions T + ðP  1Þ T + t2T a T t  that contains the length of rail freight cars

184

Optimization models for rail car fleet management

of certain types. K represents a column vector of dimensions [N  1] with station capacities or maximum numbers of freight cars on the state in a station during every period of the planning horizon.

6.4 State vector estimation Two components are crucial for the estimation of the state vector—the dynamic model of the freight car fleet system (6.14) as well as the forecasted number of cars of each type t, Zt(n), at discrete time instants. Planning of the freight car inventory level in every station on the network requires forecasting of their levels in future time periods. These values represent specific state observations, X(n). The observation equation in analytical form looks like the following: Zit ðnÞ ¼ Sit ðnÞ + Δit ðnÞ, i 2 N , t 2 T , n ¼ 0, 1,…,P

(6.19)

The first term Sit(n) denotes the forecasted number of freight cars of type t at the station i in the time period n whereas the second term Δit(n) is a random error of the forecasts for the same station and period. Various methods can be applied for deriving the forecasts of the freight car levels. In this chapter, the ARIMA-Kalman forecasting model is developed and applied. The details of this technique follow after this section. The forecasting error is determined as the mean value of the differences between the actual and forecasted value of the number of freight cars by stations. The forecasted number of cars in the initial period is assumed to be equal to the mean value of the initial state vector which is calculated by averaging over all periods for which historical data are available. Having in mind that the process of variation in the status of freight cars by periods is relatively stable, this estimation enables an appropriate initial value. A forecasted value or an observation can be represented as a linear function of the state: Z ðnÞ ¼ HX ðnÞ + vðnÞ, n ¼ 0,1,…, P

(6.20)

where follows that 2

Z ðnÞ ¼ ½Z1 ðnÞ, …, ZT ðnÞT

H ¼ 4 diag½Ht O

 2N ðN 1Þ T a + ðP1Þ

 5, t 2 T

X t2T

T a

vðnÞ ¼ ½v1 ðnÞ, …, vT ðnÞT

3

(6.21) (6.22)

t

(6.23)

Stochastic model for heterogeneous rail freight car fleet management

185

Ht is a matrix of dimensions N  [N(2N  1) + 2N(N  1)(P  1)] and vt(n) is a N-component vector of the error in the forecasted number of freight cars by station during period n. vt(n) is a N-component vector of the error in the forecasted number of freight cars by station during period n. Vector vt(n) is the white noise, which represents the forecasting errors for cars of type t. X(0) and {d(n), v(n)} are independent Gauss’s vectors with statistics: E ½xð0Þ ¼ x^ð0j 0Þ ¼ xð0Þ;E ½dðnÞ ¼ E ½vðnÞ ¼ 0

 cov½xð0Þ ¼ E ½xð0Þ  xð0Þ½xð0Þ  xð0ÞT ¼ M ð0Þ

(6.24) (6.25)

cov½d ðnÞ ¼ R1 ðnÞ

(6.26)

cov½vðnÞ ¼ R2 ðnÞ

(6.27)

E[x(0)], E[d(n)], and E[v(n)] are mathematical expectations of the initial state, state, and measure disturbances, respectively. cov[x(0)], cov[d(n)], and cov[v(n)] are covariances of the initial state, state, and measure disturbances, respectively.  R1 ðnÞ ¼ diag R1t ðnÞ ON ðN 1ÞT a N ðN 1ÞT a R1t0 ðnÞ X  T t , t 2 T , t0 2 T a O2N ðN 1ÞðP1Þ t2T a

is a symmetric, positive-semidefinite matrix of covariances of the random component of demand for the original car types t and aggregate car     P types t0 with NT ð2N  1Þ + 2N ðN  1Þ T a + ðP  1Þ T + t2T a T t     P dimensions.  NT ð2N  1Þ + 2N ðN  1Þ T a + ðP  1Þ T + t2T a T t R2(n) ¼ diag [R2t(n)], t 2 T is a symmetric, positive-semidefinite matrix of covariances of white noise in the forecast with dimensions TN  TN.

6.5 Forecasting the state of rail freight cars by ARIMA-Kalman method The forecasts of the number of rail freight cars on state in each station are obtained by a sophisticated mathematical model based on a combination of the ARIMA forecasting technique and Kalman recursions. More precisely, based on the available historical data ARIMA processes for modeling the state of freight cars have been proposed. The estimated ARIMA model is incorporated into the state-space framework and classical Kalman recursions were applied for forecasting the future values.

186

Optimization models for rail car fleet management

The ARIMA method was introduced by Box and Jenkins (1976). The method represents one of the most frequently used univariate time series modeling tools. ARIMA models are based on the autoregressive (AR) model, moving average (MA), and the combination of AR and MA models (Milenkovi
c et al., 2016). The AR model involves lagged terms of the time series itself, MA model includes lagged terms on the noise or residuals. The first requirement for ARIMA modeling is that the time series data to be modeled is stationary or can be transformed into stationary. The letter “I” (integrated) means that the first-order difference is applied in order to stationarize given time series. The generalized form of an ARIMA(p, d, q) model for a series yl can be written as χ p ðωÞð1  ωÞd yl ¼ ψ q ðωÞεl

(6.28)

where εl is a white-noise sequence χ p ðωÞ ¼ 1  χ 1 ω  χ 2 ω2  ⋯  χ p ωp

(6.29)

is the AR polynomial term of order p ψ q ðωÞ ¼ 1  ψ 1 ω  ψ 2 ω2  ⋯  ψ q ωq

(6.30)

is the MA part of order q. (1  ω)d is the differencing operator of order d used to eliminate polynomial trends. ω is the backshift operator, whose effect on a time series yl can be summarized as ωdyl ¼ yld. The Box-Jenkins method includes three steps for fitting the ARIMA models. These are the identification, estimation, and validation of the model (Box et al., 2008). The first step is generally based on an analysis of the autocorrelation function (ACF) and partial autocorrelation function (PACF) and their comparison with theoretical profiles of these functions in AR, MA, and ARMA processes. The output of the identification step is determined as an appropriate (p, d, q) model structure. The model structure selected in the previous step has to be fitted to the time series and its parameters have to be estimated. This is the essence of the second step and it is done by using the conditional sum of squares or maximum likelihood method. Validation of the selected model is performed within the diagnostic check with the analysis of stationarity, invertibility as well as the presence of redundancy in model parameters. If a selected model fails in a diagnostic check, it is necessary to repeat the whole procedure again. After an appropriate model is found, it can be used for the forecasting purpose (Box et al., 2008).

Stochastic model for heterogeneous rail freight car fleet management

187

Subjectivity is mostly present in the model identification step. The reason for this lies in the fact that this step is essentially based on graphical interpretations of ACF/PACF estimates. To cope with this subjectivity and also to improve the determination the final orders of the ARMA processes, there are a lot of model selection criteria proposed in the literature (De Gooijer et al., 1985). The most frequently used are information criteria, such as Akaike information criterion (AIC), Bayesian information criterion (BIC), and the normalized version of the BIC. These information criteria consist of the natural log of the mean square error (MSE), plus a penalty for the number of parameters being estimated (Yaffee and McGee, 1999): AIC ¼ T ln ðMSEÞ + 2k

(6.31)

BIC ¼ T ln ðMSEÞ + k ln ðT Þ

(6.32)

lnðT Þ (6.33) T where T is the number of observations, k is the number of parameters in the model k ¼ p + q + P + Q + 1, and MSE is the mean squared error. The common procedure often involves estimating the largest model which is assumed to correctly capture the dynamics of a time series and then decreasing its size (dropping the lags) until the minimum value of AIC, BIC, or the normalized BIC is reached. The state-space form provides a unified representation of a wide range of linear Gaussian time series models including ARMA and unobserved component (UC) models (Hindrayanto et al., 2010). Statespace representations of ARMA/ARIMA models are presented in Box et al. (2008) and Brockwell and Davis (2002). The Gaussian state-space model for the Po-dimensional observation sequence y1, …, yPo is given as follows (Durbin and Koopman, 2001): Normalized BIC ¼ ln ðMSEÞ + k

ρl + 1 ¼ Yl ρl + I l ηl , ηl N ð0, M l Þ,l ¼ 1,…,Po

(6.34)

yl ¼ Sl ρl + εl , εl N ð0, V l Þ,ρ1 N ða1 , ΠÞ

(6.35)

where ρl is the state vector, ηl and εl are disturbance vectors, and the system matrices Y l , I l , Sl , M l , and V l are fixed and known but a selection of elements may depend on an unknown parameter vector. Eq. (6.34) is the state transition equation, while Eq. (6.35) represents the observation or measurement equation. The transition matrix Y l determines the dynamic evolution of the state vector. I l represents the state disturbance transform matrix. The

188

Optimization models for rail car fleet management

disturbance vector ηl for the state vector update has the zero mean and variance matrix M l . Matrix Sl links the observation vector yl with the unobservable state vector ρl. The observation disturbance vector εl has the zero mean and variance matrix V l . The observation and state disturbances, ηl and εl, are assumed to be serially independent and independent of each other in all time instances. The Kalman filtering (Kalman, 1960) allows a unified approach to prediction and estimation for all processes that can be given by state-space representation. Therefore, the state-space model can be treated by standard time series methods based on the Kalman filter (Durbin and Koopman, 2001; Anderson and Moore, 1979). The variance matrix Π of the initial state vector ρ1 may contain diffuse elements when nonstationary components are included in ρo. In this case, diffuse initialization methods for the Kalman filter exist to evaluate the exact or diffuse likelihood function (De Jong, 1991; Koopman, 1997). After the initialization of transition and measurement equations, the Kalman model can be used for forecasting the number of rail freight cars by applying the Kalman recursive prediction (KRP). Recursive steps of KRP are presented in Durbin and Kopman (2001).

6.5.1 Components of weighted matrices A, B, and L This section describes the procedure for determining the elements of weighted matrices as well as their proper structuring. The importance of this process for the quadratic optimality criterion and recommendations for their structuring were presented in Athans and Falb (1966). Comparing to the fuzzy rail freight car inventory model, in this case, the basic stochastic version is, however, extended to include each freight car type independently. Therefore, for components of matrices A, B, and L, the multiperiod EOQ model (Taha, 2003) was applied. Parameter γ it(n) represents the average freight car of type t (t 2 T) inventory level and shortage in station i for the period n of the planning horizon. Parameter ηit(n) represents the number of cars ordered for a station and a period created whenever the inventory of rail freight cars of type t drops to some reorder level. Determining the optimal number of freight cars and the optimal ordering quantity is based on a cost minimum consisted from a sum of fleet ownership and allocation costs. Dit represents the total demand for freight cars of type t in station i. Dit includes the demand specified in terms of the original car type t extended

Stochastic model for heterogeneous rail freight car fleet management

189

for the demand for the aggregate car type t0 (Dit0 ) which will be fulfilled by the original type t (t 2 T t0 , t0 2 Ta ). This extension is based on an estimation of the past share of type t freight cars in substitutions for satisfying the demand for the aggregate car type t0 . The setup cost can be expressed as N  Dit 1 X ϑt  ejit + ð1  ϑt Þljit , t 2 T ηit N  1 j¼1

(6.36)

Setup cost represents the cost of the delivery of freight cars of type t per period, from stations where there is an excess of freight cars to the station for which the optimal delivery has been determined. The relationship Dη it is the approximate number of orders per period. it The unit average cost of supplying an order of freight cars is PN  1 j¼1 ϑt  ejit + ð1  ϑt Þljit , t 2 T, where the parameter ϑt denotes the ratio N 1 of the number of empty cars of type t and the total number of freight cars. This coefficient is determined for all cars of the original type. Unit costs of empty and loaded movement, for each type t of freight cars between all stations j and station i, are represented by ejit and ljit, respectively. In order to determine the average fleet holding cost, we assumed that the number of cars of type t will vary linearly from ηit, the value at the beginning of the planning period, to γ it, the value at the end. Therefore, the average fleet holding cost can be expressed as !! N X γ it + ηit hit  ,t 2 T (6.37) Dijt E 2 j¼1 Parameter hit is the unit holding cost for a freight car of type t in a station i. Dijt is the average daily demand between stations i and j per period defined on P  1 periods of the planning horizon. In order to formulate the cost of shortage for a rail freight car of type t, we introduced the function: δpt

XN j¼1

Dijt , γ i



8 X < 0, D < γ it j ijt X ¼ X , : D  γ , D  γ ijt ijt it it j j

t2T

(6.38)

190

Optimization models for rail car fleet management

The expected total demand for cars of type t, backordered or transferred to a future period is

X X ð∞ X δpt ¼ Dijt  γ it f Dijt d Dijt γ it j j j ð∞ X X Dijt f Dijt  γ it Y ðγ it Þ,t 2 T (6.39) ¼ γ it

j

j

Y(γ it) is a complementary cumulative distribution function. The demand for cars of type t between stations on the network comply with the normal distribution law ( Jordan, 1982). Considering the central limit theorem, the sum of two Gaussian random variables also follows the normal distribution law. Let this sum be denoted by κit: κit ¼ Di1t + Di2t + ⋯ + Dijt + ⋯ + DiNt , t 2 T

(6.40)

Considering that Dijt( j 2 N, t 2 T ), and with N(μijt, σ ijt) normally distributed, there is κ it with N(μit, σ it), where μit ¼ μi1t + μi2t + ⋯ + μiNt

(6.41)

σ 2it ¼ σ 2i1t + σ 2i2t + ⋯ + σ 2iNt

(6.42)

The total expected cost for a freight car of type t for a station i over the planning period can be formulated as N  N X Dit 1 X η + γ it E½C ðγ it , ηit Þ ¼ ϑt ejit + ð1  ϑt Þljit + hit it Dijt E ηit N  1 j¼1 2 j¼1 N Tit X + pijt ηit j¼1

"ð

N ∞X

X

γ i j¼1

j

Dijt f

!

!!

#

Dijt  γ it Y ðγ it Þ , t 2 T (6.43)

in order to find γ it and ηit which minimize the function of total estimated inventory costs for each freight car type t. The set of necessary conditions is determined based on the following two conditions: N ∂E ½C ðγ it , ηit Þ Dit 1 X hit ¼ 2 ϑt ejit + ð1  ϑt Þljit Þ + ηit N  1 j¼1 2 ∂ηit 



 N Dit X γ it  μit γ it  μit  ðμit  γ it ÞΦ ¼0 pijt σ it ϕ  2 ηit j¼1 σ it σ it

(6.44)

Stochastic model for heterogeneous rail freight car fleet management



N ∂E ½C ðγ it , ηit Þ hit Dit X γ it  μit ¼0 ¼  pijt Φ 2 ηit j¼1 σ it ∂γ it where

191

(6.45)

1 2 ϕðγ it Þ ¼ pffiffiffiffiffi eγ =2 2π

and

∞ ð

Φðγ Þ ¼

ϕðηÞdη γ

is the complementary cumulative distribution of the demand density function. Now, for a given Dit, Dijt, pijt, eijt, lijt, σ it, μit, the optimal rail freight car inventory level and the optimal reordering number of freight cars of type t in a station i during the period n can be calculated by the following algorithm: Step 0. Find the initial, smallest value of ηit which is η∗it, and let γ (0) it ¼ 0. Set k ¼ 1 and go to stepk. ffi vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi XN  u u2Dit ϑ e + ð1  ϑt Þljit t j¼1 t jit ð1Þ ηit ¼ η∗it ¼ (6.46) ðN  1Þhit (k) (k) (k1) Step k. On the base of η(k) it calculate γ it from Eq. (6.45). If ηit ηit (k) (k) ∗ stop, the optimal solution is γ ∗it ¼ γ (k) it and ηi ¼ ηi . Otherwise, use ηit in (k) Eq. (6.46) to compute γ it . Set k ¼ k + 1, and repeat step k until convergence. Since the point (γ ∗it, η∗it) represents the global minimum of the total cost function, it can be approximated by Taylor’s polynomial by an expansion to the second term:

C ðγ it , ηit Þ C ðγ it ∗ , ηit ∗ Þ + +

∂C ∂C 1 ∂2 C ðγ it  γ it ∗ Þ + ðηit  ηit ∗ Þ + ðγ  γ it ∗ Þ2 ∂γ t ∗ ∂ηt ∗ 2 ∂γ t ∗ it

∂2 C 1 ∂2 C ðγ it  γ it ∗ Þðηit  ηit ∗ Þ + ðη  ηit ∗ Þ2 ∂γ t ∂ηt 2 ∂ηt ∗ it (6.47)

The coefficients of the quadratic terms γ it and ηit can be determined based on the cost function of the following form (Monks, 1982): 1 1 Ct ¼ π st γ 2 + π ht γη + π ut η2 2 2

(6.48)

192

Optimization models for rail car fleet management

Coefficients π s, π h, and π u will be derived from expression (6.47). By analyzing the quadratic terms, the unknown coefficients for all car types are ∗

N Dit X γ it  μit s (6.49) π it ¼ ∗ pijt ϕ σ it ηit σ it j¼1 π hit

¼

∗

γ it  μit pijt Φ σ it j¼1

N Dit X

η∗it

2

(6.50)

N  N 1 X Dit X + 2 ϑ e + ð 1  ϑ Þl pijt t jit t jit 3 3 η∗it N  1 j¼1 η∗it j¼1  ∗

∗

 γ  μit γ  μit  σ it ϕ it  ðμit  γ ∗it ÞΦ it σ it σ it

π uit ¼ 2

Dit

(6.51)

The formulated procedure for finding the unknown coefficients π st, π ht , and π ut has to be applied for each station and each original freight car type. What remains to be done, is to define the quadratic cost coefficients for variables that describe the freight car flows. Namely, the state vector, beside the number of freight cars by types, undispatched car flows, and demand requests, also contains the variables that symbolize the cars in transit during the particular planning period. The quadratic cost coefficient for these variables is assumed to be the linear function of the number of cars on the route considered: ! P X xjitt0 m ðnÞQ ¼ xjitt0 m ðnÞ qt0 + qt xjitt0 m ðnÞ sajits ðn + s  mÞ , t 2 T, t 0 2 T t s¼m

yjitt0 m ðnÞQ ¼ yjitt0 m ðnÞ qt0 + qt yjitt0 m ðnÞ

P X

!

(6.52)

sbjits ðn + s  mÞ , t 2 T, t 0 2 T t

s¼m

(6.53) It is now possible to determine the coefficients in matrices A(n), B, and L of the generalized quadratic performance index.

Stochastic model for heterogeneous rail freight car fleet management

193

Therefore, matrix A(n) has the structure as follows: 2 3 A1 ðnÞ O T ðN ð2N  1Þ + 2N ðN  1ÞðP  1ÞÞ X 2N ðN  1ÞðT a + ðP  1Þ Tt Þ 7 6 6 7 t2T AðnÞ ¼ 6 7 4 O ð2N ðN  1ÞðT a + ðP  1Þ X Tt Þ 5 A2 ðnÞ a

t2T a

T ðN ð2N  1Þ + 2N ðN  1ÞðP  1ÞÞ

(6.54) A1 ðnÞ ¼ diagðAt ðnÞÞ, t 2 T At ðnÞ ¼ diagðπ s1t ,…,π s1t ,π s2t ,…,π s2t , π sNt , …, π sNt ,

π s1t ,…,π s1t ,…,π s2t ,…,π s2t ,…,π sNt ,…,π sNt ,2qt P P X X  sa21ts ðn + s  1Þ,…,2qt saN 1ts ðn + s  1Þ,…, 2qt s¼1 P X



s¼1

sa21ts ðn + s  P + 1Þ,…,2qt

s¼P1 P X

saN 1ts ðn + s  P + 1Þ,…,2qt

P X

sb21ts ðn + s  1Þ,…,2qt s¼P1 s¼1 P X sbN 1ts ðn + s  1Þ,…,2qt s¼1 P P X X  sb21ts ðn + s  P + 1Þ,…,2qt sbN 1ts ðn + s  P + 1Þ,…,2qt s¼P1 s¼P1 P P X X  sa1Nts ðn + s  1Þ,…,2qt sbN 1, Nts ðn + s  P + 1ÞÞ, t 2 T s¼1 s¼P1   A2 ðnÞ ¼ diag At 0 ðnÞ , t0 2 T a

At 0 ðnÞ ¼ diagð2qt



s¼1

sa21ts ðn + s  P + 1Þ, …, 2qt Þ

s¼P1 P X

saN 1ts ðn + s  P + 1Þ, …, 2qt

s¼P1 P X

P X sb21ts ðn + s  1Þ,…,2qt s¼1

(6.57)

sbN 1ts ðn + s  1Þ,…,2qt

 

s¼1 P X

sb21ts ðn + s  P + 1Þ, …, 2qt

s¼P1 P X

sa1Nts ðn + s  1Þ…2qt

s¼1

(6.56)

P P X X sa21ts ðn + s  1Þ,…, 2qt saN 1ts ðn + s  1Þ, …,2qt s¼1

P X

(6.55)

P X

P X

sbN 1ts ðn + s  P + 1Þ, …,2qt

s¼P1

sbN 1, Nts ðn + s  P + 1ÞÞ, t 0 2 T a , t 2 T t’

s¼P1

194

Optimization models for rail car fleet management

matrix A(n) is a ½NT ð2N  1Þ + 2N ðN  1ÞðT a + ðP  1Þ    P T + t2T a T t Þ T + t2T a T t Þ  ½NT ð2N  1Þ + 2N ðN  1ÞðT a + ðP  1Þ P diagonal square matrix and has the components π sit, qt Ps¼m saijts ðn + s  mÞ, PP and qt s¼m sbijts ðn + s  mÞ. B is a ½2NT ðN  1Þ + h i P P 2N ðN  1Þ t0 2T a T t0   2NT ðN  1Þ + 2N ðN  1Þ t’2T a T t’ diagonal square Thus, 

P

the

matrix which is a positive-definite matrix with components π uit and π uit’, t 2 T, t 0 2 Ta. B ¼ diag(B1, …, Bt, …, BT) where

  Bt ¼ diag π u1t , …, π u1t , π u2t , …, π u2t , …, π uNt , …, π uNt , π u1t , …, π u1t , π u2t , …, π u2t , …, π uNt , …, π uNt (6.58)

P a Matrix L is a ½NT ð2N  1Þ + 2N  PðN  1ÞðT + ðP  1ÞðT + t2T a T t ÞÞ h T ð2N ðN  1ÞÞ + 2N ðN  1Þ t0 2T a T t0 matrix with components π i and rijtt0 2

L1

X t 3 T 7 t2T a 7 7 L2 X X t 7 5 O 2N ðN 1ÞðP1Þ Tt  2N ðN  1Þ T

O

6 6 L¼6 ðN 1Þ 6O2N ðN1ÞTa 2NT X 4O 2N ðN 1ÞðP1Þ Tt  2NT ðN  1Þ t2T a

T ðN ð2N 1Þ + 2N ðN 1ÞðP1ÞÞ2N ðN1Þ

t2T a

L1 ¼ diag (Lt), t 2 T where 2 h M1 6 6 O1 6 6⋯ 6 6 6 O1 6 h Lt ¼ 6 6 N1 6 6 O2 6 6⋯ 6 6 4 O2

t2T a

(6.59) O1 O1 M1h O1 O1

3

7 M2h O1 O1 M2h O1 7 7 ⋯ ⋯ ⋯ ⋯ ⋯7 7 h h 7 O1 MN O1 O1 MN 7 7 O2 O2 N1h O2 O2 7 7 7 N2h O2 O2 N2h O2 7 7 ⋯ ⋯ ⋯ ⋯ ⋯7 7 7 O2 NNh O2 O2 NNh 5

O3 O3 O3 O3 O3 O3 Mh1 ¼ [π h1t]NN1 Nh1 ¼ [π h1t]N1N1 Mh2 ¼ [π h2t]NN1 Nh2 ¼ [π h2t]N1N1 MhN ¼ [π hNt]NN1 NhN ¼ [π hNt]N1N1 O1 is a N  N  1 zero matrix O2 is a N  1  N  1 zero matrix O3 is a [2N(N  1)(P  1)]  [(N  1)] zero matrix

(6.60)

6 6 6 6 6 6 6 L2 ¼ 6 6 6 6 6 6 4

diagðrij1t10 Þ

⋯

0

  diag rijT t10 t10

diagðrijt1 t10 Þ

⋯

diagðrij1t20 Þ

⋯

diagðrijt2 t20 Þ

⋯

⋮

⋮

⋮

⋮

0

0

diagðrij1T a Þ

⋯

⋮  diag rijt3 T a

⋯

⋮   diag rijT T a T a

⋮

⋮

⋮

⋮

⋮

⋮

⋮

0

0

0

0

0

0

0

i, j 2 N , t1 ¼ 1,…, T t10 ,t2 ¼ 1,…, T t20 , t3 ¼ 1, …,T T a , t10 , t20 2 T a

0

0

  diag rijT t20 t20

0

⋮

t10

3 7 7 7 7 7 7 7 7 7 7 7 7 7 5

(6.61)

Stochastic model for heterogeneous rail freight car fleet management

2

195

196

Optimization models for rail car fleet management

6.5.2 The components of matrix Γ(P) The number of cars in the last period of the planning horizon influences the sensitivity of decisions about their routing within the given planning horizon. If these decisions are not of critical importance for the supply of freight cars by stations, it is not necessary to take this value into consideration. In the case of models proposed in this book, this results from the assumed cyclic nature of demand on the period basis. For determining the components for the last term, Γ(P) a random newsboy model is applied. It includes a single variable (system state) and a single time period and it is presented by a second-order gradient method through quadratic and linear terms in state variables. The last period of the planning horizon lasts one day. This model determines the optimal level of freight cars for a single period in a particular station in case of stochastic demand and when the objective is to minimize the expected total cost. Based on a single period newsboy inventory model for each freight car type, it is possible to determine the coefficients of the weighted matrix of the P terminal P member. Therefore, in case if N D j¼1 ijt < γ it holds, then the quanN tity γ  D is held during the period. Otherwise, a shortage amount j¼1 ijt PN it PN j¼1 Dijt  γ it will result if j¼1 Dijt > γ it .The expected cost for the period, E[C(γ it)], is expressed as ! ð γit N X γ it  Dijt f ðDit Þd ðDit Þ E ½C ðγ it Þ ¼ hit 0 j¼1 ! (6.62) ð∞ X N N X Dijt  γ it f ðDit ÞdðDit Þ + pijt j¼1

γi

j¼1

P where N j¼1Dijt N(μit, σ it). A necessary condition for determining the optimal number of cars in a station is ∂E½C ðγ it Þ ¼0 ∂γ it ∗

γ it  μit hit  ¼ Φ XN σ it hit + p j¼1 ijt

(6.63) (6.64)

where Φ(k) is the complementary cumulative distribution function. Considering that the function of total costs is strictly concave, its absolute minimum is unique. The second derivative of this function is

Stochastic model for heterogeneous rail freight car fleet management N X ∂2 E½C ðγ it Þ ¼ h + pijt it ∂γ 2it j¼1

!

 γ it  μit ϕ >0 σ it

197

(6.65)

Solving Eq. (6.64) will result in the optimal rail freight car inventory level γ ∗it for the last day of the planning horizon. In the vicinity of the minimum this function can be approximated by the Taylor’s polynomial by expansion to the second term: C ðγ it Þ C ðγ ∗it Þ +

∂C 1 ∂2 C | ðγ it  γ ∗it Þ + ðγ it  γ ∗it Þ 2 | ðγ t  γ ∗t Þ ∂γ it γ ∗ 2 ∂γ t γ ∗t

(6.66)

it

The coefficient of the quadratic term is   γ ∗  μ XN it ,i 2 N ,t 2 T p ϕ it git ¼ hit + j¼1 ijt σ it

(6.67)

    P Γ(P) is a NT ð2N  1Þ + 2N ðN  1Þ T a + ðP  1Þ T + t2T a T t      P NT ð2N  1Þ + 2N ðN  1Þ T a + ðP  1Þ T + t2T a T t positivesemidefinite matrix in the quadratic term of the state vector for the last period of the planning horizon. 2 6 6 Γ ðP Þ ¼ 6 4O

Γ 1 ðP Þ

O T ðN ð2N  1Þ + 2N ðN  1ÞðP  1ÞÞ 2N ðN  1ÞðT a + ðP  1Þ

ð2N ðN  1ÞðT a + ðP  1Þ

X

Γ2 ðP Þ

Tt Þ

X

t2T a

Tt Þ

3 7 7 7 5

(6.68)

t2T a

T ðN ð2N  1Þ + 2N ðN  1ÞðP  1ÞÞ

Γ 1(P) ¼ diag(Γ t(P)), t 2 T where Γ t ðP Þ ¼ diagðg1t , …, g1t , …, gNt ,…, gNt , g1t , g2t , …,gNt ,2qt

P X

sa21ts ðP + s  1Þ,  , 2qt

s¼1



P P P X X X saN 1ts ðP + s  1Þ,2qt sa21ts ðs + 1Þ,  , 2qt saN 1ts ðs + 1Þ, 2qt s¼1

s¼1

s¼P1

P P P X X X  sb21ts ðP + s  1Þ,  , 2qt sbN 1ts ðP + s  1Þ, 2qt sbN 1ts ðs + 1Þ,  ,2qt s¼1

s¼1

s¼1

s¼P1

s¼P1

P P X X  sa1Nts ðP + s  1Þ,2qt sbN 1, Nts ðs + 1ÞÞ, t 2 T

(6.69) 0

Γ 2 ðP Þ ¼ diagðΓ t0 ðP ÞÞ, t 2 T

a

(6.70)

198

Optimization models for rail car fleet management P P X X sa21ts ðP + s  1Þ,  ,2qt saN 1ts ðP + s  1Þ,2qt sa21ts ðs + 1Þ,2qt s¼1 s¼1 s¼1 P P P X X X  saN 1ts ðs + 1Þ,2qt  sb21ts ðP + s  1Þ,  ,2qt sbN 1ts ðP + s  1Þ,2qt s¼P1 s¼1 s¼1 P P P X X X  sbN 1ts ðs + 1Þ,  ,2qt sa1Nts ðP + s  1Þ,2qt sbN 1, Nts ðs + 1ÞÞ,t 0 2 T a ,t 2 T t0 s¼P1 s¼1 s¼P1 (6.71)

Γ t 0 ðP Þ ¼ diagð2qt

P X

6.6 MPC controller MPC is a multivariable control method usually applied for the optimal control of complex systems with an explicit consideration of constraints. The dynamics of the system is represented by a discrete time predictive model by which the state of the system in the next period is a function of the actual state, actual demands (disturbance), and the actual control vector with constraints. In each period, the optimal control problem is solved online based on the actual state (at time n ¼ n0) and the predictive demands over a N step finite horizon (at time n ¼ n0, …, n0 + N  1). The result of the optimization is a sequence of control vectors overtime n ¼ n0, …, n0 + N  1 but only the first control vector (at time n ¼ n0) is applied to the system. In the next state (at time n ¼ n0 + 1), the optimal control problem is solved again for the time horizon n ¼ n0, …, n0 + N, then only the control vector at time n ¼ n0 + 1 is applied and so forth.

6.6.1 Optimization problem The stochastic rail freight car fleet sizing and allocation model is represented by a discrete linear model in which the state X(n), control U(n), and disturbance are related in accordance to the following relation: X ðn + 1Þ ¼ Λ½n, aðnÞ, bðnÞX ðnÞ + GU ðnÞ + ΦdðnÞ, n ¼ 0,1,…,P  1 (6.72) with X(0) given and subject to the following system constraints: U ðnÞ  O2

(6.73)

X ðnÞ  O3

(6.74)

DX ðnÞ  K

(6.75)

Based on the system dynamics (6.72) and constraints (6.73)–(6.75), the problem is to determine the set of control actions U(n) in order to minimize the objective functional of the form:

199

Stochastic model for heterogeneous rail freight car fleet management

ð6:76Þ In order to apply the MPC approach, the objective functional has been transformed into the standard form (Anderson and Moore, 2007): P1   1 1X X ðnÞT Q1 ðnÞX ðnÞ + Ue ðnÞT BU e ðnÞ J ¼ X ðP ÞT ΓðP ÞX ðP Þ + 2 2 n¼0

(6.77) where Q1(n) ¼ A(n)  LB1LT and Ue(n) ¼ U(n) + B1LTX(n). The system dynamics constraint (6.72) now takes the form X ðn + 1Þ ¼ N ðnÞX ðnÞ + GU e ðnÞ + ΦdðnÞ, n ¼ 0,1,…, P  1

(6.78)

where N(n) ¼ Λ[n, a(n), b(n)]  GB1LT. Consequently, the problem of the minimization of the criterion (6.77) subject to dynamics (6.78) and the set of constraints (6.73)–(6.75) can be solved as an MPC problem.

6.6.2 Detailed description of the MPC approach In this part, the MPC approach is clarified. The proposed controller is parameterized by a set of discrete decision periods 0,1, …, m, …, n, …, P, a positive-definite matrix Ψ of [T(2N(N  1) + 2N(N  1)(P  1))]  [T(2N(N  1) + 2N(N  1)(P  1))] dimensions and a vector Ω of dimension [T(2N(N  1) + 2N(N  1)(P  1))]. At any time n, the controller applies the optimal solution of a quadratic programming (QP) problem in which the decision variables represent the empty and loaded flow control actions over the time horizon n, n + 1, …, n + P  1 for each type of freight car. The following state-space model represents the considered rail freight car fleet system: X ðn + 1Þ ¼ N ðnÞX ðnÞ + GU e ðnÞ + Φd ðnÞ, n ¼ 0,1,…, P  1 Z ðnÞ ¼ HX ðnÞ + vðnÞ, n ¼ 0,1,…, P

(6.79) (6.80)

The state noise d(n) and measurement noise v(n) are assumed to be Gaussian distributed with zero mean and respective covariances of R1(n) and R2(n) with cross-covariance R3(n):       R1 ðnÞ R3 ðnÞ dðnÞ 0 (6.81) N , R3T ðnÞ R2 ðnÞ vðnÞ 0

200

Optimization models for rail car fleet management

The defined model is used for predictions about the rail freight car system in the following periods of the prediction horizon, using information (measurements of inputs and outputs) up to and including the current period n (Wills, 2004). Now it is possible to make optimal (from the aspect of minimum variance) predictions of state and output by applying the Kalman filter (Athans and Falb, 1966; Kalman, 1960). Accordingly, let X ðij jÞ represents an estimate of the state at time i given information up to and including time j, where j  i. Then, X ðn + 1j nÞ ¼ N ðnÞX ðnj n  1Þ + GU e ðnÞ + K ðnÞ½Z ðnÞ  Z ðnj n  1Þ (6.82) Z ðnj n  1Þ ¼ HX ðnj n  1Þ  1 K ðnÞ ¼ N ðnÞP ðnÞH T + R3 ðnÞ HP ðnÞH T + R2 ðnÞ 

P ðnÞ ¼ R1 ðnÞ + N ðnÞP ðnÞΛT ðnÞ  ðN ðnÞP ðnÞH T + R3 ðnÞÞ 1  ðHP ðnÞH T + R2 ðnÞÞ R3T ðnÞ + HP ðnÞN T ðnÞ

(6.83) (6.84) (6.85)

The Kalman filter is represented by recurrent relations (6.82) and (6.83). In fact, these two relations are the model of the system extended for a correction term, which is proportional to the difference between the values of Z(n) and the values estimated on the basis of the fuzzy state function H  X ðnÞ. Matrix (6.84), K(n), is the Kalman filter gain. The discrete-time algebraic Riccati equation (DARE) for P(n) is given by Eq. (6.85). To apply the MPC approach, the estimates of the state over the whole prediction horizon need to be determined. These predictions can be made only based on the information up to and including the current time n. Eqs. (6.82)–(6.85) can be used to obtain X ðn + 1j nÞ, and optimal estimates over the prediction horizon can be obtained as follows: X ðn + i + 1j nÞ ¼ N ðnÞX ðn + ij nÞ + GU e ðn + ij nÞ

(6.86)

Z ðn + ij nÞ ¼ HX ðn + ij nÞ

(6.87)

Here, i ¼ n, n + 1, …, n + P  1 and notation Ue(n + ij n) is applied to distinguish the actual input at time t + i from the input used for prediction purposes, namely Ue(n + i j n). Based on the estimated current state X ðnÞ and the controls Ue(i), i ¼ n, n + 1, …, n + P  1, the state predictions are generated and used for the

Stochastic model for heterogeneous rail freight car fleet management

201

formulation of the objective and constraints of the QP problem (Le et al., 2013). The detailed description of the QP formulation is given below. The main components of the prediction, optimization, and receding horizon implementation constitute the MPC law. The dynamic model (6.79) is used for the prediction of the future response of the controlled system. For the given predicted input sequence, the corresponding sequence is generated by simulating the model forward over the prediction horizon of P  1 periods. These predicted sequences are places in vectors U, X defined by 2 3 2 3 Ue ðnj nÞ X ðn + 1j nÞ 6 Ue ðn + 1j nÞ 7 6 7 7,X ðnÞ ¼ 6 X ðn + 2j nÞ 7 Ue ðnÞ ¼ 6 (6.88) 4 5 4 5 ⋮ ⋮ Ue ðn + P  1j nÞ X ðn + Pj nÞ Ue(n + ij n) and X ðn + ij nÞ are the input and state vectors at time n + i which are predicted at time n. Then, X ðn + ij nÞ evolves according to the prediction model: X ðn + i + 1Þ ¼ N ðnÞX ðn + ij nÞ + GU e ðn + ij nÞ + Φd ðnÞ, i ¼ 0,1, …,P  1 (6.89) The initial condition at the beginning of the prediction horizon is defined as X ðnj nÞ ¼ X ðnÞ In compact form, the predicted state sequence generated by Eq. (6.89) looks as follows: X ðn + ij nÞ ¼

0 Y

N ðn + jÞX ðnÞ + Ci ðnÞUe ðnÞ + Ei ðnÞdðnÞ,i ¼ 0,…,P  1

j¼i1

(6.90) or X ðnÞ ¼ N X ðnÞ + CU e ðnÞ + Ed ðnÞ where 2

3 N ðnÞ 6 7 N ðn + 1ÞN ðnÞ 7 N ¼6 4 5 ⋮ N ðP  1ÞN ðP  2Þ⋯N ðnÞ

(6.91)

202

Optimization models for rail car fleet management

C is the (convolution) matrix with rows Ci defined by 2 3 G 0 ⋯ 0 6 7 6 7 N ð n + 1 ÞG G ⋯ 0 6 7 C¼6 7 6 7 ⋮ ⋮ ⋮ ⋮ 4 5 N ðP  1Þ⋯N ðnÞG N ðP  2Þ⋯N ðnÞG ⋯ G C0 ¼ 0, Ci ¼ ith block row of C 2 Φ 6 N ðn + 1ÞG 6 E¼6 4 ⋮

0 Φ ⋮

3 ⋯ 0 ⋯ 07 7 7 ⋮ ⋮5

(6.92)

(6.93)

N ðP  1Þ⋯N ðnÞΦ N ðP  2Þ⋯N ðnÞΦ ⋯ Φ E0 ¼ 0, Ei ¼ ith block row of E We obtain the performance criteria as a quadratic function of Ue(n) by a substitution for X ðn + ij nÞ in Eq. (4-77): 1 J ðnÞ ¼ Ue ðnÞT Ψ Ue ðnÞ + ΩUe ðnÞ + Υ 2 e1C + B e Ψ ¼ CT Q

(6.94) (6.95)

e 1 C + d ðnÞT E T Q e1C Ω ¼ X ðnÞT N T Q (6.96)  1 e X ðnÞ + 1 X ðnÞT N T Q e1N + Γ e 1 EdðnÞ Υ ¼ X ðnÞT N T Q 2 2 (6.97) 1 1 T T e T T e + dðnÞ E Q1 N X ðnÞ + dðnÞ E Q1 EdðnÞ 2 2 Matrix Ψ is a constant positive-definite matrix, Ω is a vector of corresponding dimensions, and Υ is a scalar which depends on X ðnÞ and d(n) with 2

3 2 3 2 3 Q1 0 … 0 B 0 … 0 Γ 0 … 0 6 0 ⋱ 6 ⋮ 7 0 ⋱ ⋮7 ⋮7 7 e 6 7 ,and Γ 7 e¼60 ⋱ e1 ¼ 6 Q 6 7,B ¼ 6 4 5 4 ⋮ B 0 ⋮ Γ 05 4 ⋮ Q1 0 5 0 ⋯ 0 B 0 ⋯ 0 Γ 0 ⋯ 0 Q1 Quadratic cost is now defined as a function of input predictions. The next step is to transform the set of constraints (6.73)–(6.75) in the form AUe  B which can be included in solving the problem of optimizing Eq. (6.94).

Stochastic model for heterogeneous rail freight car fleet management

203

Input constraints (6.73) can be expressed as Ue(n)  O2. Applied to predictions, Ue(n + ij n), i ¼ 0, …, P  1 these constraints can be expressed in terms of Ue(n) as IU e ðnÞ  O2

(6.98)

where I represents a unit matrix. In the same way, the state constraints (6.74) applied to predictions X ðn + ij nÞ, i ¼ 1, …,P are equivalent to Ci Ue ðnÞ  O3 + Ni X ðnÞ + Ei dðnÞ, i ¼ 0,…, P  1

(6.99)

Capacity constraints (4-75) are formulated as follows: DC i Ue ðnÞ  DRi  DSi + K, i ¼ 0,…,P  1 where

2 6 6 R¼6 4

N ðnÞ

(6.100)

3

7 N ðn + 1ÞN ðnÞ 7 7X ðnÞ ¼ N X ðnÞ 5 ⋮ N ðP  1ÞN ðP  2Þ⋯N ðnÞ

Ri is the ith block row of matrix R 2 ΦðnÞdðnÞ 6 N ðn + 1ÞΦðnÞdðnÞ + Φðn + 1Þd ðn + 1Þ 6 S¼6 4 ⋮

(6.101)

3 7 7 7 5

N ðP  1Þ⋯N ðn + 1ÞΦðnÞdðnÞ + ⋯ + Φðn + P  1Þd ðn + P  1Þ (6.102) Si is the ith block row of matrix S, n ¼ 0, …, P  1. In summary, the MPC formulation for the stochastic rail freight car fleet sizing and allocation problem can be stated as the minimization of a quadratic objective over Ue(n) subject to a set of linear constraints:

where

1 minimize Ue ðnÞT Ψ Ue ðnÞ + ΩUe ðnÞ 2 Ue ðnÞ

(6.103)

subject to Ae ðnÞUe ðnÞ  Be ðnÞ

(6.104)

2 3 3 I O2 Ae ðnÞ ¼ 4 Ci 5, Be ðnÞ ¼ 4 O3 + Ni X ðnÞ + Ei dðnÞ 5, i ¼ 0,…, P  1 DRi  DSi + K DC i (6.105) 2

204

Optimization models for rail car fleet management

The stated problem is a QP problem. Ψ matrix is a positive-definite matrix and the constraints are linear. Therefore, this problem is a convex problem which is solved by the QP solver (quadprog) in the Matlab (MathWorks Inc., 2012) for each period n. As it has been stated, only the first element of the optimal predicted input sequence is selected as an input to the rail freight car fleet system and therefore, for updating the state X(n) and the right side of the inequalities, Be(n) in each time period within the planning horizon.

6.6.3 Numerical experiments In this section, a set of numerical experiments has been carried out in order to demonstrate the validity of the developed MPC approach. The approach is explained through an example based on a network of four stations (N ¼ 4). The results of other test instances are presented at the end of the section. All test cases are based on the real rail freight car operation practice of Serbian railways. Four original car types (T ¼ 1, 2, 3, 4) are used for freight transport on the considered rail network. These types also satisfy the demand for two aggregate car types (Ta ¼ 5, 6). This assumption is based on the existence of four basic types (open, closed, flat, and other) on most railways and the common interchangeability between open and flat as well as closed and other car types. Fig. 6.2 illustrates possible substitution possibilities. As it can be noticed, the demand for the aggregate car type 5 will be fulfilled by either an original car of type 1 (open) or 3 (flat), while the demand for cars of type 1 or 3 can be fulfilled only by a corresponding type. The same situation holds for the aggregate car type 6 and original car types 2 (closed) and 4 (other). The planning horizon is composed of P ¼ 4 days where each day is one decision period. The length of 4 days is considered having in mind the longest rail freight car traveling time on Serbian railways.

Fig. 6.2 Substitution possibilities between car types.

Stochastic model for heterogeneous rail freight car fleet management

205

Table 6.2 contains the input data on unit costs of empty trips (ejit), unit costs of loaded trips (ljit), and unit car shortage costs (pijt) in monetary units per car for all routes and rail freight car types. Unit costs of empty and loaded trips are estimated and expressed in national currency. Expected travel times and variances for loaded and empty freight cars are given in Tables 6.3 and 6.4. Proportions of arrivals for loaded rail freight car flows together with the quadratic term in the unit ownership cost, as well as the proportions of arrivals for empty freight car flows of different types are given in Tables 6.5 and 6.6. The coefficient of empty trips is 0.35 for all freight car types. (See Table 6.2.) The demand is assumed to have a stochastic nature. Table 6.7 contains the mean value and standard deviation of daily demand for all periods, types of freight cars, and all o-d pairs of nodes. Unit holding costs (hit) per car per period, the estimate of initial number of cars (Sit(0)) per type per station, available capacities for the storage of freight cars by types (CAPi), and the lengths of rail cars are presented in Table 6.8. Table 6.9 contains unmet freight car requests from the preceding period. The numbers of loaded and empty cars that were initially in transit at the beginning of the horizon are given in Tables 6.10 and 6.11, respectively. Table 6.12 contains the cost of substitution which represents a penalty for using one of the original car types for satisfying the demand for aggregate car types 5 and 6. Random inventory models have been used for filling in the matrices A(0), A(1), A(2), A(3), B, and L. Optimal values for the freight car inventory level (γ it) and the ordering quantity (ηit) for each type have been determined after three to five iterations by using the approach described in Section 6.5.1 (Table 6.13). The random newsboy inventory model (Taha, 2003) has been used for determining the values of the car inventory level (γ Pit ) by stations in the last period of the planning horizon. Based on these values, the matrix Γ(3) has been defined (Table 6.12). For the estimation of the state vector, the ARIMA-Kalman method is applied to forecast the number of rail freight cars in each station. The method is applied for each car type t and station i. Forecasting results are given for the first station and first rail freight car type. A data sample of 100 days before the beginning of the planning horizon was used. The time series representing the state of rail freight cars in the station during the previous 100 days is given in Fig. 6.3. We used the first 85 days to fit the ARIMA models and the last 15 days as a holdout period to evaluate the forecasting performance.

206

Freight car type I

II

III

IV

Route

ejit

ljit

pijt

ejit

ljit

pijt

ejit

ljit

pijt

ejit

ljit

pijt

A-B A-C A-D B-A B-C B-D C-A C-B C-D D-A D-B D-C

300 200 200 100 125 150 150 200 100 200 150 125

100 70 70 85 100 120 100 120 70 120 100 85

300 200 200 200 250 300 300 400 200 400 300 250

300 200 200 100 70 200 100 100 150 150 150 150

100 70 70 70 70 100 70 70 100 100 100 100

300 200 200 200 200 300 200 200 300 300 300 300

200 200 300 200 100 125 150 200 100 150 125 200

70 70 100 120 70 85 100 120 70 100 85 120

200 200 300 400 200 250 300 400 200 300 250 400

200 200 300 150 150 125 200 200 150 200 200 300

70 70 100 100 100 85 120 120 100 140 140 200

200 200 300 300 300 250 400 400 300 400 400 500

Optimization models for rail car fleet management

Table 6.2 Cost parameters for rail freight car fleet sizing and the allocation problem

Table 6.3 Expected travel times and variances for loaded freight cars II

III

IV

Route

Expected travel time (days)

Travel time variance (days2)

Expected travel time (days)

Travel time variance (days2)

Expected travel time (days)

Travel time variance (days2)

Expected travel time (days)

Travel time variance (days2)

A-B A-C A-D B-A B-C B-D C-A C-B C-D D-A D-B D-C

1.3 1.2 2.0 1.0 1.5 1.3 1.5 1.5 1.0 2.0 2.0 1.0

0.30 0.40 0.50 0.50 0.30 0.30 0.30 0.30 0.50 0.50 0.50 0.60

1.2 1.0 1.0 1.4 1.5 1.5 1.0 1.5 1.0 1.0 2.0 1.5

0.50 0.40 0.50 0.40 0.50 0.50 0.40 0.50 0.40 0.40 0.60 0.50

1.0 1.5 2.5 1.0 2.0 1.0 1.0 2.0 1.4 2.5 1.0 1.0

0.55 0.50 0.50 0.55 0.65 0.55 0.55 0.65 0.45 0.50 0.55 0.55

1.0 1.5 1.5 1.0 2.0 1.5 1.0 2.0 1.0 2.0 1.5 1.0

0.50 0.50 0.50 0.50 0.65 0.50 0.50 0.65 0.50 0.65 0.50 0.50

Stochastic model for heterogeneous rail freight car fleet management

I

207

208

I

II

III

IV

Route

Expected travel time (days)

Travel time variance (days2)

Expected travel time (days)

Travel time variance (days2)

Expected travel time (days)

Travel time variance (days2)

Expected travel time (days)

Travel time variance (days2)

A-B A-C A-D B-A B-C B-D C-A C-B C-D D-A D-B D-C

1.3 1.2 2.0 1.0 1.5 1.3 1.5 1.5 1.0 2.0 2.0 1.0

0.30 0.40 0.50 0.50 0.30 0.30 0.30 0.30 0.50 0.50 0.50 0.60

1.2 1.0 1.0 1.4 1.5 1.5 1.0 1.5 1.0 1.0 2.0 1.5

0.50 0.40 0.50 0.40 0.50 0.50 0.40 0.50 0.40 0.40 0.60 0.50

1.0 1.5 2.5 1.0 2.0 1.0 1.0 2.0 1.4 2.5 1.0 1.0

0.55 0.50 0.50 0.55 0.65 0.55 0.55 0.65 0.45 0.50 0.55 0.55

1.0 1.5 1.5 1.0 2.0 1.5 1.0 2.0 1.0 2.0 1.5 1.0

0.50 0.50 0.50 0.50 0.65 0.50 0.50 0.65 0.50 0.65 0.50 0.50

Optimization models for rail car fleet management

Table 6.4 Expected travel times and variances for empty freight cars

Table 6.5 Proportion of loaded cars arrived during period n and unit car ownership costs Qt Freight car type I C-A

D-A

A-B

C-B

D-B

A-C

B-C

D-C

A-D

B-D

C-D

0 1 2 3

0 0.5695 0.3516 0.0353

0 0.5267 0.2805 0.0781

0 0.1733 0.4040 0.3175

0 0.5267 0.2805 0.0781

0 0.1733 0.4040 0.3175

0 0.5695 0.3516 0.0353

0 0.5695 0.3516 0.0353

0 0.1733 0.4040 0.3175

0 0.5267 0.2805 0.0781

0 0.1436 0.3821 0.2945

0 0.5267 0.2805 0.0781

0 0.5695 0.3516 0.0353

0 0.5070 0.2836 0.0900

0 0.1554 0.3750 0.3163

0 0.5070 0.2836 0.0900

0 0.1589 0.3808 0.3168

0 0.5070 0.2836 0.0900

0 0.5070 0.2836 0.0900

0 0.6422 0.3507

0 0.1613 0.4130 0.2928

0 0.6422 0.3507

0 0.1613 0.4130 0.2928

0 0.6422 0.3507

0 0.5267 0.2805 0.0781

0 0.5695 0.3516 0.0353

0 0.1444 0.3543 0.3168 15

0 0.5267 0.2805 0.0781

0 0.1436 0.3821 0.2945

0 0.5267 0.2805 0.0781 10

0 0.5267 0.2805 0.0781

II

0 1 2 3

0 0.5362 0.3484 0.0572

0 0.5362 0.3484 0.0572

0 0.1589 0.3808 0.3168

0 0.5362 0.3484 0.0572

0 0.1589 0.3808 0.3168

0 0.5362 0.3484 0.0572 III

0 1 2 3

0 0.6422 0.3507

0 0.5267 0.2805 0.0781

0 0.1512 0.3661 0.3189

0 0.6422 0.3507

0 0.1613 0.4130 0.2928

0 0.6422 0.3507 IV

0 1 2 3 Qt

0 0.5267 0.2805 0.0781

0 0.1436 0.3821 0.2945 10

0 0.1436 0.3821 0.2945

0 0.5695 0.3516 0.0353

0 0.1444 0.3543 0.3168 15

0 0.5695 0.3516 0.0353

209

B-A

Stochastic model for heterogeneous rail freight car fleet management

n

210

Table 6.6 Proportion of empty cars arrived during period n Freight car type

n

B-A

C-A

D-A

A-B

C-B

D-B

A-C

B-C

D-C

A-D

B-D

C-D

0 1 2 3

0 0.5695 0.3516 0.0353

0 0.5457 0.2761 0.0659

0 0.1733 0.4040 0.3175

0 0.5363 0.2784 0.0720

0 0.5457 0.2761 0.0659

0 0.5457 0.2761 0.0659

0 0.6345 0.3513

0 0.5457 0.2761 0.0659

0 0.5362 0.3484 0.0572

0 0.1733 0.4040 0.3175

0 0.5363 0.2784 0.0720

0 0.5362 0.3484 0.0572

0 0.6072 0.3525 0.0085

0 0.6072 0.3525 0.0085

0 0.6345 0.3513

0 0.5070 0.2836 0.090

0 0.1344 0.3364 0.3125

0 0.6072 0.3525 0.0085

0 0.5070 0.2836 0.090

0 0.5070 0.2836 0.090

0 0.5695 0.3516 0.0353

0 0.5070 0.2836 0.090

0 0.6072 0.3525 0.0085

0 0.6345 0.3513

0 0.1583 0.3787 0.3200

0 0.6345 0.3513

0 0.1312 0.3304 0.3108

0 0.6345 0.3513

0 0.5070 0.2836 0.090

0 0.1312 0.3304 0.3108

0 0.6345 0.3513

0 0.1588 0.3787 0.3200

0 0.6345 0.3513

0 0.5095 0.2833 0.0885

0 0.5695 0.3516 0.0353

0 0.1312 0.3304 0.3108

0 0.5695 0.3516 0.0353

0 0.1312 0.3304 0.3108

0 0.5070 0.2836 0.090

0 0.5070 0.2836 0.090

0 0.1312 0.3304 0.3108

0 0.5695 0.3516 0.0353

0 0.5070 0.2836 0.090

0 0.5070 0.2836 0.090

0 0.5695 0.3516 0.0353

II

0 1 2 3

0 0.5169 0.2822 0.0842 III

0 1 2 3

0 0.6345 0.3513 IV

0 1 2 3

0 0.5695 0.3516 0.0353

Optimization models for rail car fleet management

I

Table 6.7 Daily demand between stations on the considered part of the rail network Daily demand (SD)

Daily demand (SD)

1

2

3

4

Route

1

2

3

4

A-B

I II III IV V VI I II III IV V VI I II III IV V VI I II III IV V VI

10 (2) 16 (4) 5 (1) 4 (2) 9 (1) 10 (2) 8 (2) 13 (3) 8 (3) 7 (3) 10 (2) 9 (1) 12 (3) 17 (4) 6 (2) 5 (1) 10 (2) 10 (2) 3 (1) 7 (2) 11 (1) 6 (1) 7 (1) 9 (2)

5 (1) 9 (1) 11 (3) 5 (1) 8 (2) 11 (2) 7 (2) 11 (2) 10 (2) 3 (1) 12 (3) 11 (2) 9 (1) 15 (3) 7 (1) 6 (2) 11 (2) 13 (2) 6 (1) 11 (3) 12 (2) 3 (1) 5 (1) 6 (1)

4 (1) 10 (3) 2 (1) 3 (1) 12 (2) 13 (3) 5 (1) 6 (2) 3 (1) 5 (2) 13 (4) 12 (3) 7 (1) 3 (1) 1 (1) 3 (1) 7 (1) 15 (2) 11 (2) 8 (2) 3 (1) 6 (1) 6 (1) 10 (2)

7 (2) 12 (3) 5 (1) 7 (2) 11 (1) 12 (2) 9 (2) 5 (1) 9 (2) 6 (1) 11 (2) 10 (2) 11 (3) 7 (2) 13 (3) 9 (2) 9 (2) 9 (1) 12 (2) 7 (1) 9 (2) 9 (2) 11 (2) 9 (1)

C-A

7 (1) 17 (3) 12 (3) 6 (2) 11 (1) 13 (2) 8 (2) 9 (2) 18 (3) 9 (3) 10 (2) 8 (1) 11 (3) 5 (1) 12 (2) 4 (1) 10 (2) 11 (2) 4 (1) 9 (2) 5 (1) 7 (2) 6 (1) 7 (1)

13 (2) 5 (1) 14 (3) 6 (2) 9 (1) 10 (2) 12 (2) 7 (2) 19 (3) 11 (3) 11 (2) 9 (1) 6 (1) 9 (2) 11 (2) 18 (4) 12 (2) 13 (2) 6 (1) 11 (3) 19 (3) 5 (1) 10 (1) 12 (2)

9 (1) 6 (2) 19 (4) 4 (1) 7 (1) 6 (1) 13 (2) 8 (2) 7 (1) 8 (2) 13 (1) 6 (1) 21 (4) 15 (3) 12 (3) 5 (1) 9 (1) 8 (1) 3 (1) 15 (3) 17 (3) 13 (2) 11 (2) 9 (1)

18 (4) 11 (3) 10 (2) 16 (4) 8 (1) 9 (1) 22 (5) 6 (2) 12 (3) 11 (2) 12 (2) 7 (1) 29 (5) 9 (2) 8 (3) 6 (1) 9 (1) 6 (1) 9 (2) 14 (3) 12 (3) 12 (3) 12 (2) 9 (1)

A-C

A-D

B-A

C-B

C-D

D-A

Continued

211

Freight car type

Stochastic model for heterogeneous rail freight car fleet management

Route

212

Daily demand (SD)

Daily demand (SD)

Route

Freight car type

1

2

3

4

Route

1

2

3

4

B-C

I II III IV V VI I II III IV V VI

6 (2) 9 (2) 13 (2) 8 (2) 8 (1) 7 (2) 2 (1) 12 (3) 8 (2) 4 (1) 10 (1) 11 (2)

9 (2) 15 (3) 6 (1) 9 (2) 9 (1) 8 (1) 13 (3) 6 (1) 15 (3) 12 (3) 13 (2) 9 (1)

15 (3) 3 (1) 5 (1) 7 (2) 10 (2) 9 (1) 8 (2) 5 (1) 4 (1) 9 (2) 11 (2) 10 (1)

14 (4) 9 (2) 13 (4) 11 (3) 10 (2) 11 (2) 11 (2) 13 (3) 12 (3) 15 (4) 12 (2) 11 (2)

D-B

7 (2) 13 (3) 9 (2) 3 (1) 5 (1) 4 (1) 5 (1) 5 (1) 15 (3) 10 (3) 6 (1) 9 (1)

9 (3) 5 (2) 27 (4) 7 (2) 10 (1) 12 (2) 5 (2) 4 (1) 5 (1) 3 (1) 8 (1) 9 (1)

3 (1) 6 (1) 5 (1) 18 (3) 11 (1) 10 (1) 4 (1) 7 (1) 11 (2) 5 (1) 11 (2) 5 (1)

11 (3) 11 (3) 14 (4) 19 (5) 13 (2) 12 (2) 8 (2) 16 (4) 11 (4) 9 (3) 14 (2) 10 (2)

B-D

D-C

Optimization models for rail car fleet management

Table 6.7 Daily demand between stations on the considered part of the rail network—cont’d

Station A

B

C

D

Freight car type

hit

S0(SD)

hit

S0(SD)

hit

S0(SD)

hit

S0(SD)

Lt(m)

I II III IV

20 15 15 20

30 39 38 36

15 15 15 15

28 30 32 35

20 15 15 20

26 37 33 32

15 20 15 20

38 28 24 23

14 15 15 16

CAPi(m)

(3.00) (4.00) (2.00) (4.00)

2000

(4.00) (3.00) (2.00) (3.00)

1200

(3.00) (2.00) (3.00) (2.00)

1800

(5.00) (3.00) (4.00) (4.00)

1200

Stochastic model for heterogeneous rail freight car fleet management

Table 6.8 Unit holding cost (hit), initial number of cars in stations with standard disturbance (S0(SD)), maximum station’s capacities (CAPi(m)) and lengths of freight cars (Lt(m))

213

214

Optimization models for rail car fleet management

Table 6.9 Unmet demand from the preceding period Route

Freight car type

A-B

A-C

A-D

B-A

B-C

B-D

I II III IV V VI

4 1 2 3 2 5

3 5 3 4 4 4

9 10 4 7 3 2

3 2 11 5 5 5

5 3 3 4 7 6

7 11 6 8 6 1

C-A

C-B

C-D

D-A

D-B

D-C

12 13 2 6 1 1

4 5 3 7 2 3

5 4 5 9 1 2

3 3 4 11 3 4

10 4 6 6 4 3

2 5 2 5 5 3

I II III IV V VI

Table 6.10 Number of loaded cars dispatched before the beginning of the cycle Route

A-B

A-C

A-D

B-A

B-C

B-D

Freight car type

One period

Two periods

I II III IV I II III IV I II III IV I II III IV I II III IV I II III IV

12 12 14 14 18 21 19 20 21 12 14 14 10 11 13 14 11 12 22 22 23 14 11 13

16 21 8 9 15 9 14 11 10 11 24 26 10 5 15 11 13 15 6 19 13 9 13 14

Route

C-A

C-B

C-D

D-A

D-B

D-C

One period

Two periods

8 19 10 25 14 14 15 25 14 15 5 25 14 13 23 13 25 15 6 11 13 5 8 13

9 13 7 16 21 19 9 21 9 7 19 11 17 14 16 7 13 17 5 23 11 16 12 26

Stochastic model for heterogeneous rail freight car fleet management

215

Table 6.11 Number of empty cars dispatched before the beginning of the cycle Route

A-B

A-C

A-D

B-A

B-C

B-D

Freight car type

One period

Two periods

I II III IV I II III IV I II III IV I II III IV I II III IV I II III IV

3 4 2 6 3 3 3 5 3 2 3 3 3 6 3 3 4 4 4 4 2 6 9 9

4 4 5 4 3 2 5 6 5 4 5 6 2 5 5 9 5 4 4 7 6 5 7 6

Route

C-A

C-B

C-D

D-A

D-B

D-C

One period

Two periods

4 4 6 2 5 5 6 3 4 9 4 2 2 3 11 4 2 3 2 6 6 5 2 6

3 7 9 15 3 5 7 4 4 3 7 11 4 5 3 17 4 6 5 9 2 6 4 7

The best-fitting model was found based on the assessment of the alternative ARIMA models. For the first rail car type in the first station, the most appropriate model was found with the Box-Jenkins approach and AIC criteria. This approach involved (1) the selection of the candidate model set, (2) the estimation of the model and determination of AIC, and (3) a diagnostic check. According to the Box-Jenkins approach (Section 6.5), selection of the candidate models first starts with the evaluation of sample estimates of the ACF and the PACF, in order to determine three major orders of the ARIMA models. As it can be noticed from Fig. 6.4, there is a gradual decrease in autocorrelation values which indicates a long-term trend. Also, there is one very large spike on the PACF plot. This reflects the need for

216

Optimization models for rail car fleet management

Table 6.12 Substitution cost rijtt0 Route

t 5 1, t0 5 5

t 5 3, t0 5 5

t 5 2, t0 5 6

t 5 4, t0 5 6

A-B A-C A-D B-A B-C B-D C-A C-B C-D D-A D-B D-C

10 20 15 15 20 25 15 25 10 10 20 15

15 25 20 20 30 35 20 35 15 15 30 10

10 15 30 15 20 10 10 25 15 25 20 10

15 20 35 20 25 15 15 20 20 35 25 15

including a first-order difference term in the model structure (d ¼ 1). Therefore, the basic structure of the alternative ARIMA models has the form ARIMA(p, 1, q). Now, by evaluating the different AR and MA orders (p and q), ARIMA(2, 1, 14) was selected as the most appropriate model with the lowest AIC (AIC ¼5.3562) with R-squared equal to 92.6%. Model equation has the following form:  1 + 0:577ω + 0:284ω2 ð1  ωÞyl ¼ð1  0:798ω  0:264ω2 + 0:087ω3  0:002ω4 + 0:136ω5 + 0:187ω6 + 0:008ω7  0:108ω8 + 0:088ω9 + 0:006ω10  0:013ω11  0:679ω12  0:401ω13  0:145ω14 Þεl (6.106)

Diagnostic checking proves the stationarity and invertibility of the selected model without redundant parameters. The absence of autocorrelation between residuals in different lags (Fig. 6.5) indicated that the residuals were white noise (Ljung-Box Q ¼ 6.295, P-value >0.05). Forecasting was performed by incorporating the identified ARIMA model in a state-space form. As explicit equation of the identified ARIMA(2,1,14) model was obtained as follows. First, the state vector ρl of ARIMA(2,1,14) model was defined as ρl ¼ ðyl1 , y∗l ,  0:577y∗l + 0:798εl + 0:264εl1 + ⋯ + 0:145εl13 , 0:264εl  0:087εl1 + ⋯ + 0:145εl12 , …,0:401εl + 0:145εl1 , 0:145εl Þ (6.107)

γ it

γ Pit

γ it

Station 1

Freight car type

1 2 3 4

ηit

γ 11 ¼ 71.70 γ 12 ¼ 94.84 γ 13 ¼ 70.03 γ 14 ¼ 62.72

η11 ¼ 50.69 η12 ¼ 64.11 η13 ¼ 38.97 η14 ¼ 40.24

γ 31 ¼ 75.05 γ 32 ¼ 69.79 γ 33 ¼ 89.03 γ 43 ¼ 70.44

η31 ¼ 63.49 η32 ¼ 50.97 η33 ¼ 66.92 η43 ¼ 49.95

γ Pit

Station 2

γ 311 ¼ 49.70 γ 312 ¼ 63.01 γ 313 ¼ 38.57 γ 314 ¼ 39.55

γ 21 ¼ 74.39 γ 22 ¼ 69.41 γ 23 ¼ 67.99 γ 24 ¼ 68.63

Station 3

1 2 3 4

ηit

η21 ¼ 50.80 η22 ¼ 44.99 η23 ¼ 48.97 η24 ¼ 44.05

γ 321 ¼ 49.92 γ 322 ¼ 44.34 γ 323 ¼ 47.99 γ 324 ¼ 43.35

Station 4

γ 331 ¼ 62.06 γ 332 ¼ 50.01 γ 333 ¼ 65.66 γ 334 ¼ 49.02

γ 41 ¼ 67.69 γ 42 ¼ 64.16 γ 43 ¼ 84.32 γ 44 ¼ 78.75

η41 ¼ 39.60 η42 ¼ 47.97 η43 ¼ 59.14 η44 ¼ 46.00

γ 341 ¼ 39.16 γ 342 ¼ 46.97 γ 343 ¼ 58.20 γ 344 ¼ 45.49

Stochastic model for heterogeneous rail freight car fleet management

Table 6.13 Rail freight car inventory level and ordering quantity

217

218

Optimization models for rail car fleet management

Number of rail freight cars on state

100.00

80.00

60.00

40.00

100 97 94 91 88 85 82 79 76 73 70 67 64 61 58 55 52 49 46 43 40 37 34 31 28 25 22 19 16 13 10 7 4 1

Time (days)

Fig. 6.3 Time series plot: daily state of rail freight cars of type 1 in station 1.

Fig. 6.4 Sample autocorrelation function (ACF) and partial autocorrelation function (PACF) of the time series data.

With the corresponding disturbance vector I l ηl ¼ðO11 , εl + 1 , 0:798εl + 1 ,0:264εl + 1 ,  0:087εl + 1 , 0:002εl + 1 ,  0:136εl + 1 ,  0:187εl + 1 ,  0:008εl + 1 ,0:108εl + 1 ,  0:089εl + 1 ,  0:006εl + 1 ,0:013εl + 1 , 0:679εl + 1 ,0:401εl + 1 ,0:145εl + 1 Þ0 (6.108)

Stochastic model for heterogeneous rail freight car fleet management

Lag

Residual ACF

219

Residual PACF

20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1

–1.0

–0.5

0.0

0.5

1.0

–1.0

–0.5

0.0

0.5

1.0

Residual

Fig. 6.5 Autocorrelation function (ACF) and partial autocorrelation function (PACF) of residuals.

The transition matrix Y l is, therefore, 16  16 and Y l and Sl are given by

ð6:109Þ Y d ¼ ½I14 , Y c ¼ ½0:5769  0:2842 O141 0 , Y e ¼ ½O115 , Sl ¼ ½1 1 O114  Then the variance-covariance matrix of the state disturbances can be defined as   O11 O115 0 (6.110) I lI l ¼ 0 O151 I ∗l I ∗l 0

The 15  15 stationary part of this variance matrix is given by I ∗l I ∗l , Variance matrices M l and V l have dimensions (11), values 9.54 and 0, respectively. The mean of the initial state vector is given by a ¼ E(α1) ¼ O161 and the corresponding variance matrix is given by   kI 1 O115 with k ! ∞ (6.111) Π¼ O151 Π ∗1515 Matrix Π ∗1515 represents the unconditional variance matrix for the stationary part of the state vector.

2

15:22 6 10:67 6 6 3:46 6 6 21:37 6 6 0:22 6 6 21:06 6 6 22:15 6 ∗ Π 1515 ¼ 6 6 21:07 6 0:58 6 6 21:21 6 6 21:12 6 6 0:90 6 6 7:12 6 4 4:13 1:38

10:67 12:92 5:58 0:14 20:41 20:40 21:01 20:99 20:66 21:15 21:01 20:54 4:89 3:02 1:10

3:46 5:58 7:57 3:18 0:87 20:73 0:01 0:11 20:76 21:82 20:79 20:71 1:38 0:89 0:36

21:37 0:14 3:18 6:91 3:4 0:86 20:39 0:48 0:13 21:03 21:60 20:77 20:74 20:33 20:12

0:22 20:41 0:87 3:4 6:84 3:40 0:75 20:55 0:47 0:22 21:11 21:60 20:76 20:18 0:03

21:06 20:40 20:73 0:86 3:40 6:84 3:40 0:75 20:54 0:47 0:22 21:11 21:60 20:78 20:18

22:15 21:01 0:01 20:39 0:75 3:40 6:66 3:16 0:74 20:40 0:36 0:21 21:09 20:72 20:25

21:07 20:99 0:11 0:48 20:54 0:75 3:16 6:33 3:14 0:94 20:56 0:34 0:29 0:11 20:01

0:58 20:66 20:76 0:13 0:47 20:54 0:74 3:14 6:32 3:15 0:93 20:56 0:35 0:28 0:14

21:21 21:14 21:82 21:03 0:22 0:47 20:40 0:94 3:15 6:21 3:24 0:94 20:58 20:34 20:12

21:12 21:01 20:79 21:60 21:11 0:22 0:36 20:56 0:93 3:24 6:14 3:24 0:95 20:00 20:00

0:90 20:54 20:71 20:77 21:60 21:11 0:21 0:34 20:56 0:94 3:24 6:14 3:24 0:99 0:01

7:12 4:89 1:38 20:74 20:76 21:60 21:09 0:23 0:35 20:58 0:95 3:24 6:14 3:15 0:94

4:13 3:02 0:89 20:33 20:18 20:78 20:72 0:11 0:28 20:34 20:00 0:99 3:15 1:73 0:55

3 1:38 1:10 7 7 0:36 7 7 20:12 7 7 0:00 7 7 20:18 7 7 20:25 7 7 20:01 7 7 0:14 7 7 20:12 7 7 20:00 7 7 0:01 7 7 0:94 7 7 0:55 5 0:20

Stochastic model for heterogeneous rail freight car fleet management

221

The obtained ARIMA-Kalman hybrid forecasting model was incorporated in the Matlab (MathWorks Inc., 2012). For determining the most suitable ARIMA model IBM SPSS Statistics 19 (trial version) was used. Conversion into the state-space model and forecasting were performed with the SSM Toolbox (Peng and Aston, 2006) and SSMMATLAB (Gomez, 2012). The last 15 observations were used for the assessment of the model ARIMA-Kalman model performances. This was performed by comparison of the hybrid ARIMA-Kalman model with the ARIMA model, and evaluating the daily forecast errors based on the following accuracy measures: root mean square error (RMSE), mean error (ME), and mean absolute percentage error (MAPE). According to the results of these accuracy measures, the ARIMAKalman model reports better performances in comparison to ARIMA. It results in a 50.6 reduction in RMSE, 57.1 reduction in MAE, and 55.5% reduction in MAPE. Fig. 6.6 presents the fitting and forecasting results. The same procedure was repeated for other freight car types and all stations. Model equations and short-term forecasts are summarized in Table 6.14. The proposed solution procedure was applied for determining the optimal control law for rail the freight car fleet size and allocation system. The computer algorithm for performing all computations was coded in the Matlab software (MathWorks Inc., 2012). Experiments were carried out on the Intel Core i3-2350M (2.30 GHz) personal computer. Figs. 6.7–6.10 summarize the obtained results.

Fig. 6.6 ARIMA-Kalman: fitting and forecasting results.

Table 6.14 Models and forecasts of cars at stations (SD) per day Freight car type

Time period Forecasting model

1

2

3

61.76 (4.64)

64.04 (5.28)

65.36 (7.96)

55.06 (8.9) 46.72 (9.81)

55.18 (6.2) 46.19 (10.25)

52.65 (8.95) 47.23 (10.31)

40.26 (4.35)

40.96 (4.62)

41.34 (5.96)

Station A

I

II III IV

ð1 + 0:577B + 0:284B2 Þð1  BÞyn ¼ ð1  0:798B  0:264B2 + 0:087B3  0:002B4 + 0:136B5 + 0:187B6 + 0:008B7  0:108B8 + 0:088B9 + 0:006B10  0:013B11  0:679B12  0:401B13  0:145B14 Þεn (1  0.100B  0.517B2)(1  B)yn ¼ (1  329B  0.491B2 + 0.169B3)εn ð1 + 0:035B + 0:202B2 + 0:845B3 + 0:154B4 Þð1  BÞyn ¼ ð1  0:193B + 0:263B2 + 0:769B3 Þεn ð1 + 0:598BÞð1  BÞyn ¼ ð1 + 0:078B + 0:078B2 + 0:022B3  0:091B4 + 0:160B5 + 0:097B6  0:0053B7  0:219B8 0:140B9 + 0:216B10  0:056B11 + 0:555B12 Þεn Station B

I II III

IV

(1  1.079B + 0.355B )(1  B)yn ¼ (1  1.000B)εn (1  0.739B + 0.186B2)(1  B)yn ¼ (1  0.542B  0.064B2)εn ð1  BÞyn ¼ ð1 + 0:164B  0:180B2 +0:118B3  0:177B4 + 0:151B5  0:388B6 + 0:005B7 + 0:039B8 0:444B9  0:144B10 + 0:422B11 + 0:522B12 Þεn (1 + 0.935B + 0.126B2)(1  B)yn ¼ (1 + 0.997B)εn 2

39.77 (3.13) 41.66 (5.91) 46.38 (7.19)

40.62 (2.41) 41.58 (6.25) 47.10 (8.21)

40.90 (4.72) 41.41 (6.11) 49.76 (8.50)

61.37 (9.12)

60.65 (8.20)

61.28 (7.91)

Station C

I II III IV

(1  0.597B + 0.039B + 0.214B  0.059B )(1  B)yn ¼ (1  0.807B)εn ð1 + 0:843B  0:227B2  0:378B3 + 0:288B4 Þð1  BÞyn ¼ ð1 + 0:684B  0:698B2  0:968B3  0:025B4 Þεn ð1 + 0:918B + 0:107B2 + 0:161B3 + 0:568B4 Þð1  BÞyn ¼ ð1 + 0:674B  0:518B2  0:643B3  0:001B4  0:204B5 Þεn ð1 + 0:184B  0:092B2 + 0:130B3 + 0:324B4  0:423B5 Þð1  BÞyn ¼ ð1  0:212B  0:530B2  0:008B3 + 0:065B4 0:820B5 + 0:299B6 + 0:734B7 Þεn 2

3

4

42.42 (5.21) 55.66 (4.90)

45.13 (6.35) 52.41 (3.26)

46.29 (5.01) 53.12 (4.05)

52.20 (5.01)

52.11 (3.92)

57.42 (4.65)

49.58 (7.29)

48.02 (7.35)

45.31 (7.85)

47.43 (3.23)

48.02 (3.56)

49.93 (3.12)

59.43 (4.51)

59.16 (6.20)

57.95 (5.36)

50.57 (6.36)

51.46 (4.89)

52.92 (7.20)

60.00 (5.25)

60.46 (3.92)

61.59 (6.00)

Station D

I

II

III IV

ð1 + 0:270B  0:410B2 + 0:415B3 + 0:478B4  0:319B5 Þð1  BÞyn ¼ ð1  0:059B  0:597B2 + 0:486B3 + 0:212B4 0:163B5  0:023B6  0:546B7 Þεn ð1  0:931B + 0:688B2  0:551B3 Þð1  BÞyn ¼ ð1  0:560B  0:321B2  0:356B3 0:076B4 + 0:471B5  0:158B6 Þεn ð1  0:091B + 0:657B2 + 0:538B3 Þð1  BÞyn ¼ ð1  0:078B + 0:486B2 + 0:646B3  0:186B4 Þεn ð1  1:123B + 0:530B2 + 0:276B3  0:326B4 + 0:159B5 Þð1  BÞyn ¼ ð1  0:945B  0:176B2 + 1:175B3 0:372B4  0:746B5 + 0:741B6 Þεn

224 Optimization models for rail car fleet management

Fig. 6.7 Control actions (loaded cars).

E

E

E E E

E E

E E

Stochastic model for heterogeneous rail freight car fleet management

Control action

E

E

E

225

Fig. 6.8 Control actions (empty cars).

E

E E

E

E E E

E E

E E E

E E

E

E E

E E

E E

E E

E

226

Unmet demand

V

V

V

Fig. 6.9 Unmet demands.

V

V

V

V

V

V

V

V

V

V

V

V V

V

V

V

V

V

V

V

V

V

V

V

V

V V

V V

V V

Optimization models for rail car fleet management

V

V

Number of cars

Stochastic model for heterogeneous rail freight car fleet management

Fig. 6.10 Number of cars in stations by periods.

227

228

Optimization models for rail car fleet management

The model provides rail network information related to the number of cars (Sit(n)) in all stations by type and period, unmet orders (Uijt(n)) and empty and loaded freight car movements for every period, and between all pairs of stations (Xijtt’(n), Yijtt’(n)). Now, it is possible to determine the fleet size required for the proper functioning of the rail freight car system. The fleet size actually represents the average number of rail cars available within the planning period. To get this result, the following expression is applied: FS ¼ where Að1Þ ¼

X XX

P 1X AðnÞ P n¼1

(6.112)

ðEijtt0 ð1Þαijt ð0, 1Þ + Fijtt0 ð1Þθijt ð0, 1Þ + Eijtt0 ð2Þαijt ð0, 3Þ

i, j2N t2T t0 2T t

+ Fijtt0 ð2Þθijt ð0, 3Þ + Eijtt0 ð3Þαijt ð0, 2Þ + Fijtt0 ð3Þθijt ð0, 2ÞÞ +

XX

Sit ð1Þ

i2N t2T

Að2Þ ¼

X XX

(6.113) ðEijtt0 ð2Þαijt ð0, 1Þ + Fijtt0 ð2Þθijt ð0, 1Þ + Eijtt0 ð3Þαijt ð0, 3Þ

i, j2N t2T t0 2T t

+ Fijtt0 ð3Þθijt ð0, 3Þ + Eijtt0 ð1Þαijt ð0, 2Þ + Fijtt0 ð1Þθijt ð0, 2ÞÞ +

XX

Sit ð2Þ

i2N t2T

Að3Þ ¼

X XX

(6.114) ðEijtt0 ð3Þαijt ð0, 1Þ + Fijtt0 ð3Þθijt ð0, 1Þ + Eijtt0 ð1Þαijt ð0, 3Þ

i, j2N t2T t0 2T t

+ Fijtt0 ð1Þθijt ð0, 3Þ + Eijtt0 ð2Þαijt ð0, 2Þ + Fijtt0 ð2Þθijt ð0, 2ÞÞ +

XX

Sit ð3Þ

i2N t2T

(6.115) The optimal fleet size (FS) is 6120 freight cars. This fleet generates the total cost of J ¼ 7.46  106 monetary units. Table 6.15 summarizes the results for other tests that involve an increased number of railway stations. The CPU time for the performed tests is insignificant.

Stochastic model for heterogeneous rail freight car fleet management

229

Table 6.15 Alternate example problems MPC

Problem

Number of stations

Number of periods P

Rail freight car type

FS

I II III IV I II III IV I II III IV I II III IV I II III IV

272.6 298.9 181.2 210.2 1555.1 1599.1 1695.5 1270.6 1925.3 1258.6 1100.9 1736.3 2929.2 2027.1 1645.6 1472.3 3834.2 3439.9 3952.3 3179.1

1

2

4

2

4

4

3

5

4

4

10

4

5

20

4

J (106 m.u.)

CPU time (s)

2.49

00:00:01

7.46

00:00:02

8.95

00:00:04

18.3

00:00:12

38.2

00:00:21